Home > voicebox > imagehomog.m

imagehomog

PURPOSE ^

IMAGEHOMOG Apply a homography transformation to an image with bilinear interpolation

SYNOPSIS ^

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

DESCRIPTION ^

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:

SOURCE CODE ^

0001 function [ih,xa,ya]=imagehomog(im,h,m,clip)
0002 %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: imagehomog.m 1641 2012-03-16 16:20:40Z 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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003