Code covered by the BSD License  

Highlights from
Affine optic flow

image thumbnail

Affine optic flow

by

 

26 Mar 2010 (Updated )

Estimates the parameters of an affine (first-order) optic flow model from two images.

affine_flow
classdef affine_flow
    % AFFINE_FLOW estimates affine (first-order) optic flow.
    %
    % Given two images, this class computes estimates of the parameters of
    % an affine (first-order) flow model which describes the overall
    % deformation between the two images.
    %
    % Affine flow
    % -----------
    %
    % An affine flow field has six parameters. One way to express these is
    % as follows. If the optic flow vector at an image location (x,y) is
    % (vx,vy), the first-order model is:
    %
    %           [vx vy] = [x y 1] [  d+s1  s2+r  ]
    %                             [  s2-r  d-s1  ]
    %                             [   vx0   vy0  ]
    %
    % Note that when x and y refer to pixels in the image array, the
    % conventional order is reversed, and the image intensity at (x,y) is
    % given by IMAGE(y,x) - that is, y is the row and x is the column.
    %
    % The parameters may be thought of thus: (vx0, vy0) is the optic flow
    % at the origin (one pixel outside the top left corner of the image);
    % d is the rate of dilation; r is the rate of rotation; s1 is the
    % shear along the main image axes; s2 is the shear along the diagonal
    % axes.
    %
    % Algorithm
    % ---------
    %
    % The flow is estimated by carrying out a least-squares fit of the
    % model to the spatial and temporal gradients of the image pair; this
    % is an extension of the Lucas-Kanade method. Gradients are estimated
    % with simple symmetrical differences and optional smoothing.
    %
    % Optionally, the images may be resampled on a log-polar grid, as
    % described here: http://www.bmva.org/bmvc/1994/bmvc-94-057.pdf. This
    % may increase the estimation accuracy because the grid spacing can be
    % matched to the likely magnitude of the optic flow. When log-polar
    % sampling is used, an iterative scheme to refine the estimates by
    % moving the sampling centre with the flow may be activated.
    %
    % The simplest use of the class, with no properties set explicitly, is
    % unlikely to yield good results on real images. At a minimum,
    % smoothing of the images is needed during computation of the
    % gradients. For more information, see the help for the individual
    % properties.
    %
    % Usage
    % -----
    %
    % This is a value class, so any update to the state of an affine_flow
    % object requires an assignment.
    %
    % A new affine_flow object is set up using parameter/value arguments,
    % for example:
    %
    %   af = affine_flow('sigmaXY', 10, 'sampleStep', 2);
    %
    % Parameters may also be updated subsequently. Parameters that are not
    % set are given defaults. To see the property names and default values
    % for rectilinear sampling, do
    %
    %   disp(affine_flow);
    %
    % To see the property names and default values for log-polar sampling,
    % do
    %
    %   disp(affine_flow('sampleMethod', 'logpolar'));
    %
    % The images to be processed may be supplied to the constructor, or
    % subsequently:
    %
    %   af.image1 = im1;
    %   af.image2 = im2;
    %
    % The flow is computed by a call to the findFlow method:
    %
    %   af = af.findFlow;
    %
    % and the results obtained from the flowStruct property:
    %
    %   flow = af.flowStruct;
    %
    % See the individual properties and methods for more details, and
    % information about the results structure.
    %
    % Efficiency notes
    % ----------------
    %
    % If many pairs of images are to be processed, it is most efficient to
    % create a single affine_flow object outside the loop and use it to
    % process each pair of images.
    %
    % If multiple images are to be matched to a given initial image, only
    % the image2 property should be updated within the loop, to avoid
    % refiltering the first image.
    %
    % If the flow between successive pairs of a sequence of images is
    % required, it is most efficient to call the advance method within the
    % loop to move image2 to image1, and only to update the image2 property
    % directly.
    %
    % See also: AFFINE_FLOWDEMO, LOGSAMPLE, GRADIENTS_XYT, 
    % AFFINE_FLOWDISPLAY, AFFINE_FLOWEDGEDISPLAY

    % Copyright David Young 2010

    
    properties (Dependent)
        
        % Sampling method
        %
        % Values are 'rectilinear' or 'logpolar'; the default is
        % 'rectilinear'. The image gradients are sampled either on a
        % regular rectilinear grid aligned with the image axes, or on a
        % log-polar grid. The possible advantage of a log-polar grid is
        % that the sample spacing increases with distance from the centre,
        % so provided the optic flow is small at the grid centre, or the
        % grid is moved to track the motion, the sample spacing will be
        % better matched to the flow velocity variation across the image.
        %
        % See also logsample
        sampleMethod
        
        % General region of interest
        %
        % A binary array (double or logical) the same size as the images,
        % used as a mask, or the empty matrix. Gradients are only computed
        % for pixels where there is a non-zero value in the roi array. The
        % smoothing process may access pixels outside the ROI mask. Ignored
        % if set to the empty matrix (the default).
        %
        % There are simpler ways to restrict processing to part of the
        % image: for rectilinear sampling, a rectangular region may be
        % specified with the rectRegion property; for log-polar sampling, a
        % circular region may be specified with logCentre and logRmax.
        %
        % Setting regionOfInterest will discard any previous value for the
        % rectRegion property. These two properties should not both be set
        % in a call to the class constructor.
        regionOfInterest
        
        % Rectangular region of interest (rectilinear sampling only)
        %
        % A 4-element row vector of integers, or the empty matrix. This
        % specifies a rectangular region of the image in which to compute
        % the flow, in the form [ROWMIN ROWMAX COLMIN COLMAX]. Because this
        % resembles array indexing, (row,column) rather than (x,y) order is
        % used. Gradients are computed only for pixels in the region, but
        % the smoothing process may access pixels just outside the region. 
        % 
        % If set to the empty matrix (the default), the region used is that
        % for which the smoothing and differencing convolutions, applied to
        % the whole image, produce 'valid' results (i.e. results which do
        % not require a guess at what lies outside the image boundary).
        %
        % Setting rectRegion will discard any previous value for the
        % regionOfInterest property. These two properties should not both
        % be set in a call to the class constructor.
        rectRegion
        
        % Smoothing constant for rectilinear sampling
        %
        % A non-negative real scalar; default 0. The parameter for the
        % Gaussian mask used to smooth the images for gradient estimation.
        %
        % See also GSMOOTH2
        sigmaXY
        
        % Sampling interval (rectilinear sampling only)
        %
        % A positive integer; default 1. The system of equations that is
        % solved is heavily overdetermined and it is not generally
        % necessary to use gradients for every pixel for the least-squares
        % solution. Gradients can therefore be sampled on a grid with a
        % spacing of sampleStep pixels in x and y. The sampling is done
        % after smoothing. It is usually reasonable to set sampleStep
        % roughly equal to sigmaXY.
        sampleStep

        % Smoothing constants for log-polar sampling
        %
        % Either a non-negative real scalar or a 2-element row vector of
        % non-negative real numbers; default 0. If a vector, the first
        % element is the parameter for the Gaussian mask used to smooth the
        % log-polar sampled image across rings (i.e. for a mask oriented
        % along the ring axis); the second element is the parameter for
        % smoothing across wedges. A scalar sets both parameters.
        sigmaRW
     
        % Centre of log-polar sampling pattern
        %
        % A 2-element positive real row vector of the form [X0 Y0], or the
        % empty matrix. X0 specifies the x (column) coordinate of the
        % centre of the log-polar sampling grid. Y0 specifies the y (row)
        % coordinate of the centre. If empty (the default) the centre of
        % the image is used.
        %
        % See also logsample
        logCentre
        
        % Inner radius of log-polar sampling pattern
        %
        % A positive integer; default 1. The radius of the smallest ring of
        % the log-polar sampling grid.
        %
        % See also logsample
        logRmin

        % Outer radius of log-polar sampling pattern
        %
        % A positive integer or the empty matrix, specifying the radius of
        % the largest ring of the log-polar sampling grid. If empty (the
        % default) and maxIter is 1, the largest value that keeps the outer
        % ring entirely within the image is used; if maxIter is greater
        % than 1, defaults to 75% of this.
        %
        % See also logsample
        logRmax

        % Number of wedges in the log-polar sampling pattern
        %
        % A positive integer or the empty matrix, specifying the number of
        % wedges in the log-polar sampling grid. If empty (the default)
        % round(logRmax) is used; this gives an arbitrary sample spacing of
        % 2*pi pixels in the outer ring.
        %
        % The number of rings is chosen automatically to satisfy the
        % "circular samples" constraint.
        %
        % See also logsample
        logWedges

        % Maximum number of iterations (log-polar sampling only)
        %
        % A positive integer; default 1. If greater than 1 the log-polar
        % grid applied to image2 is moved to track the current estimated
        % optic flow for the grid centre, and the flow is then reestimated.
        % This can effectively null out the flow at the centre of the grid,
        % allowing more accurate estimation. The flow parameters returned
        % are corrected for the grid motion.
        %
        % Iteration stops when the number of iterations reaches maxIter, or
        % when the residual velocity drops below velTol, or when the outer
        % ring of the log-polar grid would go outside the image.
        %
        % If regionOfInterest has been specified, the ROI for the second
        % image moves with the log-polar grid during tracking.
        maxIter

        % Speed estimation tolerance for tracking (log-polar sampling only)
        %
        % A positive real scalar; default 1. If tracking is in effect,
        % iteration will halt when the magnitude of the residual flow at
        % the centre of the grid has dropped below velTol. The residual
        % flow is the estimated true flow minus the offset of the log-polar
        % grid between IM1 and IM2.
          velTol
        
        % The initial image
        %
        % Must be set by assignment or as an argument to the constructor
        % before findFlow can be called. Must be set to a 2D array.
        image1
        
        % The image after motion has taken place
        %
        % Must be set by assignment or as an argument to the constructor
        % before findFlow can be called. Must be set to a 2D array. Must
        % have the same size as image1.
        image2
        
    end
    
    properties (Dependent, SetAccess = private)
        
        % The affine flow estimate as a structure
        %
        % The estimated parameters. The structure has fields vx0, vy0, d,
        % r, s1 and s2. For the meanings of the parameters, see the
        % introductory section of the help for the class.
        flowStruct
        
        % The affine flow estimate as a matrix
        %
        % The estimated parameters, stored as a matrix M such that the
        % estimated flow at a position (x, y) is given by [x y 1]*M.
        flowMatrix
        
        % The estimate of the image deformation matrix
        %
        % The estimated parameters, stored as a matrix W such that if (x,
        % y) is a position in image1, the estimate of the corresponding
        % position in image2 is given by [x y 1] * W.
        flowWarp
        
        % Tracking information (log-polar sampling only)
        %
        % The effect of tracking during the last call to findFlow is
        % reported in the trackInfo structure, whose fields are:
        %
        %   numIter - the number of iterations that took place;
        %   lastCentre - a 2-element row vector containing the x and y
        %       coordinates of the centre of the final sampling pattern
        %       applied to image2;
        %   lastspeed, the final residual flow speed at the log-polar
        %       centre after tracking has been applied.
        trackInfo
        
    end
    
    properties (Access = private)
        
        % Shadows of dependent variables
        
        im1 = []
        im2 = []
        rectmeth = true
        roi = []
        region = []
        sigmaxy = 0
        sampstep = 1
        sigmarw = [0 0]
        x0 = []
        y0 = []
        rmin = 1
        rmax = []
        nw = []
        maxiter = 1
        veltol = 1
        
        % Intermediate results, held to possibly reduce later computation
        
        smth1 = []
        smth2 = []
        roireg = []
        roitrimmed = []
        
        logsmth1 = []
        logsmth2 = []
        logreg = []
        logroi = []
        
        % Results of main computation
        
        flow = []
        resOK = false
        iterations = []
        lastx0 = []
        lasty0 = []
        lastspeed = []
    end
    
    
    methods         % constructor and get/set methods
        
        function a = affine_flow(varargin)
            % Constructor for affine flow objects
            %
            % af = affine_flow('prop1', val1, 'prop2', val2, ...)
            % constructs an affine flow object, setting property prop1 to
            % value val1 etc. See the class information for details.
            % Properties not set are given default values; to see these
            % execute disp(affine_flow) and
            % disp(affine_flow('sampleMethod', 'logpolar')).
            %
            % See also: affine_flowdemo
            a = setProps(a, varargin{:});
        end
        
        % Set/get methods with argument checking and state maintenance.
        % Care is needed to ensure that intermediate results are correctly
        % (but not unnecessarily) invalidated if properties are changed.
        
        function a = set.image1(a, im)
            % Could validate image1, but don't for speed
            a.im1 = im;
            a.smth1 = [];
            a.logsmth1 = [];
            a.resOK = false;
        end
        function im = get.image1(a)
            im = a.im1;
        end
        
        function a = set.image2(a, im)
            % Could validate image2, but don't for speed
            a.im2 = im;
            a.smth2 = [];
            a.logsmth2 = [];
            a.resOK = false;
        end
        function im = get.image2(a)
            im = a.im2;
        end
        
        function a = set.sampleMethod(a, s)
            if ~ismember(s, {'rectilinear' 'logpolar'})
                error('affine_flow:badmethod', ...
                    'Method must be rectilinear or logpolar');
            end
            a.rectmeth = isequal(s, 'rectilinear');
            a.resOK = false;
        end
        function s = get.sampleMethod(a)
            if a.rectmeth
                s = 'rectilinear';
            else
                s = 'logpolar';
            end
        end
        
        function a = set.regionOfInterest(a, roi)
            if ~isempty(roi)
                validateattributes(roi, {'double' 'logical'}, {'binary'});
            end
            a.roi = roi; a.roireg = []; a.logroi = [];
            % rect preprocessing needs to be redone because bounding box
            % may have changed, but lp preprocessing doesn't
            a.smth1 = []; a.smth2 = [];
            a.region = [];
            a.resOK = false;
        end
        function r = get.regionOfInterest(a)
            r = a.roi;
        end
        
        function a = set.rectRegion(a, reg)
            if ~isempty(reg)
                validateattributes(reg, {'double'}, {'integer' 'size' [1 4]});
            end
            a.roi = []; a.roireg = []; a.logroi = [];
            a.smth1 = []; a.smth2 = [];
            a.region = reg;
            a.resOK = false;
        end
        function r = get.rectRegion(a)
            a.checkrelevance('rectRegion', true);
            r = a.region;
        end
        
        function a = set.sigmaXY(a, sigma)
            validateattributes(sigma, {'double'}, ...
                {'nonnegative' 'real' 'finite' 'scalar'});
            a.sigmaxy = sigma;
            a.smth1 = []; a.smth2 = [];
            a.resOK = false;
        end
        function s = get.sigmaXY(a)
            a.checkrelevance('sigmaXY', true);
            s = a.sigmaxy;
        end
        
        function a = set.sampleStep(a, step)
            validateattributes(step, {'double'}, ...
                {'positive' 'integer' 'scalar'});
            a.sampstep = step;
            a.resOK = false;
        end
        function s = get.sampleStep(a)
            a.checkrelevance('sampleStep', true);
            s = a.sampstep;
        end
        
        function a = set.sigmaRW(a, sigma)
            validateattributes(sigma, {'double'}, ...
                {'nonnegative' 'real' 'finite'});
            if isscalar(sigma)
                a.sigmarw = [sigma sigma];
            else
                if ~isequal(size(sigma), [1 2])
                    error('affine_flow:badsigmaRW', ...
                        'sigmaRW must be scalar or 1x2 vector');
                end
                a.sigmarw = sigma;
            end
            a.logsmth1 = []; a.logsmth2 = []; a.logroi = [];
            a.resOK = false;
        end
        function s = get.sigmaRW(a)
            a.checkrelevance('sigmaRW', false);
            s = a.sigmarw;
        end
        
        function a = set.logCentre(a, c)
            if isempty(c)
                a.x0 = [];
                a.y0 = [];
            else
                validateattributes(c, {'double'}, ...
                    {'real' 'finite' 'scalar' 'size' [1 2]});
                a.x0 = c(1);
                a.y0 = c(2);
            end
            a.logsmth1 = []; a.logsmth2 = []; a.logroi = [];
            a.resOK = false;
        end
        function c = get.logCentre(a)
            a.checkrelevance('logCentre', false);
            c = [a.x0 a.y0];
        end
        
        function a = set.logRmin(a, r)
            validateattributes(r, {'double'}, ...
                {'positive' 'real' 'finite' 'scalar'});
            a.rmin = r;
            a.logsmth1 = []; a.logsmth2 = []; a.logroi = [];
            a.resOK = false;
        end
        function r = get.logRmin(a)
            a.checkrelevance('logRmin', false);
            r = a.rmin;
        end
        
        function a = set.logRmax(a, r)
            if ~isempty(r)
                validateattributes(r, {'double'}, ...
                    {'positive' 'real' 'finite' 'scalar'});
            end
            a.rmax = r;
            a.logsmth1 = []; a.logsmth2 = []; a.logroi = [];
            a.resOK = false;
        end
        function r = get.logRmax(a)
            a.checkrelevance('logRmax', false);
            r = a.rmax;
        end
        
        function a = set.logWedges(a, n)
            validateattributes(n, {'double'}, {'positive' 'integer' 'scalar'});
            a.nw = n;
            a.logsmth1 = []; a.logsmth2 = []; a.logroi = [];
            a.resOK = false;
        end
        function n = get.logWedges(a)
            a.checkrelevance('logWedges', false);
            n = a.nw;
        end
        
        function a = set.maxIter(a, n)
            validateattributes(n, {'double'}, {'positive' 'integer' 'scalar'});
            a.maxiter = n;
            a.resOK = false;
        end
        function n = get.maxIter(a)
            a.checkrelevance('maxIter', false);
            n = a.maxiter;
        end
        
        function a = set.velTol(a, v)
            validateattributes(v, {'double'}, ...
                {'positive' 'real' 'finite' 'scalar'});
            a.veltol = v;
            a.resOK = false;
        end
        function v = get.velTol(a)
            a.checkrelevance('velTol', false);
            v = a.veltol;
        end
        
        function f = get.flowStruct(a)
            if a.resOK
                f = a.flow;
            else
                error('affine_flow:noresult', 'No flow result available');
            end
        end
        
        function m = get.flowMatrix(a)
            m = affine_flow.matrix(a.flowStruct);
        end
        
        function w = get.flowWarp(a)
            w = affine_flow.warp(a.flowStruct);
        end
        
        function t = get.trackInfo(a)
            if a.resOK && ~a.rectmeth
                t.lastCentre = [a.lastx0 a.lasty0];
                t.lastSpeed = a.lastspeed;
                t.numIter = a.iterations;
            else
                error('affine_flow:notrack', 'No tracking result available');
            end
        end
            
    end
    
    
    methods         % public computational methods
        
        function a = findFlow(a)
            % Computes affine flow between image1 and image2
            %
            % Estimates flow using extended Lucas-Kanade method (see
            % algorithm description). Results are accessed via the
            % flowStruct/flowMatrix/flowWarp and trackInfo properties.
            
            if a.rectmeth
                a = a.affine_flowxy;
            else
                a = a.affine_flowrw;
            end
        end
        
        function a = advance(a)
            % advance advances image2 to image1.
            %
            % This method is for computation of successive sets of
            % flow fields in image sequences, by moving image2 to
            % image1 and updating the feature positions to the new
            % positions in image2. Its main effect is thus
            %
            %   af.image1 = af.image2;
            %
            % However, the method is more efficient than the assignment
            % would be. It could be used like this, assuming we have a cell
            % array of images and we want a cell array of flow field
            % results:
            %
            %   af = affine_flow('image1', images{1});
            %   for k = 2:length(images)
            %       af.image2 = images{k};
            %       af = af.findFlow;
            %       flowfields{i} = af.flowStruct;
            %       af = af.advance;
            %   end
            %
            % If log-polar sampling is used and tracking is in effect, the
            % initial sampling centre also moves to the final centre after
            % tracking, in effect:
            %
            %   af.logCentre = af.trackInfo.lastCentre;
            
            a.im1 = a.im2;
            a.smth1 = a.smth2;
            a.logsmth1 = a.logsmth2;
            a.x0 = a.lastx0;
            a.y0 = a.lasty0;
            a.resOK = false;
        end
        
    end
    
    
    methods (Access = private)
        
        function a = affine_flowxy(a)
            % Main computational function for rectilinear sampling.
            
            roil = a.roi;
            regionl = a.region;
            
            if ~isempty(roil)
                if isempty(a.roireg)
                    % Trim all arrays to the bounding box to avoid
                    % unnecessary work. Convert to double for regionprops
                    % as region may be unconnected
                    bbox = regionprops(double(a.roi), 'BoundingBox');
                    bbox = bbox.BoundingBox;
                    regionl = [ceil(bbox(2)) floor(bbox(2)+bbox(4)) ...
                        ceil(bbox(1)) floor(bbox(1)+bbox(3))];
                    roil = logical(a.getreg(a.roi, regionl));
                    a.roireg = regionl;
                    a.roitrimmed = roil;
                else
                    regionl = a.roireg;
                    roil = a.roitrimmed;
                end
            elseif isempty(regionl)
                regionl = gradients_xyt(a.im1, a.im2, a.sigmaxy, ...
                    'region', false);
            end
            
            % Image gradients and x,y coordinate arrays
            if isempty(a.smth1)
                a.smth1 = gsmooth2(a.im1, a.sigmaxy, ...
                    regionl + [-1 1 -1 1], false);
            end
            if isempty(a.smth2)
                a.smth2 = gsmooth2(a.im2, a.sigmaxy, ...
                    regionl + [-1 1 -1 1], false);
            end
            [xg, yg, tg] = gradients_xyt(a.smth1, a.smth2, 0, [], false);
            
            % Use origin at centre of image for now - more stable
            [x, y, cx, cy] = a.centredxy(xg);
            
            % Cut down typing by assembling the arguments
            flowargs = {x, y, xg, yg, tg};
            
            % Resample on bigger grid if required to cut down least-squares
            % solver effort
            if a.sampstep > 1
                steps = [a.sampstep a.sampstep];
                flowargs = cellfun(@(z) {affine_flow.sample(z, steps)}, ...
                    flowargs);
            end
            
            % Get values in the ROI if there is one, otherwise just convert
            % to column vectors
            if ~isempty(roil)
                if a.sampstep > 1
                    roil = affine_flow.sample(roil, steps);
                end
                flowargs = cellfun(@(z) {z(roil)}, flowargs);
            else
                flowargs = cellfun(@(z) {z(:)}, flowargs);
            end
            
            % Assemble and solve the least squares system
            flowl = affine_flow.solvexy(flowargs{:});
            
            % Move origin back to origin of image coordinates
            a.flow = affine_flow.shift(flowl, ...
                1-regionl(3)-cx, 1-regionl(1)-cy);
            a.resOK = true;
            
        end
        
        
        function a = affine_flowrw(a)
            % Main computational function for log-polar sampling.
            
            centre = (size(a.im1)+1)/2;
            xc = a.x0;
            if isempty(xc)
                xc = centre(2);
            end
            yc = a.y0;
            if isempty(yc)
                yc = centre(1);
            end
            rmx = a.rmax;
            if isempty(rmx)
                % keep samples just within image
                rmx = min([centre-1 size(a.im1)-centre]);
                % but if we are iterating, allow room for movement
                if a.maxIter > 1
                    rmx = 0.75 * rmx;   % arbitrary
                end
            end
            nwed = a.nw;
            if isempty(nwed)
                nwed = round(rmx);     % quite arbitrary
            end
            
            [maxy, maxx] = size(a.im1);
            
            % Sample first image onto log-polar grid
            if isempty(a.logsmth1)
                logim1 = logsample(a.im1, a.rmin, rmx, xc, yc, [], nwed);
                [a.logsmth1, regions] = gsmooth2(logim1, a.sigmarw, [], ...
                    [false true]);
                a.logreg = regions + [0 0 1 -1];
            end
            
            [r, w] = affine_flow.logcoords(a.logreg);
            
            % and roi
            roil = a.roi;
            if isempty(a.logroi) && ~isempty(roil)
                % imtransform (and so logsample) is OK with logical images
                logroil = logsample(logical(roil), a.rmin, rmx, ...
                    xc, yc, [], nwed);
                a.logroi = affine_flow.getreg(logroil, a.logreg);
            end
            
            xc2 = xc;
            yc2 = yc;
            iter = 0;
            
            while true      % iterate to track motion
                iter = iter + 1;
                
                % Sample second image onto current centre for second grid
                if isempty(a.logsmth2)
                    logim2 = logsample(a.im2, a.rmin, rmx, xc2, yc2, [], nwed);
                    a.logsmth2 = gsmooth2(logim2, a.sigmarw, [], ...
                        [false true]);
                end
                
                % Get gradients
                [gr, gw, gt] = gradients_xyt(a.logsmth1, a.logsmth2, ...
                    0, [], [false true]);
                
                % Cut down typing by assembling the arguments
                flowargs = {r, w, gr, gw, gt};
                
                % Get values in the ROI if there is one, otherwise just
                % convert to column vectors
                if ~isempty(roil)
                    flowargs = cellfun(@(z) {z(a.logroi)}, flowargs);
                else
                    flowargs = cellfun(@(z) {z(:)}, flowargs);
                end
                
                % Build and solve the least-squares system
                flowl = affine_flow.solverw(a.rmin, nwed, flowargs{:});
                
                % Prepare for tracking if required. save most recent sample
                % centre
                a.lastx0 = xc2;
                a.lasty0 = yc2;
                % update centre position to follow flow
                xc2 = xc2 + flowl.vx0;
                yc2 = yc2 + flowl.vy0;
                % test whether finished
                speed = norm([flowl.vx0 flowl.vy0]);
                if iter >= a.maxiter || speed <= a.veltol || ...
                        min([xc2-1 yc2-1 maxx-xc2 maxy-yc2]) < rmx
                    break;
                end
                
                a.logsmth2 = [];
            end
            
            % Correct centre velocity for sample centre shift
            flowl.vx0 = xc2 - xc;
            flowl.vy0 = yc2 - yc;
            % Move origin back to origin of image coordinates
            a.flow = affine_flow.shift(flowl, -xc, -yc);
            a.resOK = true;
            
            % Record extra information
            a.iterations = iter;
            a.lastspeed = speed;
            
        end
        
        function checkrelevance(a, propname, rectprop)
            if rectprop ~= a.rectmeth
                error('affine_flow:badfieldrequest', ...
                    ['Property ' propname ...
                    ' not used for current sampling method']);
            end
        end
        
    end
    
    
    methods (Static, Access = private)
        
        function [x, y, cx, cy] = centredxy(im)
            % Makes x and y arrays centred on the centre of the image
            [s1, s2] = size(im);
            cy = (s1-1)/2;
            cx = (s2-1)/2;
            [x, y] = meshgrid(linspace(-cx, cx, s2), linspace(-cy, cy, s1));
        end
        
        
        function [r, w] = logcoords(reg)
            % makes r and w index arrays for the specified region of the
            % log-sampled image
            [r, w] = meshgrid(reg(3)-1:reg(4)-1, reg(1)-1:reg(2)-1);
        end
        
        
        function arr = getreg(arr, reg)
            % trim array to specified region
            arr = arr(reg(1):reg(2), reg(3):reg(4));
        end
        
        
        function arr = sample(arr, steps)
            % sample every sx'th column and sy'th row, centering the grid
            sz = size(arr);
            nm1 = floor((sz-1)./steps);
            starts = floor((sz - (steps.*nm1 + 1))/2) + 1;
            arr = arr(starts(2):steps(2):end, starts(1):steps(1):end);
        end
        
        
        function affineparams = solvexy(x, y, gx, gy, gt)
            % AFFINE.SOLVEXY first order optic flow parameters -rectilinear
            %   AFFINEPARAMS = AFFINE.SOLVEXY(X, Y, GXS, GYS, GTS) returns
            %   the affine flow parameters vx0, vy0, d, r, s1, s2, as a
            %   struct with those f ields. The inputs are the x and y
            %   sample positions and the x, y and t gradients, all as
            %   column vectors.
            
            xgx = x .* gx;
            ygx = y .* gx;
            xgy = x .* gy;
            ygy = y .* gy;
            A = [gx, gy, xgx+ygy, xgy-ygx, xgx-ygy, ygx+xgy];
            
            a = -(A \ gt);
            
            affineparams = struct('vx0', a(1), 'vy0', a(2), ...
                'd', a(3), 'r', a(4), 's1', a(5), 's2', a(6));
        end
        
        
        function affineparams = solverw(rmin, nw, r, w, gr, gw, gt)
            % AFFINE.SOLVERW first order optic flow parameters - log-polar
            %   AFFINEPARAMS = AFFINE.SOLVERW(RMIN, NW, R, W, GR, GW, GT)
            %   returns the affine flow parameters vx0, vy0, d, r, s1, s2,
            %   as a struct with those fields. The inputs are the
            %   parameters of the log-polar grid and the gradients measured
            %   on it.
            
            kinv = 2*pi/nw;
            rhoinv = exp(-kinv*r)/rmin;
            theta = w * kinv;
            c = cos(theta);
            s = sin(theta);
            grc = gr .* c;
            grs = gr .* s;
            gwc = gw .* c;
            gws = gw .* s;
            
            gx = (grc - gws) .* rhoinv;
            gy = (grs + gwc) .* rhoinv;
            gs1 = grc.*c - grs.*s - 2*gwc.*s;
            gs2 = 2*grc.*s + gwc.*c - gws.*s;
            
            A = [gx, gy, gr, gw, gs1, gs2];
            
            a = -kinv * (A \ gt);
            
            affineparams = struct('vx0', a(1), 'vy0', a(2), ...
                'd', a(3), 'r', a(4), 's1', a(5), 's2', a(6));
        end
        
    end
    
    
    methods (Static)
        
        function m = matrix(f)
            %AFFINE_FLOW.MATRIX converts affine flow structure to matrix
            %   M = AFFINE_FLOW.MATRIX(F) takes an affine flow structure as
            %   returned by AFFINE_FLOW and returns a matrix M such that if
            %   P is a row vector with components [X, Y, 1] representing a
            %   position, P*M is a vector representing the optic flow at
            %   that position.
            m = [f.d+f.s1,  f.s2+f.r;
                f.s2-f.r,  f.d-f.s1;
                f.vx0,     f.vy0];
        end
        
        function w = warp(f)
            %AFFINE_FLOW.WARP converts a flow to a warp
            %   W = AFFINE_FLOW.WARP(F) takes either a flow structure as
            %   returned by AFFINE_FLOW or a flow matrix as returned by
            %   AFFINE_FLOW2MATRIX and returns a matrix W such that if P is
            %   a row vector with components [X, Y, 1] representing a
            %   position, then P*M is a vector representing the new
            %   position of the vector after one frame of the flow has
            %   occured.
            %
            %   Uses a very simple approximation!
            if isstruct(f)
                f = affine_flow.matrix(f);
            end
            w = f + [eye(2); 0 0];
        end
        
        function f = shift(f, x0, y0)
            %AFFINE_FLOW.SHIFT shifts the origin of affine flow
            %   F = AFFINE_FLOW.SHIFT(F, X0, Y0) returns affine flow
            %   parameters as returned by AFFINE_FLOW, with the origin
            %   shifted to X0, Y0 relative to the current origin.
            f.vx0 = f.vx0 + x0*(f.d+f.s1) + y0*(f.s2-f.r);
            f.vy0 = f.vy0 + x0*(f.s2+f.r) + y0*(f.d-f.s1);
        end
        
    end
    
end

Contact us