[Vision] Depth estimation from Blur estimation

Introduction When images are captured with a small depth of field, objects that are away from the focal plane are out of focus and are perceived as blurry. This effect usually occurs at relatively large apertures or when the focal plane is close to the lens. An image that includes objects in focus and out of focus could therefore be segmented in terms of depth. If we are able to measure the size of this blur, we can also produce a three dimensional depth map of the same image. In this brief overview, blur estimation will be performed using methods proposed by Elder and Zucker (1998) and Hu and Haan (2006).

[caption id="attachment_2130" align="aligncenter" width="300" caption="Test image (35mm f1.4)"][/caption]

Small videos showing the resulting depth map (nothing fancy)

The processing speed was about 10fps for the first algorithm, and 1fps for the second. So the first can be done realtime with some optimization. http://www.youtube.com/watch?v=68i5wgVC0OY http://www.youtube.com/watch?v=w4kySlumyUs

Hu and Haan algorithm

  1. %Tim Zaman, TU Delft ME, 01-01-2012
  2. %This code is licensed under a
  3. %Creative Commons Attribution 3.0 Unported License.
  4. clc;clear all;close all
  5. im=imread('imnikon.jpg');
  6. im=im(:,:,1);
  7. sigmaA=8;
  8. sigmaB=10;
  9. sigmaMax=max([sigmaA sigmaB]);
  10. fsz=[sigmaMax, sigmaMax];
  11. kernelA=fspecial('gaussian', fsz, sigmaA);
  12. kernelB=fspecial('gaussian', fsz, sigmaB);
  13. %imA = filter2(kernelA, im);
  14. %imB = filter2(kernelB, im);
  15. imA = imfilter(im, kernelA, 'symmetric', 'conv');
  16. imB = imfilter(im, kernelB, 'symmetric', 'conv');
  17. R1=single(im)-single(imA);
  18. %R1=((single(imA)./single(im))-1).*single(imB);
  19. R2=single(imA)-single(imB);
  20. %R1=abs(R1);
  21. %R2=abs(R2);
  22. %R2(R2<1)=1;
  23. R=R1./R2;
  24. R(isnan(R))=0;
  25. R(isinf(R))=0;
  26. Rf=maxfilt2(R,8);
  27. %Rf=max(max(Rf))-Rf;
  28. blurmap=(sigmaA*sigmaB)./((sigmaB-sigmaA)*Rf+sigmaB);
  29. figure;
  30. subplot(3,3,1), imshow(im), title('original')
  31. subplot(3,3,2), imshow(uint8(imA)), title('imA')
  32. subplot(3,3,3), imshow(uint8(imB)), title('imB')
  33. subplot(3,3,4), imagesc(R1), title('R1')
  34. set(gca,'xticklabel',[],'yticklabel',[])
  35. set(gca,'xtick',[],'ytick',[]),axis image
  36. subplot(3,3,5), imagesc(R2), title('R2')
  37. set(gca,'xticklabel',[],'yticklabel',[])
  38. set(gca,'xtick',[],'ytick',[]),axis image
  39. subplot(3,3,6), imagesc(R), title('R')
  40. set(gca,'xticklabel',[],'yticklabel',[])
  41. set(gca,'xtick',[],'ytick',[]),axis image
  42. subplot(3,3,7), imagesc(blurmap), title('blurmap')
  43. set(gca,'xticklabel',[],'yticklabel',[])
  44. set(gca,'xtick',[],'ytick',[]),axis image
  45. subplot(3,3,8), imagesc(blurmap>1 & blurmap<4), title('sigma [2:4]')
  46. set(gca,'xticklabel',[],'yticklabel',[])
  47. set(gca,'xtick',[],'ytick',[]),axis image
  48. subplot(3,3,9), imagesc(blurmap>0 & blurmap<0.4), title('sigma [0:0.5]')
  49. set(gca,'xticklabel',[],'yticklabel',[])
  50. set(gca,'xtick',[],'ytick',[]),axis image

