New Templatematch

demobatura.m

@@ -17,9 +17,9 @@
 
 A=imread(fullfile(datafolder,'batura_2001.tif')); 
 B=imread(fullfile(datafolder,'batura_2002.tif')); %normally you would use geotiffread
-x=(0:size(A,2)-1)*15+451357.50; %if you have mapping toolbox then use pixcenters here. 
+x=(0:size(A,2)-1)*15+451357.50; %if you have mapping toolbox then you would use pixcenters here. 
 y=(0:size(A,1)-1)*15+4060432.50;
-dx=15;%m/pixel
+deltax=15;%m/pixel
 
 
 %make regular grid of points to track:
@@ -31,12 +31,10 @@
      1061 1168;663 718;456 686;25 877;28 627;407 465];
 
 ix=find(inpolygon(pu,pv,roi(:,1),roi(:,2)));
-uvA=[pu(ix) pv(ix)]; %these are the points we want to track.
+pu=pu(ix); pv=pv(ix);
 
-whtemplate=10; %template/chip size
-whsearch=30; %search window size
-super=1; %No super sampling 
-[dxy,C]=templatematch(A,B,uvA,whtemplate,whsearch,super,[0 0],{'2001' '2002'});
+t
+[du,dv,C,Cnoise,pu,pv]=templatematch(A,B,pu,pv,'showprogress',{'2001' '2002'});
 close all
 
 %visualize the results
@@ -46,17 +44,13 @@
 axis equal off tight ij
 hold on
 
-signal2noise=C(:,1)./C(:,2);
-keep=(signal2noise>2.3)&(C(:,1)>.7);
-V=dxy*dx; 
-Vn=sqrt(sum(V.^2,2));
-scatter(uvA(keep,1),uvA(keep,2),200,Vn(keep),'.') %colors speed
-quiver(uvA(keep,1),uvA(keep,2),V(keep,1)./Vn(keep),V(keep,2)./Vn(keep),0.2,'k') %arrows show direction. 
+signal2noise=C./Cnoise;
+keep=(signal2noise>2.3)&(C>.65);
+V=(du+dv*1i)*deltax; %m/yr
+V(~keep)=nan;
+Vn=abs(V);
+scatter(pu(keep),pv(keep),300,Vn(keep),'.') %colors speed
+quiver(pu,pv,real(V)./Vn,imag(V)./Vn,0.2,'k') %arrows show direction. 
 colormap jet
 caxis([0 200])
 colorbar('southoutside');
-% if ~verLessThan('matlab', '8.4.0')
-%     set(hcb,'limits',caxis) %workaround for a critical colorbar bug in Matlab 2014b preview... TODO: check if bug present in final release
-% end
-
-

demobindschadler.m

@@ -13,25 +13,23 @@
 B=imread(fullfile(datafolder,'conv_89.png'));
 
 %regular grid of points to track:
-[pu,pv]=meshgrid(10:10:size(A,2),10:10:size(A,1));
 
-uvA=[pu(:) pv(:)]; %these are the points we want to track.
+tic
 
-whtemplate=10; %template size half-width
-whsearch=30;   %search template half-width
-super=1;
-[dxy,C]=templatematch(A,B,uvA,whtemplate,whsearch,super,[0 0],true);
+[du,dv,C,Cnoise,pu,pv]=templatematch(A,B);
+toc
 
 close all
 image(repmat(A,[1 1 3]),'CDataMapping','scaled') %the cdatamapping is a workaround for a bug in R2014+)
 axis equal off tight ij
 hold on
-signal2noise=C(:,1)./C(:,2);
-keep=(signal2noise>2)&(C(:,1)>.5);
-V=dxy*28.5/2; %m/yr
-Vn=sqrt(sum(V.^2,2));
-scatter(uvA(keep,1),uvA(keep,2),300,Vn(keep),'.') %colors speed
-quiver(uvA(keep,1),uvA(keep,2),V(keep,1)./Vn(keep),V(keep,2)./Vn(keep),0.2,'k') %arrows show direction. 
+signal2noise=C./Cnoise;
+keep=(signal2noise>2)&(C>.5);
+V=(du+dv*1i)*28.5/2; %m/yr - Using imaginary values to indicate direction (for convenience).
+V(~keep)=nan;
+Vn=abs(V);
+scatter(pu(keep),pv(keep),300,Vn(keep),'.') %colors speed
+quiver(pu,pv,real(V)./Vn,imag(V)./Vn,0.2,'k') %arrows show direction. 
 colormap jet
 caxis([0 400])
 colorbar('southoutside');

