[Vision] Depth estimation from Blur estimation

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

%Tim Zaman, TU Delft ME, 01-01-2012
%This code is licensed under a 
%Creative Commons Attribution 3.0 Unported License.
 
clc;clear all;close all
 
im=imread('imnikon.jpg');
im=im(:,:,1);
 
sigmaA=8;
sigmaB=10;
sigmaMax=max([sigmaA sigmaB]);
fsz=[sigmaMax, sigmaMax];
 
kernelA=fspecial('gaussian', fsz, sigmaA);
kernelB=fspecial('gaussian', fsz, sigmaB);
%imA = filter2(kernelA, im);
%imB = filter2(kernelB, im);
imA = imfilter(im, kernelA, 'symmetric', 'conv');
imB = imfilter(im, kernelB, 'symmetric', 'conv');
 
 
R1=single(im)-single(imA);
%R1=((single(imA)./single(im))-1).*single(imB);
R2=single(imA)-single(imB);
 
%R1=abs(R1);
%R2=abs(R2);
%R2(R2<1)=1;
 
R=R1./R2;
 
R(isnan(R))=0;
R(isinf(R))=0;
Rf=maxfilt2(R,8);
%Rf=max(max(Rf))-Rf;
 
blurmap=(sigmaA*sigmaB)./((sigmaB-sigmaA)*Rf+sigmaB);
 
figure;
 
subplot(3,3,1), imshow(im), title('original')
subplot(3,3,2), imshow(uint8(imA)), title('imA')
subplot(3,3,3), imshow(uint8(imB)), title('imB')
subplot(3,3,4), imagesc(R1), title('R1')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
subplot(3,3,5), imagesc(R2), title('R2')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
subplot(3,3,6), imagesc(R), title('R')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
subplot(3,3,7), imagesc(blurmap), title('blurmap')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
subplot(3,3,8), imagesc(blurmap>1 & blurmap<4), title('\sigma [2:4]')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
subplot(3,3,9), imagesc(blurmap>0 & blurmap<0.4), title('\sigma [0:0.5]')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image

Elder and Zucker algorithm

%Tim Zaman, TU Delft ME, 01-01-2012
%This code is licensed under a 
%Creative Commons Attribution 3.0 Unported License.
 
clc;close all;clear all
 
%Read image
im = imread('imnikon.jpg');
im = im(:,:,1);
 
%Some variables
sn = 3; %image noise [sigma]
alphap = 2e-7; %2e-7
 
%Apply white noise sn to image (optional)
%im = imnoise(im,'gaussian',0,(sn/255)^2);
imd=double(im);
 
sigma1 = [8 4 2 1 0.5]; %Scale space octaves for 1st deriv
sigma2 = sigma1;
 
%Preallocate some matrices
mrsmap1=ones(size(im))*4;
mrsmap2=ones(size(im))*4;
edge2nd=zeros(size(im));
edge3rd=zeros(size(im));
 
%Make steerable gaussian 1st order derivative filters
for i=1:length(sigma1)
    %First derivative part
 
    %Create normal gaussian filter kernel
    h = fspecial('gaussian', round(10*sigma1(i)), sigma1(i));
    %Create gaussian 1st order gradient in X and Y direction
    [g1x,g1y] = gradient(h);
    %Apply the filters on the image
    r1x=imfilter(imd,g1x,'replicate');
    r1y=imfilter(imd,g1y,'replicate');
    thetam=atan2(r1y, r1x);
    r1theta=cos(thetam).*r1x + sin(thetam).*r1y;
    r1theta_or=r1theta;
    %Compute critical value
    g1theta2=1/(2*sqrt(2*pi)*sigma1(i)^2);
    s1=g1theta2*sn;
    c1=s1*sqrt(-2*log(alphap));
    %Make the binary 1st deriv reliable map
    reliable1=(abs(r1theta) >= c1);
 
    %Put into the Minimum Reliable Scale map (1st)
    mrsmap1(reliable1)=sigma1(i);
    %Apply threshold on r1theta and thetam
    thetam=thetam.*reliable1;
    r1theta=r1theta.*reliable1;
 
 
    %Second derivative part
 
    %Make 2nd order gaussian filter kernels
%    [g2x,g2xy] = gradient(gradient(h));
%    g2y=rot90(g2x);
    %Apply 2nd order gaussian filter kernels to image