Elder and Zucker algorithm

  1. %Tim Zaman, TU Delft ME, 01-01-2012
  2. %This code is licensed under a
  3. %Creative Commons Attribution 3.0 Unported License.
  4. clc;close all;clear all
  5. %Read image
  6. im = imread('imnikon.jpg');
  7. im = im(:,:,1);
  8. %Some variables
  9. sn = 3; %image noise [sigma]
  10. alphap = 2e-7; %2e-7
  11. %Apply white noise sn to image (optional)
  12. %im = imnoise(im,'gaussian',0,(sn/255)^2);
  13. imd=double(im);
  14. sigma1 = [8 4 2 1 0.5]; %Scale space octaves for 1st deriv
  15. sigma2 = sigma1;
  16. %Preallocate some matrices
  17. mrsmap1=ones(size(im))*4;
  18. mrsmap2=ones(size(im))*4;
  19. edge2nd=zeros(size(im));
  20. edge3rd=zeros(size(im));
  21. %Make steerable gaussian 1st order derivative filters
  22. for i=1:length(sigma1)
  23.     %First derivative part
  24.     %Create normal gaussian filter kernel
  25.     h = fspecial('gaussian', round(10*sigma1(i)), sigma1(i));
  26.     %Create gaussian 1st order gradient in X and Y direction
  27.     [g1x,g1y] = gradient(h);
  28.     %Apply the filters on the image
  29.     r1x=imfilter(imd,g1x,'replicate');
  30.     r1y=imfilter(imd,g1y,'replicate');
  31.     thetam=atan2(r1y, r1x);
  32.     r1theta=cos(thetam).*r1x + sin(thetam).*r1y;
  33.     r1theta_or=r1theta;
  34.     %Compute critical value
  35.     g1theta2=1/(2*sqrt(2*pi)*sigma1(i)^2);
  36.     s1=g1theta2*sn;
  37.     c1=s1*sqrt(-2*log(alphap));
  38.     %Make the binary 1st deriv reliable map
  39.     reliable1=(abs(r1theta) >= c1);
  40.     %Put into the Minimum Reliable Scale map (1st)
  41.     mrsmap1(reliable1)=sigma1(i);
  42.     %Apply threshold on r1theta and thetam
  43.     thetam=thetam.*reliable1;
  44.     r1theta=r1theta.*reliable1;
  45.     %Second derivative part
  46.     %Make 2nd order gaussian filter kernels
  47. %    [g2x,g2xy] = gradient(gradient(h));
  48. %    g2y=rot90(g2x);
  49.     %Apply 2nd order gaussian filter kernels to image
  50. %    r2x=imfilter(imd,g2x,'replicate');
  51. %    r2y=imfilter(imd,g2y,'replicate');
  52. %    r2xy=imfilter(imd,g2xy,'replicate');
  53.     [r2x,r2y] = gradient(r1theta_or);
  54.     %Compute 2nd derivative in maginitude direction
  55.     r2theta=cos(thetam).*r2x + sin(thetam).*r2y;
  56.     r2theta_or=r2theta;
  57.     %r2theta =(cos(thetam).^2 .* r2x)+(sin(thetam).^2 .* r2y)...
  58.     %    -(2.*cos(thetam).*sin(thetam).*r2xy);
  59.     g22=1/(4*sqrt(pi/3)*sigma2(i)^3);
  60.     s2=g22*sn;
  61.     c2=sqrt(2)*s2*erfinv(1-alphap); %paper says (alpha)???
  62.     %Make the binary 2nd deriv reliable map
  63.     reliable2=(abs(r2theta) >= c2);
  64.     %Put into the Minimum Reliable Scale map (2nd)
  65.     mrsmap2(reliable2)=sigma2(i);
  66.     %Apply threshold on r2theta
  67.     r2theta=r2theta.*reliable2;
  68.     %Now remove the new lower sigma part for filling in
  69.     edge2nd=edge2nd.*(1-reliable1.*reliable2);
  70.     %Locate the edges at zero crossings in
  71.     % binary image (discrete for simplicity)
  72.     edge2nd=edge2nd+edge(r2theta>0,'canny')*sigma1(i);
  73.     %LOOK FOR ZERO CROSSINGS in [3x3] window
  74. %     for x=1:(size(im,2)-2);
  75. %         for y=1:(size(im,1)-2);
  76. %             %[3x3] window
  77. %             imw=r2theta(y:y+2,x:x+2);
  78. %             %Check for sign changes within window
  79. %             sp=sum(sum(imw>=0));
  80. %             sm=sum(sum(imw<=0));
  81. %             %Apply to edgemap
  82. %             if (sp~=9 && sm~=9)
  83. %                 edge2nd(y,x)=sigma1(i);
  84. %             end
  85. %         end
  86. %     end
  87.     %Now remove the new lower sigma part for filling in
  88.     edge3rd=edge3rd.*(1-reliable1.*reliable2);
  89.     %Make the third derivative
  90.     [r3x,r3y]=gradient(r2theta_or);
  91.     r3theta=cos(thetam).*r3x + sin(thetam).*r3y;
  92.     edge3rd=edge3rd+edge(r3theta>0,'canny')*sigma1(i);
  93.     %Now ignore the edges at faulty locations
  94.     %edge3rd=edge3rd.*imerode(reliable1,ones(10));
  95.     %Now, for every point in edge2nd, run along the direction
  96.     % of the gradient (thetam) and find the first occurancy of
  97.     % a zero in the third derivative map.
  98.     %for ii=1:sum(sum(edge2nd>0)) %For all edgepoints..
  99. %     for y=2:255;
  100. %         x=r2theta(y,:);
  101. %         yy=1:1/5:length(x);
  102. %         %Interpolate the signal (simple solution)
  103. %         p = spline([1:length(x)],x,yy);
  104. %         [~,loc]=findpeaks(abs(p),'MINPEAKDISTANCE',50,'MINPEAKHEIGHT',max(p)*0.95);
  105. %         d=max(abs(diff(loc)/5));
  106. %         sigmab(y)=sqrt(((d/2)^2)-sigma2(i)^2);
  107. %         %plot(x)
  108. %         %pause(.01)
  109. %     end
  110.     %imagesc(r2theta)
  111.     %colormap gray
  112. %     figure
  113. %     plot(sigmab,'-b')
  114. %     hold on
  115. %     plot([1:26.6/length(x):26.6],'--r')
  116. %     hold off, grid on
  117. %     legend('estimated','actual')
  118. %     ylabel('Blur scale (px)')
  119. %     xlabel('Arclength(px)')
  120. %     title(['sigma=' num2str(sn)])
  121.     %pause(2)
  122.     %figure
  123.     %imagesc(r1theta)
  124.     %colormap gray
  125.     %figure
  126.     %plot(r2theta(:,550))
  127.     %pause(1)
  128. end
  129. figure
  130. subplot(2,3,1), imshow(im)
  131. title('original')
  132. subplot(2,3,2), imagesc(sqrt(mrsmap1))
  133. set(gca,'xticklabel',[],'yticklabel',[])
  134. set(gca,'xtick',[],'ytick',[]),axis image
  135. %colormap hot
  136. title('mrs map 1st deriv')
  137. subplot(2,3,3), imagesc(sqrt(mrsmap2))
  138. set(gca,'xticklabel',[],'yticklabel',[])
  139. set(gca,'xtick',[],'ytick',[]),axis image
  140. %colormap hot
  141. title('mrs map 2nd deriv')
  142. subplot(2,3,6), imagesc(mrsmap2==0.5)
  143. title('sigma=0.5')
  144. set(gca,'xticklabel',[],'yticklabel',[])
  145. set(gca,'xtick',[],'ytick',[]),axis image
  146. subplot(2,3,5), imagesc(mrsmap2==1)
  147. title('sigma=1')
  148. set(gca,'xticklabel',[],'yticklabel',[])
  149. set(gca,'xtick',[],'ytick',[]),axis image
  150. subplot(2,3,4), imagesc(mrsmap2==8)
  151. title('sigma=8')
  152. set(gca,'xticklabel',[],'yticklabel',[])
  153. set(gca,'xtick',[],'ytick',[]),axis image
  154. %subplot(1,3,1), [f1]=imagesc(r2theta,[-1 1]);
  155. %axis image, colormap gray
  156. %title('2nd derivative map')
  157. %set(gca,'xticklabel',[],'yticklabel',[])
  158. %set(gca,'xtick',[],'ytick',[])
  159. %subplot(1,3,2), [f2]=imagesc(r3theta,[-0.01 0.01]);
  160. %axis image, colormap gray
  161. %title('3rd derivative map')
  162. %set(gca,'xticklabel',[],'yticklabel',[])
  163. %set(gca,'xtick',[],'ytick',[])
  164. %imrgb(:,:,1)=im+uint8(edge2nd*100);
  165. %imrgb(:,:,2)=im+uint8(edge3rd*100);
  166. %imrgb(:,:,3)=im;
  167. %subplot(1,3,3)
  168. %imshow(imrgb)
  169. %title('2nd and 3rd order zero crossing overlay')
  170. figure
  171. surf(mrsmap2,'edgecolor','none')
  172. colormap gray

Downloads

Depth-blur-Website.rar DEPTH-FROM-BLUR_ZAMAN_TIM_WEB.pdf]]>

You may also like...

2 Responses

  1. Rahma says:

    Hi, I’m interested to try your code above, but there is some error when I run your Elder Zucker algorithm. Here is the error:
    ??? Error using ==> iptcheckinput
    Function EDGE expected its first input, I,
    to be one of these types:
    double, single, uint8, uint16, uint32, int8, int16, int32
    Instead its type was logical.
    Error in ==> edge>parse_inputs at 564
    iptcheckinput(I,{‘numeric’},{‘nonsparse’,’2d’},mfilename,’I’,1);
    Error in ==> edge at 197
    [a,method,thresh,sigma,thinning,H,kx,ky] = parse_inputs(varargin{:});
    Error in ==> elderzucker_tim at 86
    edge2nd=edge2nd+edge(r2theta>0,’canny’)*sigma1(i);
    I wonder if you could help me with this trouble, thank you very much 🙂

  2. nor says:

    Do you have the any source code for the video?