[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 […]

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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
%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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
%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

Downloads


About Tim Zaman

There's page about myself with some babble and pics on this site Here!. // R&D : Vision/Imaging | 2D/3D Image Processing | Robotics/Automation | Programming | Electronics | Open Source | And, obviously, Photography.

Creative Commons License
[Vision] Depth estimation from Blur estimation by Tim Zaman, unless otherwise expressly stated, is licensed under a Creative Commons Attribution 3.0 Unported License.