Code covered by the BSD License  

Highlights from
Light Field Toolbox v0.2

image thumbnail

Light Field Toolbox v0.2

by

 

26 Apr 2013 (Updated )

A set of tools for working with light field (aka plenoptic) imagery in Matlab

LFBuildLenseletGridModel( WhiteImg, GridModelOptions )
% LFBuildLenseletGridModel - builds a lenselet grid model from a white image
%
% Usage:
% 
%     [LenseletGridModel, GridCoords] = LFBuildLenseletGridModel( WhiteImg, GridModelOptions )
%
%
% Inputs: 
% 
%     WhiteImg : path of an image taken through a diffuser, or of an entirely white scene
%
% 
%     GridModelOptions : struct controlling the model-bulding options
% 
%          .ApproxLenseletSpacing : A rough initial estimate to initialize the lenselet spacing
%                                   computation
% 
%           .FilterDiskRadiusMult : Filter disk radius for prefiltering white image for locating
%                                   lenselets; expressed relative to lenselet spacing; e.g. a
%                                   value of 1/3 means a disk filte with a radius of 1/3 the
%                                   lenselet spacing
% 
%                        .CropAmt : Image edge pixels to ignore when finding the grid
% 
%                       .SkipStep : As a speed optimization, not all lenselet centers contribute
%                                   to the grid estimate; <SkipStep> pixels are skipped between
%                                   lenselet centers that get used; a value of 1 means use all
%
%           [optional] .Precision : 'single' or 'double'
% 
% 
% Outputs:
% 
%     LenseletGridModel : struct describing the lenselet grid
% 
%         .HSpacing, .VSpacing : Spacing between lenselets, in pixels
% 
%          .HOffset, .VOffset  : Offset of first lenselet, in pixels
% 
%                         .Rot : Rotation of lenselets, in radians
%
% 
%     GridCoords : a list of N x M x 2 pixel coordinates generated from the estimated
%                  LenseletGridModel, where N and M are the estimated number of lenselets in the
%                  horizontal and vertical directions, respectively.
%
% See also:  LFUtilProcessWhiteImages

% Part of LF Toolbox v0.2 released 27-May-2013
% Copyright (c) 2013, Donald G. Dansereau

function [LenseletGridModel, GridCoords] = LFBuildLenseletGridModel( WhiteImg, GridModelOptions )

%---Defaults---
GridModelOptions = LFDefaultField( 'GridModelOptions', 'Precision', 'single' );


% Try locating lenselets by convolving with a disk
% Also tested Gaussian... the disk seems to yield a stronger result
h = zeros( size(WhiteImg), GridModelOptions.Precision );
hr = fspecial( 'disk', GridModelOptions.ApproxLenseletSpacing * GridModelOptions.FilterDiskRadiusMult );
hr = hr ./ max( hr(:) ); 
hs = size( hr,1 );
DiskOffset = round( (size(WhiteImg) - hs)/2 );
h( (1:hs) + DiskOffset(1), (1:hs) + DiskOffset(2) ) = hr;

fprintf( 'Filtering...\n' );
% Convolve using fft
WhiteImg = fft2(WhiteImg);
h = fft2(h);
WhiteImg = WhiteImg.*h;
WhiteImg = ifftshift(ifft2(WhiteImg));
WhiteImg = (WhiteImg - min(WhiteImg(:))) ./ abs(max(WhiteImg(:))-min(WhiteImg(:)));
WhiteImg = cast(WhiteImg, GridModelOptions.Precision);
clear h

fprintf('Finding Peaks...\n');
% Find peaks in convolution... ideally these are the lenselet centers
Peaks = imregionalmax(WhiteImg);
PeakIdx = find(Peaks==1);
[PeakIdxY,PeakIdxX] = ind2sub(size(Peaks),PeakIdx);
clear Peaks

% Crop to central peaks; eliminates edge effects
InsidePts = find(PeakIdxY>GridModelOptions.CropAmt & PeakIdxY<size(WhiteImg,1)-GridModelOptions.CropAmt & ...
    PeakIdxX>GridModelOptions.CropAmt & PeakIdxX<size(WhiteImg,2)-GridModelOptions.CropAmt);
PeakIdxY = PeakIdxY(InsidePts);
PeakIdxX = PeakIdxX(InsidePts);

%---Form a Delaunay triangulation to facilitate row/column traversal---
Triangulation = DelaunayTri(PeakIdxY, PeakIdxX);

%---Traverse rows and columns of lenselets, collecting stats---

fprintf('Vertical fit...\n');
%--Traverse vertically--
YStart = GridModelOptions.CropAmt*2;
YStop = size(WhiteImg,1)-GridModelOptions.CropAmt*2;

XIdx = 1;
for( XStart = GridModelOptions.CropAmt*2:GridModelOptions.SkipStep:size(WhiteImg,2)-GridModelOptions.CropAmt*2 )
    CurPos = [XStart, YStart];
    YIdx = 1;
    while( 1 )
        ClosestLabel = nearestNeighbor(Triangulation, CurPos(2), CurPos(1));
        ClosestPt = [PeakIdxX(ClosestLabel), PeakIdxY(ClosestLabel)];

        RecPtsY(XIdx,YIdx,:) = ClosestPt;
        
        CurPos = ClosestPt;
        CurPos(2) = round(CurPos(2) + GridModelOptions.ApproxLenseletSpacing * sqrt(3));
        if( CurPos(2) > YStop )
            break;
        end
        YIdx = YIdx + 1;
    end
    %--Estimate angle for this most recent line--
    LineFitY(XIdx,:) = polyfit(RecPtsY(XIdx, 3:end-3,2), RecPtsY(XIdx, 3:end-3,1), 1);
    XIdx = XIdx + 1;