%    r2x=imfilter(imd,g2x,'replicate');
%    r2y=imfilter(imd,g2y,'replicate');
%    r2xy=imfilter(imd,g2xy,'replicate');
    [r2x,r2y] = gradient(r1theta_or);
 
    %Compute 2nd derivative in maginitude direction
    r2theta=cos(thetam).*r2x + sin(thetam).*r2y;
    r2theta_or=r2theta;
    %r2theta =(cos(thetam).^2 .* r2x)+(sin(thetam).^2 .* r2y)...
    %    -(2.*cos(thetam).*sin(thetam).*r2xy);
    g22=1/(4*sqrt(pi/3)*sigma2(i)^3);
    s2=g22*sn;
    c2=sqrt(2)*s2*erfinv(1-alphap); %paper says (alpha)???
    %Make the binary 2nd deriv reliable map
    reliable2=(abs(r2theta) >= c2);
    %Put into the Minimum Reliable Scale map (2nd)
    mrsmap2(reliable2)=sigma2(i);
    %Apply threshold on r2theta
    r2theta=r2theta.*reliable2;
    %Now remove the new lower sigma part for filling in
    edge2nd=edge2nd.*(1-reliable1.*reliable2);
 
    %Locate the edges at zero crossings in 
    % binary image (discrete for simplicity)
    edge2nd=edge2nd+edge(r2theta>0,'canny')*sigma1(i);
 
    %LOOK FOR ZERO CROSSINGS in [3x3] window
%     for x=1:(size(im,2)-2);
%         for y=1:(size(im,1)-2);
%             %[3x3] window
%             imw=r2theta(y:y+2,x:x+2);
%             %Check for sign changes within window
%             sp=sum(sum(imw>=0));
%             sm=sum(sum(imw<=0));
%             %Apply to edgemap
%             if (sp~=9 && sm~=9)
%                 edge2nd(y,x)=sigma1(i);
%             end
%         end
%     end
 
 
 
 
 
 
    %Now remove the new lower sigma part for filling in
    edge3rd=edge3rd.*(1-reliable1.*reliable2);
 
    %Make the third derivative
    [r3x,r3y]=gradient(r2theta_or);
    r3theta=cos(thetam).*r3x + sin(thetam).*r3y;
    edge3rd=edge3rd+edge(r3theta>0,'canny')*sigma1(i);
    %Now ignore the edges at faulty locations
    %edge3rd=edge3rd.*imerode(reliable1,ones(10));
 
 
    %Now, for every point in edge2nd, run along the direction
    % of the gradient (thetam) and find the first occurancy of
    % a zero in the third derivative map.
 
 
 
 
 
    %for ii=1:sum(sum(edge2nd>0)) %For all edgepoints..
%     for y=2:255;
%         x=r2theta(y,:);
%         yy=1:1/5:length(x);
%         %Interpolate the signal (simple solution)
%         p = spline([1:length(x)],x,yy);
%         [~,loc]=findpeaks(abs(p),'MINPEAKDISTANCE',50,'MINPEAKHEIGHT',max(p)*0.95);
%         d=max(abs(diff(loc)/5));
%         sigmab(y)=sqrt(((d/2)^2)-sigma2(i)^2);
%         %plot(x)
%         %pause(.01)
%     end
 
 
 
    %imagesc(r2theta)
    %colormap gray
 
%     figure
%     plot(sigmab,'-b')
%     hold on
%     plot([1:26.6/length(x):26.6],'--r')
%     hold off, grid on
%     legend('estimated','actual')
%     ylabel('Blur scale (px)')
%     xlabel('Arclength(px)')
%     title(['\sigma=' num2str(sn)])
    %pause(2)
    %figure
    %imagesc(r1theta)
    %colormap gray
    %figure
    %plot(r2theta(:,550))
    %pause(1)
end
 
figure
subplot(2,3,1), imshow(im)
title('original')
subplot(2,3,2), imagesc(sqrt(mrsmap1))
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
%colormap hot
title('mrs map 1st deriv')
subplot(2,3,3), imagesc(sqrt(mrsmap2))
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
%colormap hot
title('mrs map 2nd deriv')
 
subplot(2,3,6), imagesc(mrsmap2==0.5)
title('\sigma=0.5')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
subplot(2,3,5), imagesc(mrsmap2==1)
title('\sigma=1')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
subplot(2,3,4), imagesc(mrsmap2==8)
title('\sigma=8')
set(gca,'xticklabel',[],'yticklabel',[])
set(gca,'xtick',[],'ytick',[]),axis image
 
 
%subplot(1,3,1), [f1]=imagesc(r2theta,[-1 1]);
%axis image, colormap gray
%title('2nd derivative map')
%set(gca,'xticklabel',[],'yticklabel',[])
%set(gca,'xtick',[],'ytick',[])
%subplot(1,3,2), [f2]=imagesc(r3theta,[-0.01 0.01]);
%axis image, colormap gray
%title('3rd derivative map')
%set(gca,'xticklabel',[],'yticklabel',[])
%set(gca,'xtick',[],'ytick',[])
 
%imrgb(:,:,1)=im+uint8(edge2nd*100);
%imrgb(:,:,2)=im+uint8(edge3rd*100);
%imrgb(:,:,3)=im;
%subplot(1,3,3)
%imshow(imrgb)
%title('2nd and 3rd order zero crossing overlay')
 
figure
surf(mrsmap2,'edgecolor','none')
colormap gray