demoengabreen.m

@@ -79,24 +79,23 @@
 
 %First get an approximate estimate of the image shift using a single large
 %template
-points=[3000, 995];
-xyoffset=templatematch(A,B,points,200,260,0.5,[0 0],false,'PC') %Note PC=PhaseCorrelation is still experimental. 
+[duoffset,dvoffset]=templatematch(A,B,3000,995,'templatewidth',261,'searchwidth',400,'supersample',0.5,'showprogress',false,'method','PC') %Note PC=PhaseCorrelation is still experimental. 
 
 %Get a whole bunch of image shift estimates using a grid of probe points.
 %Having multiple shift estimates will allow us to determine camera
 %rotation.
-[pX,pY]=meshgrid(200:700:4000,100:400:1000);
-points=round([pX(:) pY(:)+pX(:)/10]); 
-[dxy,C]=templatematch(A,B,points,30,40,3,xyoffset,[idA idB]);
+[pu,pv]=meshgrid(200:700:4000,100:400:1000);
+pu=pu(:); pv=pv(:)+pu/10; 
+[du,dv,C]=templatematch(A,B,pu,pv,'templatewidth',61,'searchwidth',81,'supersample',3,'initialdu',duoffset,'initialdv',dvoffset,'showprogress',[idA idB]);
 
 
 %Determine camera rotation between A and B from the set of image
 %shifts.
 
 %find 3d coords consistent with the 2d pixel coords in points.
-xyz=camA.invproject(points);
+xyz=camA.invproject([pu pv]);
 %the projection of xyz has to match the shifted coords in points+dxy:
-[camB,rmse]=camA.optimizecam(xyz,points+dxy,'00000111000000000000'); %optimize 3 view direction angles to determine camera B. 
+[camB,rmse]=camA.optimizecam(xyz,[pu+du pv+dv],'00000111000000000000'); %optimize 3 view direction angles to determine camera B. 
 rmse
 
 
@@ -170,7 +169,7 @@
 xyzA=[X(keepers) Y(keepers) interp2(dem.X,dem.Y,dem.filled,X(keepers),Y(keepers))];
 [uvA,~,inframe]=camA.project(xyzA); %where would the candidate points be in image A
 xyzA=xyzA(inframe,:); %cull points outside the camera field of view.
-uvA=round(uvA(inframe,:)); %round because template match only works with integer pixel coords
+uvA=uvA(inframe,:); %round because template match only works with integer pixel coords
 uvA(end+1,:)=[2275 1342]; %add a non-glaciated point to test for residual camera movement (here a tunnel entrance)
 %Note xyzA no longer corresponds exactly to uvA because of the rounding.
 
@@ -181,15 +180,21 @@
 % ( i.e. accounting only for camera shake)
 camshake=camB.project(camA.invproject(uvA))-uvA;
 
-showprogress=[idA idB];
-wsearch=40; 
-wtemplate=10;
-super=5; %supersample the input images
+options=[];
+options.pu=uvA(:,1);
+options.pv=uvA(:,2);
+options.showprogress=[idA idB];
+options.searchwidth=81; 
+options.templatewidth=21;
+options.supersample=5; %supersample the input images
+options.initialdu=camshake(:,1);
+options.initialdv=camshake(:,2);
 
-[dxy,C]=templatematch(A,B,uvA,wtemplate,wsearch,super,camshake,showprogress); %myNCC is faster than NCC in my tests
+[du,dv,C,Cnoise,pu,pv]=templatematch(A,B,options); 
+uvA=[pu pv]; %the centers of the templates may have been rounded to nearest pixel. 
 
-uvB=uvA+dxy;
-signal2noise=C(:,1)./C(:,2);
+uvB=uvA+[du dv];
+signal2noise=C./Cnoise;
 
 %% Georeference tracked points
 % ... and calculate velocities
@@ -204,7 +209,7 @@
 axis equal xy off tight
 hold on
 Vn=sqrt(sum(V(:,1:2).^2,2));
-keep=signal2noise>2&C(:,1)>.7;
+keep=signal2noise>2&C>.7;
 scatter(xyzA(keep,1),xyzA(keep,2),100,Vn(keep),'.')
 quiver(xyzA(keep,1),xyzA(keep,2),V(keep,1)./Vn(keep),V(keep,2)./Vn(keep),.2,'k')
 caxis([0 1])