[Vision] Depth estimation from Blur estimation

header_blurdepth

In this post I will briefly demonstrate how we can compute the depth estimation from ‘focal’ blur estimation. I will give a brief introduction, show some formula’s, and then give the Matlab sample code and ofcourse some nice downloads to try yourself (in Matlab).

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).

Test image (35mm f1.4)

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.

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

Downloads

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

Tim Zaman

Engineer from The Netherlands (1988). Living in Delft. MSc degree in Robotics. Is on a PhD research position on development of imaging devices. For contact see 'About' in the above menu.

You may also like...

1 Response

  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 :)

Leave a Reply

Your email address will not be published. Required fields are marked *


× six = 18

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre lang="" line="" escaped="" cssfile="">