Code covered by the BSD License  

Highlights from
deprecated -- Light Field Toolbox v0.2 -- v0.3 now available

image thumbnail

deprecated -- Light Field Toolbox v0.2 -- v0.3 now available

by

 

26 Apr 2013 (Updated )

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

LFCalRectifyLF( LF, CalInfo, RectOptions )
% LFCalRectifyLF - rectify a light field using a calibrated camera model
%
% Usage: 
%     [LF, RectOptions] = LFCalRectifyLF( LF, CalInfo, RectOptions )
%     [LF, RectOptions] = LFCalRectifyLF( LF, CalInfo )
% 
%
% This function is called by LFUtilDecodeLytroFolder to rectify a light field. It follows the
% rectification procedure described in:
%
% D. G. Dansereau, O. Pizarro, and S. B. Williams, "Decoding, calibration and rectification for
% lenselet-based plenoptic cameras," in Computer Vision and Pattern Recognition (CVPR), IEEE
% Conference on. IEEE, Jun 2013.
%
% Minor differences from the paper: light field indices [i,j,k,l] are 1-based in this
% implementation, and not 0-based as described in the paper.
% 
% Note that a calibration should only be applied to a light field decoded using the same lenselet
% grid model. This is because the lenselet grid model forms an implicit part of the calibration,
% and changing it will invalidate the calibration.
% 
% Inputs:
% 
%         LF : The light field to rectify; should be floating point
%
% 
%     CalInfo struct contains a calibrated camera model, see LFUtilCalLenseletCam:
%
%            .EstCamIntrinsicsH : 5x5 homogeneous matrix describing the lenselet camera intrinsics
% 
%            .EstCamDistortionV : Estimated distortion parameters
%  
% 
%    [optional] RectOptions struct (all fields are optional) :
% 
%        .NInverse_Distortion_Iters : Number of iterations in inverse distortion estimation
% 
%                        .Precision : 'single' or 'double'
% 
%               .RectCamIntrinsicsH : Requests a specific set of intrinsics for the rectified light
%                                     field. By default the rectified intrinsic matrix is
%                                     automatically constructed from the calibrated intrinsic
%                                     matrix, but this process can in some instances yield poor
%                                     results: excessive black space at the edges of the light field
%                                     sample space, or excessive loss of scene content off the edges
%                                     of the space. This parameters allows you to fine-tune the
%                                     desired rectified intrinsic matrix.
%
% Outputs :
% 
%     LF : The rectified light field
% 
%     RectOptions : The rectification options as applied, including any default values employed.
% 
% 
% See also:  LFUtilCalLenseletCam, LFUtilDecodeLytroFolder, LFUtilProcessWhiteImages

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

function [LF, RectOptions] = LFCalRectifyLF( LF, CalInfo, RectOptions )

%---Defaults---
RectOptions = LFDefaultField( 'RectOptions', 'Precision', 'single' );
RectOptions = LFDefaultField( 'RectOptions', 'NInverse_Distortion_Iters', 2 );

% The default rectified intrinsic matrix
if( ~isfield(RectOptions, 'RectCamIntrinsicsH') )
    fprintf('Generating rectified intrinsic matrix...\n');
    RectOptions.RectCamIntrinsicsH = CalInfo.EstCamIntrinsicsH;
    
    % Assure square pixels in ST and UV by taking the mean slopes in those pairs. Ideally the values
    % will use the hypercube-shaped light field space efficiently, leaving few black pixels at the
    % edges, and bumping as little information /outside/ the hypercube as possible. All this while
    % maintaining square pixels is s,t and in u,v, of course. Direct control is possible by
    % providing RectOptions.RectCamIntrinsicsH.
    ST_ST_Slope = mean([CalInfo.EstCamIntrinsicsH(1,1), CalInfo.EstCamIntrinsicsH(2,2)]);
    ST_UV_Slope = mean([CalInfo.EstCamIntrinsicsH(1,3), CalInfo.EstCamIntrinsicsH(2,4)]);
    UV_ST_Slope = mean([CalInfo.EstCamIntrinsicsH(3,1), CalInfo.EstCamIntrinsicsH(4,2)]);
    UV_UV_Slope = mean([CalInfo.EstCamIntrinsicsH(3,3), CalInfo.EstCamIntrinsicsH(4,4)]);
    
    RectOptions.RectCamIntrinsicsH(1,1) = ST_ST_Slope;
    RectOptions.RectCamIntrinsicsH(2,2) = ST_ST_Slope;
    RectOptions.RectCamIntrinsicsH(1,3) = ST_UV_Slope;
    RectOptions.RectCamIntrinsicsH(2,4) = ST_UV_Slope;
    
    RectOptions.RectCamIntrinsicsH(3,1) = UV_ST_Slope;
    RectOptions.RectCamIntrinsicsH(4,2) = UV_ST_Slope;
    RectOptions.RectCamIntrinsicsH(3,3) = UV_UV_Slope;
    RectOptions.RectCamIntrinsicsH(4,4) = UV_UV_Slope;
    
    %---force s,t translation to CENTER---
    RectOptions.RectCamIntrinsicsH = LFRecenterIntrinsics( RectOptions.RectCamIntrinsicsH, size(LF) );