end

fprintf('Horizontal fit...\n');
%--Traverse horizontally--
XStart = GridModelOptions.CropAmt*2;
XStop = size(WhiteImg,2)-GridModelOptions.CropAmt*2;
YIdx = 1;
for( YStart = GridModelOptions.CropAmt*2:GridModelOptions.SkipStep:size(WhiteImg,1)-GridModelOptions.CropAmt*2 )
    CurPos = [XStart, YStart];
    XIdx = 1;
    while( 1 )
        ClosestLabel = nearestNeighbor(Triangulation, CurPos(2), CurPos(1));
        ClosestPt = [PeakIdxX(ClosestLabel), PeakIdxY(ClosestLabel)];

        RecPtsX(XIdx,YIdx,:) = ClosestPt;
        
        CurPos = ClosestPt;
        CurPos(1) = round(CurPos(1) + GridModelOptions.ApproxLenseletSpacing);
        if( CurPos(1) > XStop )
            break;
        end
        XIdx = XIdx + 1;
    end
    %--Estimate angle for this most recent line--
    LineFitX(YIdx,:) = polyfit(RecPtsX(3:end-3,YIdx,1), RecPtsX(3:end-3,YIdx,2), 1);
    YIdx = YIdx + 1;
end

%--Trim ends to wipe out alignment, initial estimate artefacts--
RecPtsY = RecPtsY(3:end-2, 3:end-2,:);  
RecPtsX = RecPtsX(3:end-2, 3:end-2,:);

%--Estimate angle--
SlopeX = mean(LineFitX(:,1));
SlopeY = mean(LineFitY(:,1));

AngleX = atan2(-SlopeX,1);
AngleY = atan2(SlopeY,1);
EstAngle = mean([AngleX,AngleY]);

%--Estimate spacing, assuming approx zero angle--
t=squeeze(RecPtsY(:,:,2));
YSpacing = diff(t,1,2);
YSpacing = mean(YSpacing(:))/2 / (sqrt(3)/2);

t=squeeze(RecPtsX(:,:,1));
XSpacing = diff(t,1,1);
XSpacing = mean(XSpacing(:));

%--Correct for angle--
XSpacing = XSpacing / cos(EstAngle);
YSpacing = YSpacing / cos(EstAngle);

%---Estimate offset---
EstOffset = zeros(1,2);

%--Build initial grid estimate, starting with CropAmt for the offsets--
LenseletGridModel = struct('HSpacing',XSpacing, 'VSpacing',YSpacing*sqrt(3)/2, 'HOffset',GridModelOptions.CropAmt, ...
    'VOffset',GridModelOptions.CropAmt, 'Rot',-EstAngle);
UMax = ceil( (size(WhiteImg,2)-GridModelOptions.CropAmt*2)/XSpacing );
VMax = ceil( (size(WhiteImg,1)-GridModelOptions.CropAmt*2)/YSpacing/(sqrt(3)/2) );
GridCoords = LFBuildHexGrid( LenseletGridModel, UMax, VMax );

%--Find offset to nearest peak for each--
GridCoordsX = GridCoords(:,:,1);
GridCoordsY = GridCoords(:,:,2);
BuildGridCoords = [GridCoordsY(:), GridCoordsX(:)];

IdealPts = nearestNeighbor(Triangulation, round(GridCoordsY(:)), round(GridCoordsX(:)));
IdealPtCoords = [PeakIdxY(IdealPts), PeakIdxX(IdealPts)];

%--Bring back to absolute pixel coords--
AllOffsets = BuildGridCoords - IdealPtCoords;
EstOffset = EstOffset + median(AllOffsets);
LenseletGridModel.HOffset = GridModelOptions.CropAmt - EstOffset(2);
LenseletGridModel.VOffset = GridModelOptions.CropAmt - EstOffset(1);

%--Remove crop offset--
LenseletGridModel.HOffset = mod( LenseletGridModel.HOffset, LenseletGridModel.HSpacing );
% be sure to shift an even number of vertical rows to ensure the hex grid phase is maintained
VSteps = floor( LenseletGridModel.VOffset / LenseletGridModel.VSpacing );
VSteps = 2*floor(VSteps / 2);
LenseletGridModel.VOffset = LenseletGridModel.VOffset - VSteps * LenseletGridModel.VSpacing;

%---visualize result---
UMax = ceil((size(WhiteImg,2))/LenseletGridModel.HSpacing);
VMax = ceil((size(WhiteImg,1))/LenseletGridModel.VSpacing);
GridCoords = LFBuildHexGrid( LenseletGridModel, UMax, VMax );

fprintf('...Done.\n');

end


function [GridCoords] = LFBuildHexGrid( LenseletGridModel, UMax, VMax )

RotCent = eye(3);
RotCent(1:2,3) = [LenseletGridModel.HOffset, LenseletGridModel.VOffset];

ToOffset = eye(3);
ToOffset(1:2,3) = [LenseletGridModel.HOffset, LenseletGridModel.VOffset];

R = ToOffset * RotCent * LFRotz(LenseletGridModel.Rot) * RotCent^-1;


[vv,uu] = ndgrid((0:VMax-1).*LenseletGridModel.VSpacing, (0:UMax-1).*LenseletGridModel.HSpacing);
uu(1:2:end,:) = uu(1:2:end,:) + 0.5.*LenseletGridModel.HSpacing;

GridCoords = [uu(:), vv(:), ones(numel(vv),1)];
GridCoords = (R*GridCoords')';

GridCoords = reshape(GridCoords(:,1:2), [VMax,UMax,2]);

end

Contact us