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
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 General Public License as published by 0035 % the Free Software Foundation; either version 2 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 General Public License for more details. 0042 % 0043 % You can obtain a copy of the GNU General Public License from 0044 % http://www.gnu.org/copyleft/gpl.html or by writing to 0045 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 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