# v_imagehomog

## PURPOSE

V_IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation

## SYNOPSIS

function [ih,xa,ya]=v_imagehomog(im,h,m,clip)

## DESCRIPTION

    V_IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation
    Inputs: im(ny,nx,nc)  input image (uint8)
            h(3,3)        homography
            m             mode string
                             i  show image [default if no output arguments]
                             m  h uses matlab coordinates: (1,1) is top left [default]
                             c  h uses central coordinates: (0,0) is the centre = {(1+nx)/2,(1+ny)/2} in 'm'
                             k  clip to original image dimensions (equivalent to clip=1)
                             x  add blank rows and cols to extend image to the specified clipping rectangle
                             t  trim output image by re moving blank rows and columns
            clip(4)       output image bounding box [xmin xmax ymin ymax] (same coordinate system as h)
                             or [n] or [nx ny] to clip to a multiple of the original image [default clip=2]
     Outputs:
            ih(my,mx,nc)  output image (uint8)
            xa(mx)        x axis
            ya(my)        y axis

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_rectifyhomog V_RECTIFYHOMOG Apply rectifying homographies to an image set

## SOURCE CODE

    0001 function [ih,xa,ya]=v_imagehomog(im,h,m,clip)
    0002 %V_IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation
    0003 %Inputs: im(ny,nx,nc)  input image (uint8)
    0004 %        h(3,3)        homography
    0005 %        m             mode string
    0006 %                         i  show image [default if no output arguments]
    0007 %                         m  h uses matlab coordinates: (1,1) is top left [default]
    0008 %                         c  h uses central coordinates: (0,0) is the centre = {(1+nx)/2,(1+ny)/2} in 'm'
    0009 %                         k  clip to original image dimensions (equivalent to clip=1)
    0010 %                         x  add blank rows and cols to extend image to the specified clipping rectangle
    0011 %                         t  trim output image by re moving blank rows and columns
    0012 %        clip(4)       output image bounding box [xmin xmax ymin ymax] (same coordinate system as h)
    0013 %                         or [n] or [nx ny] to clip to a multiple of the original image [default clip=2]
    0014 % Outputs:
    0015 %        ih(my,mx,nc)  output image (uint8)
    0016 %        xa(mx)        x axis
    0017 %        ya(my)        y axis
    0018
    0019 % Bugs/Suggestions:
    0020 % (1) cope with non-uint8
    0021 % (2) cope with (a) multiple inputs, (b) multiple transformations
    0022 % (3) output a boundary mask as an alpha channel
    0023 % (4) do anti-aliasing along the boundary
    0024 % (5) check that origin shift is correct for central coordinates
    0025
    0026 %      Copyright (C) Mike Brookes 2010-2012
    0027 %      Version: $Id: v_imagehomog.m 10865 2018-09-21 17:22:45Z dmb $
    0028 %
    0029 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0030 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0031 %
    0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0033 %   This program is free software; you can redistribute it and/or modify
    0034 %   it under the terms of the GNU Lesser General Public License as published by
    0035 %   the Free Software Foundation; either version 3 of the License, or
    0036 %   (at your option) any later version.
    0037 %
    0038 %   This program is distributed in the hope that it will be useful,
    0039 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0040 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0041 %   GNU Lesser General Public License for more details.
    0042 %
    0043 %   You can obtain a copy of the GNU Lesser General Public License from
    0044 %   https://www.gnu.org/licenses/ .
    0045 %
    0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0047
    0048 maxby=1e7;  % maximum memory to use
    0049 [ny,nx,nc]=size(im);
    0050 if nargin<4 || ~numel(clip)
    0051     clip=2;
    0052 end
    0053 if nargin<3 || ~numel(m)
    0054     m='';
    0055 end
    0056 if nargin<2 || ~numel(h)
    0057     h=eye(3);
    0058 end
    0059
    0060 imr=reshape(im,nx*ny,nc); % vectorize the input image
    0061 t=eye(3);
    0062 if any(m=='c')   % convert homography and clipping box to matlab coordinates
    0063     t(7:8)=0.5+[nx ny]/2;  % shift origin to image centre
    0064     h=t*h/t;  % change homography so input and output use MATLAB coordinates
    0065     if numel(clip)==4
    0066         clip=clip+t([7 7 8 8]);  % make clipping use MATLAB coordinates as well
    0067     end
    0068 end
    0069 % determine clipping rectangle for output image
    0070 if any(m=='k')
    0071     clip=[1 1];
    0072 elseif numel(clip)==1
    0073     clip=clip(ones(1,2));
    0074 end
    0075 if numel(clip)==2
    0076     clip=[[1 nx]*clip(1)-(clip(1)-1)*(1+nx)/2 [1 ny]*clip(2)-(clip(2)-1)*(1+ny)/2];
    0077 end
    0078 clip=clip(:)';
    0079 clip(1:2:3)=floor(clip(1:2:3));
    0080 clip(2:2:4)=ceil(clip(2:2:4));
    0081 % now determine the image of the source image
    0082 bi=[1 1 nx nx; 1 ny ny 1; 1 1 1 1];
    0083 box=h*bi;
    0084 b3=box(3,:);
    0085 if any(b3<=0)
    0086     ib=find(b3>0);
    0087     nb=numel(ib);
    0088     bb=ones(3,nb+2);
    0089     if ~nb
    0090         error('image invisible');
    0091     end
    0092     bb(:,1:nb)=box(:,ib);
    0093     px=[4 1:3];
    0094     ib=find(b3.*b3(px)<=0); % find points at infinity
    0095     ip=px(ib);
    0096     af=b3(ip);
    0097     bf=b3(ib);
    0098     pof=[-1 0 1 0; 0 1 0 -1; 0 0 0 0]; % offsets to avoid exact infinity
    0099     w3=ones(3,1);
    0100     bb(:,nb+(1:2))=h*((bf(w3,:).*bi(:,ip)-af(w3,:).*bi(:,ib))./repmat(bf-af,3,1)+(pof(:,ib).*repmat(2*(b3(ib)>0)-1,3,1)));
    0101     box=bb;
    0102 end
    0103 box=box(1:2,:)./box([3 3],:);  % coordinates of mapped original image
    0104 box=[min(box(1,:)) max(box(1,:)) min(box(2,:)) max(box(2,:))];
    0105
    0106 box(1:2:3)=floor(max(clip(1:2:3),box(1:2:3)));  % no point in calculating non-existant points
    0107 box(2:2:4)=ceil(min(clip(2:2:4),box(2:2:4)));
    0108 g=inv(h);
    0109 mx=box(2)-box(1)+1; % number of columns in destination
    0110 my=box(4)-box(3)+1; % number of rows in destination
    0111
    0112 ih=zeros(my*mx,nc,'uint8');
    0113 ncol=max(1,min(mx,floor(maxby/(my*nc*8)))); % number of columns to do in a chunk
    0114 nloop=ceil(mx/ncol);
    0115 cmax=1+rem(mx-1,ncol); % final column of first iteration
    0116 jxinc=ncol*my;      % increment target indices each loop
    0117 ginc=g(:,1)*ncol; % increment transformed targets each loop
    0118 wj=ones(1,jxinc);   % repeat index for ginc
    0119 jx=1+cmax*my-jxinc:cmax*my;  % initial target indices (some might be negative)
    0120 gk=g*[reshape(repmat(cmax-ncol+box(1):cmax+box(1)-1,my,1),1,jxinc); repmat(box(3):box(4),1,ncol); ones(1,jxinc)];
    0121 gn=gk(1:2,:)./gk([3 3],:);   % normalize source coordinates
    0122 mn=[zeros(1,jxinc-cmax*my) ones(1,cmax*my)]; % mask for initial iteration
    0123 mn=mn & (gn(1,:)>-0.5 & gn(2,:)>-0.5 & gn(1,:)<nx+0.5 & gn(2,:)<ny+0.5); % mask valid pixels
    0124 w3=ones(nc,1);
    0125 for i=1:nloop
    0126     fn=max(floor(gn(:,mn)),1);
    0127     fn1=min(max(fn(1,:)',1),nx-1);
    0128     fn2=min(max(fn(2,:)',1),ny-1);
    0129     dn=gn(:,mn)-[fn1 fn2]';
    0130     dn1=min(max(dn(1,:)',0),1);
    0131     dn2=min(max(dn(2,:)',0),1);
    0132     dn1c=1-dn1;
    0133     dn2c=1-dn2;
    0134     ih(jx(mn),:)=uint8(dn1c(:,w3).*(dn2c(:,w3).*single(imr(fn2+ny*(fn1-1),:)) ...
    0135         +dn2(:,w3).*single(imr(fn2+1+ny*(fn1-1),:))) ...
    0136         +dn1(:,w3).*(dn2c(:,w3).*single(imr(fn2+ny*(fn1),:)) ...
    0137         +dn2(:,w3).*single(imr(fn2+1+ny*(fn1),:))));
    0138     jx=jx+jxinc;  % target indices
    0139     gk=gk+ginc(:,wj);
    0140     gn=gk(1:2,:)./gk([3 3],:);   % normalize source coordinates
    0141     mn=gn(1,:)>-0.5 & gn(2,:)>-0.5 & gn(1,:)<nx+0.5 & gn(2,:)<ny+0.5; % mask valid pixels
    0142 end
    0143 ih=reshape(ih,[my,mx,nc]);
    0144 if any(m=='x')          % extend blank area to specified clipping rectangle
    0145     ih=[zeros(box(3)-clip(3),clip(2)-clip(1)+1,nc,'uint8'); ...
    0146         zeros(my,box(1)-clip(1),nc,'uint8') ih zeros(my,clip(2)-box(2),nc,'uint8');
    0147         zeros(clip(4)-box(4),clip(2)-clip(1)+1,nc,'uint8')];
    0148     xa=(clip(1):clip(2))-t(7);
    0149     ya=(clip(3):clip(4))-t(8);
    0150 else
    0151     xa=(box(1):box(2))-t(7);
    0152     ya=(box(3):box(4))-t(8);
    0153 end
    0154 if any(m=='t')
    0155      ihs=sum(ih,3);
    0156      azx=all(~ihs,1);
    0157      ix1=find(~azx,1,'first');
    0158      if ~numel(ix1)
    0159          ih=[];
    0160          xa=[];
    0161          ya=[];
    0162      else
    0163          ix2=find(~azx,1,'last');
    0164          azx=all(~ihs,2);
    0165          iy1=find(~azx,1,'first');
    0166          iy2=find(~azx,1,'last');
    0167          ih=ih(iy1:iy2,ix1:ix2,:);
    0168          xa=xa(ix1:ix2);
    0169          ya=ya(iy1:iy2);
    0170      end
    0171 end
    0172 if ~nargout || any(m=='i')
    0173     imagesc(xa,ya,ih);
    0174     axis image
    0175 end
    0176