% written by Christopher Hahne % for the University of Bedfordshire % Last update 13/01/16 % Matlab version R2015b (8.6.0.267246) % % File description: % This function takes optical parameters of a standard plenoptic camera % and computes the distance at which a synthetically focused image % exhibits least blur. The measurement unit is millimetre. % % For further assistance please contact info[at]christopherhahne.de % % command execution example: % refocDist(.009,2.75,0,.125,111.0324,193.2935,-65.5563,193.2935,1,0,'on') % % pp := pixel pitch % fs := focal length micro lens % hh := space between micro lens principal planes % pm := microlens pitch / diameter % dA := distance from MLA to main lens exit pupil % fU := focal length main lens % HH := space between main lens principal planes % bU := separation between micro lens plane and H plane of main lens % a := refocusing parameter for respective object plane function [output] = refocDist(pp,fs,hh,pm,dA,fU,HH,bU,a,crop,draw) D = 72; % entrance pupil diameter of main lens % (s,u) coordinates for the intersecting rays umax = floor(pm/pp); smax = 2*umax+1; sc = (smax-1)/2; if mod(umax,2) ~= 1 warning('Given parameters yield even micro image resolution.'); umax = umax-1; end % analyse crop value which minimizes micro image resolution if (~exist('crop','var')) crop = 0; elseif crop > (umax-1)/2-1 || crop < 0 crop = (umax-1)/2-1; warning('Crop value is set to %s pixels.',num2str(crop)); else crop = round(crop); end % reduce micro image size if cropped umax = umax-2*crop; % warning if fU > bU warning('Image distance is smaller than focal length (fU>bU).') elseif a >= (smax-1)/2 warning('Refocusing slice is out of range.') end % align intersection to be as paraxial as possible c = (umax-1)/2; sjShift = -round(a*(umax-1-2*crop)/2); % set starting position for micro lens s and pixel u j0 = sjShift; s0 = j0*pm; mc0 = -s0/dA; uc0 = -mc0*fs+s0; i0 = c; u0 = uc0+i0*pp; mij0 = (s0-u0)/fs; Uij0 = mij0*bU+s0; Fij0 = mij0*fU; qij0 = (Fij0-Uij0)/fU; % neighbour for shift+a slice j2 = a*(umax-1)+sjShift; s2 = j2*pm; mc2 = -s2/dA; uc2 = -mc2*fs+s2; i2 = -c; u2 = uc2+i2*pp; mij2 = (s2-u2)/fs; Uij2 = mij2*bU+s2; Fij2 = mij2*fU; qij2 = (Fij2-Uij2)/fU; % intersection of ray lines from pixel centres into object space funD = @(x) (qij0*x+Uij0) - (qij2*x+Uij2); % intersection of ray lines behind image sensor funDnew = @(x) (mij0*x+s0) - (mij2*x+s2); bNew = bU - fzero(funDnew,0); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Depth of Field % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u0U = u0+pp/2; u0L = u0-pp/2; m0U = (s0+pm/2-u0U)/fs; m0L = (s0-pm/2-u0L)/fs; mij0U = (s0-u0U)/fs; mij0L = (s0-u0L)/fs; Uij0U = mij0U*bU+s0+pm/2; Uij0L = mij0L*bU+s0-pm/2; Fij0U = mij0U*fU; Fij0L = mij0L*fU; qij0U = (Fij0U - Uij0U)/fU; qij0L = (Fij0L - Uij0L)/fU; u2U = u2+pp/2; u2L = u2-pp/2; m2U = (s2+pm/2-u2U)/fs; m2L = (s2-pm/2-u2L)/fs; mij2U = (s2-u2U)/fs; mij2L = (s2-u2L)/fs; Uij2U = mij2U*bU+s2+pm/2; Uij2L = mij2L*bU+s2-pm/2; Fij2U = mij2U*fU; Fij2L = mij2L*fU; qij2U = (Fij2U - Uij2U)/fU; qij2L = (Fij2L - Uij2L)/fU; funDm = @(x) (qij2U*x+Uij2U) - (qij0L*x+Uij0L); funDp = @(x) (qij0U*x+Uij0U) - (qij2L*x+Uij2L); % intersection of ray lines behind image sensor funDmNew = @(x) (mij0L*x+s0-pm/2) - (mij2U*x+s2+pm/2); bMinusNew = bU - fzero(funDmNew,0); funDpNew = @(x) (mij0U*x+s0+pm/2) - (mij2L*x+s2-pm/2); bPlusNew = bU - fzero(funDpNew,0); % check if plane is at infinity if (fU == bU || fU > bU) && a <= 0 d = inf; dY = qij2*(d-fs-hh-bU-HH) + Uij2; dPlus = inf; dMinus = fzero(funDm,fs+hh+bU+HH+fU) + bU + HH; if dMinus < 0 dMinus = inf; end dNew = inf; dPlusNew = inf; dMinusNew = (1/fU-1/bMinusNew)^-1 + bU + HH; else if fU > bNew % check if plane goes beyond infinity warning('Refocused object plane goes beyond infinity.'); end d = fzero(funD,fs+hh+bU+HH+fU) + bU + HH; if d < 0 d = inf; end dY = qij2*(d-fs-hh-bU-HH) + Uij2; dPlus = fzero(funDp,fs+hh+bU+HH+fU) + bU + HH; if dPlus < 0 dPlus = inf; end dMinus = fzero(funDm,fs+hh+bU+HH+fU) + bU + HH; if dMinus < 0 dMinus = inf; end dNew = (1/fU-1/bNew)^-1 + bU + HH; dPlusNew = (1/fU-1/bPlusNew)^-1 + bU + HH; dMinusNew = (1/fU-1/bMinusNew)^-1 + bU + HH; end % Depth of field DoF = dPlus - dMinus; % comparison of image and object side approach if round(d)~=round(dNew) || round(dPlus)~=round(dPlusNew) ... || round(dMinus)~=round(dMinusNew) warning('Image and object side approach yield different results.'); end % output vector output = [d dMinus dPlus]; %%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % plot diagram % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%% if (~exist('draw','var')) draw = 'off'; end if strcmp(draw,'on') % plot rays on an optical bench plotLim = 5000; % set 50 meter as maximum plot distance if dPlus > 0 && dPlus ~= Inf largeVal = dPlus; elseif dMinus > 0 && dMinus < plotLim largeVal = dMinus; else largeVal = plotLim; end xmax = round( largeVal + largeVal/10 ); sampleDensity = xmax/5000; x = 0:sampleDensity:xmax; % rays from sensor to micro lens h plane y0(x<=fs) = mij0*x(x<=fs) + u0; y2(x<=fs) = mij2*x(x<=fs) + u2; y0Aux(x<=fs) = mc0*x(x<=fs) + uc0; y2Aux(x<=fs) = mc2*x(x<=fs) + uc2; y0(x>fs&x<=fs+hh) = s0; y2(x>fs&x<=fs+hh) = s2; y0Aux(x>fs&x<=fs+hh) = s0; y2Aux(x>fs&x<=fs+hh) = s2; % rays from micro lens h to H1 plane of main lens (thin lens) y0(x>fs+hh&x<=fs+hh+bU) = ... mij0*x(x>fs+hh&x<=fs+hh+bU) + s0 - mij0*(hh+fs); y2(x>fs+hh&x<=fs+hh+bU) = ... mij2*x(x>fs+hh&x<=fs+hh+bU) + s2 - mij2*(hh+fs); y0Aux(x>fs+hh&x<=fs+hh+bU) = ... mc0*x(x>fs+hh&x<=fs+hh+bU) + s0 - mc0*(hh+fs); y2Aux(x>fs+hh&x<=fs+hh+bU) = ... mc2*x(x>fs+hh&x<=fs+hh+bU) + s2 - mc2*(hh+fs); % rays in space between principal planes of main lens y0(x>fs+hh+bU&x<=fs+hh+bU+HH) = Uij0; y2(x>fs+hh+bU&x<=fs+hh+bU+HH) = Uij2; y0Aux(x>fs+hh+bU&x<=fs+hh+bU+HH) = ... mc0*x(x>fs+hh+bU&x<=fs+hh+bU+HH) + mc0*HH; y2Aux(x>fs+hh+bU&x<=fs+hh+bU+HH) = ... mc2*x(x>fs+hh+bU&x<=fs+hh+bU+HH) + mc2*HH; % rays from main lens plane to object space y0(x>fs+hh+bU+HH) = ... qij0*x(x>fs+hh+bU+HH) + Uij0 - (qij0*(fs+hh+bU+HH)); y2(x>fs+hh+bU+HH) = ... qij2*x(x>fs+hh+bU+HH) + Uij2 - (qij2*(fs+hh+bU+HH)); y0Aux(x>fs+hh+bU) = 0*x(x>fs+hh+bU) + 0 - (0*(fs+hh+bU)); y2Aux(x>fs+hh+bU) = 0*x(x>fs+hh+bU) + 0 - (0*(fs+hh+bU)); % DoF rays yUp(x<=fs) = m0U*x(x<=fs) + u0U; yUp(x>fs&x<=fs+hh) = s0+pm/2; yUp(x>fs+hh&x<=fs+hh+bU) = ... mij0U*x(x>fs+hh&x<=fs+hh+bU) + s0+pm - m0U*fs - mij0U*hh; yUp(x>fs+hh+bU&x<=fs+hh+bU+HH) = Uij0U; yUp(x>fs+hh+bU+HH) = ... qij0U*x(x>fs+hh+bU+HH) + Uij0U - (qij0U*(fs+hh+bU+HH)); yLm(x<=fs) = m0L*x(x<=fs) + u0L; yLm(x>fs&x<=fs+hh) = s0-pm/2; yLm(x>fs+hh&x<=fs+hh+bU) = ... mij0L*x(x>fs+hh&x<=fs+hh+bU) + s0-pm - m0L*fs - mij0L*hh; yLm(x>fs+hh+bU&x<=fs+hh+bU+HH) = Uij0L; yLm(x>fs+hh+bU+HH) = ... qij0L*x(x>fs+hh+bU+HH) + Uij0L - (qij0L*(fs+hh+bU+HH)); yUm(x<=fs) = m2U*x(x<=fs) + u2U; yUm(x>fs&x<=fs+hh) = s2+pm/2; yUm(x>fs+hh&x<=fs+hh+bU) = ... mij2U*x(x>fs+hh&x<=fs+hh+bU) + s2+pm - m2U*fs - mij2U*hh; yUm(x>fs+hh+bU&x<=fs+hh+bU+HH) = Uij2U; yUm(x>fs+hh+bU+HH) = ... qij2U*x(x>fs+hh+bU+HH) + Uij2U - (qij2U*(fs+hh+bU+HH)); yLp(x<=fs) = m2L*x(x<=fs) + u2L; yLp(x>fs&x<=fs+hh) = s2-pm/2; yLp(x>fs+hh&x<=fs+hh+bU) = ... mij2L*x(x>fs+hh&x<=fs+hh+bU) + s2-pm - m2L*fs - mij2L*hh; yLp(x>fs+hh+bU&x<=fs+hh+bU+HH) = Uij2L; yLp(x>fs+hh+bU+HH) = ... qij2L*x(x>fs+hh+bU+HH) + Uij2L - (qij2L*(fs+hh+bU+HH)); da=strcat('$d_{',num2str(a),'}$'); daMinus=strcat('$d_{',num2str(a),'-}$'); daPlus=strcat('$d_{',num2str(a),'+}$'); yDim = D/12; padding = D/36; % labels da, da-, da+ plot(d, dY,'ko',[0, round((d+20)/100)*100],[0, 0], ... 0,uc0,'ko',0,uc2,'ko',0,u0,'go',0,u2,'bo'); text(d+DoF/16,yDim-padding+1,da,'FontSize',16,'color','c', ... 'interpreter','latex') text(dMinus-dMinus/100,yDim-padding+1,daMinus,'FontSize',16, ... 'color','r','interpreter','latex', ... 'HorizontalAlignment','right') text(dPlus+DoF/16,yDim-padding+1,daPlus,'FontSize',16, ... 'color','k','interpreter','latex') % diagram title, labels, size and attributes h = gca; grid off, axis on xlabel('$z_U$~[mm]','FontSize',20,'interpreter','latex') ylabel('$(u,s)$~[mm]','FontSize',20,'interpreter','latex') axis([0 xmax -yDim-padding+2 yDim+padding-2]); set(h,'FontSize',18) % DoF plots line('XData', x, 'YData', yLm, 'LineStyle', '-', 'Color', 'r'); line('XData', x, 'YData', yUm, 'LineStyle', '-', 'Color', 'r'); line('XData', x, 'YData', yUp, 'LineStyle', '-', 'Color', 'k'); line('XData', x, 'YData', yLp, 'LineStyle', '-', 'Color', 'k'); % intersecting rays for center of corresponding refocusing plane line('XData', x, 'YData', y0Aux, 'LineStyle', '-', 'Color', 'y'); line('XData', x, 'YData', y2Aux, 'LineStyle', '-', 'Color', 'y'); line('XData', x, 'YData', y0, 'LineStyle', '-', 'Color', 'g'); line('XData', x, 'YData', y2, 'LineStyle', '-', 'Color', 'b'); % intersection planes line('XData', [d d], 'YData', ... [yDim+padding -yDim-padding], 'LineStyle', '--', 'Color', 'c') line('XData', [dPlus dPlus], 'YData', ... [yDim+padding -yDim-padding], 'LineStyle', '--', 'Color', 'k') line('XData', [dMinus dMinus], 'YData', ... [yDim+padding -yDim-padding], 'LineStyle', '--', 'Color', 'r') % sensor line('XData', [0 0], 'YData', [sc*pm -sc*pm], ... 'LineStyle', '-', 'LineWidth', 0.25); pixelY0 = uc0-pp/2-ceil(pm/pp)/2*pp:pp:uc0+pp/2+ceil(pm/pp)/2*pp; pixelY2 = uc2-pp/2-ceil(pm/pp)/2*pp:pp:uc2+pp/2+ceil(pm/pp)/2*pp; pixelX = zeros(size(pixelY0)); line('XData', pixelX, 'YData', pixelY0, ... 'Marker', '+', 'MarkerSize', 3); line('XData', pixelX, 'YData', pixelY2, ... 'Marker', '+', 'MarkerSize', 3); % micro lens grid lensY = -sc*pm+pm/2:pm:sc*pm+pm/2; lensf = -sc*pm:pm:sc*pm; lensX = (fs+hh) * ones(size(lensY)); line('XData', [fs fs], 'YData', [sc*pm -sc*pm], ... 'LineStyle', '--', 'LineWidth', 0.25); line('XData', [fs+hh fs+hh], 'YData', [sc*pm -sc*pm], ... 'LineStyle', '-', 'LineWidth', 0.25); line('XData', lensX, 'YData', lensY, 'Marker', '+', ... 'MarkerSize', 4); % lens borders line('XData', lensX, 'YData', lensf, 'Marker', '.', ... 'MarkerSize', 6); % optical axis % main lens principal planes line('XData', [fs+hh+bU+HH, fs+hh+bU+HH], 'YData', ... [-D/2, D/2], 'LineStyle', '--', 'Color', 'k'); line('XData', [fs+hh+bU, fs+hh+bU], 'YData', ... [-D/2, D/2], 'LineStyle', '--', 'Color', 'k'); text(fs+hh+bU+HH+2,yDim-padding+1,'$$H_{1U}$$', ... 'FontSize',12,'color','k','interpreter','latex') text(fs+hh+bU+2,yDim-padding+1,'$$H_{2U}$$', ... 'FontSize',12,'color','k','interpreter','latex') % main lens focal point line('XData', [fs+hh+bU+HH+fU, fs+hh+bU+HH+fU], 'YData', ... [-D/20, D/20], 'LineStyle', ':', 'Color', 'k'); % optical axis zU line('XData', [0, xmax], 'YData', [0, 0], ... 'LineStyle', '-', 'Color', 'k'); end end