end


%---Build interpolation indices---
fprintf('Generating interpolation indices...\n');
NChans = size(LF,5);
LF = cast(LF, RectOptions.Precision);

t_in=cast(1:size(LF,1), 'uint16'); % saving some mem by using uint16
s_in=cast(1:size(LF,2), 'uint16');
v_in=cast(1:size(LF,3), 'uint16');
u_in=cast(1:size(LF,4), 'uint16');

[tt,ss,vv,uu] = ndgrid(t_in,s_in,v_in,u_in);
% InterpIdx initially holds the index of the desired ray, and is evolved through the application
% of the inverse distortion model to eventually hold the continuous-domain index of the undistorted
% ray, and passed to the interpolation step.
InterpIdx = [ss(:)'; tt(:)'; uu(:)'; vv(:)'; ones(size(ss(:)'))];
DestSize = size(tt);
clear tt ss vv uu

%---Cast the to the required precision---
InterpIdx = cast(InterpIdx, RectOptions.Precision);

%---Convert the index of the desired ray to a ray representation using ideal intrinsics---
InterpIdx = RectOptions.RectCamIntrinsicsH * InterpIdx;

%---Apply inverse lens distortion to yield the undistorted ray---
k1 = CalInfo.EstCamDistortionV(1); % r^2
k2 = CalInfo.EstCamDistortionV(2); % r^4
k3 = CalInfo.EstCamDistortionV(3); % r^6
b1 = CalInfo.EstCamDistortionV(4); % decentering of lens distortion
b2 = CalInfo.EstCamDistortionV(5); % decentering of lens distortion

InterpIdx(3:4,:) = bsxfun(@minus, InterpIdx(3:4,:), [b1; b2]); % decentering of lens distortion

%---Iteratively estimate the undistorted direction----
DesiredDirection = InterpIdx(3:4,:);
for( InverseIters = 1:RectOptions.NInverse_Distortion_Iters )
    R2 = sum(InterpIdx(3:4,:).^2); % compute radius^2 for the current estimate
    % update estimate based on inverse of distortion model
    InterpIdx(3:4,:) = DesiredDirection ./ repmat((1 + k1.*R2 + k2.*R2.^2 + k3.*R2.^4),2,1);
end
clear R2 DesiredDirection

InterpIdx(3:4,:) = bsxfun(@plus, InterpIdx(3:4,:), [b1; b2]); % decentering of lens distortion

%---Convert the undistorted ray to the corresponding index using the calibrated intrinsics---
% todo[optimization]: The variable InterpIdx could be precomputed and saved with the calibration
InterpIdx = CalInfo.EstCamIntrinsicsH^-1 * InterpIdx;


%---Interpolate the required values---
fprintf('Interpolating...');
for( ColChan = 1:NChans )
    % todo[optimization]: use a weighted interpolation scheme to exploit the weight channel
    InterpSlice = interpn(squeeze(LF(:,:,:,:,ColChan)), InterpIdx(2,:),InterpIdx(1,:), InterpIdx(4,:),InterpIdx(3,:), 'linear');
    LF(:,:,:,:,ColChan) = reshape(InterpSlice, DestSize);
    fprintf('.');    
end

%---Clip interpolation result, which sometimes rings slightly out of range---
LF(isnan(LF)) = 0;
LF = max(0, min(1, LF));

fprintf('\nDone\n');

Contact us