模板匹配算法

function[I_SSD,I_NCC,Idata]=template_matching(T,I,IdataIn)%TEMPLATE_MATCHINGisacpuefficientfunctionwhichcalculatesmatching%scoreimagesbetweentemplateandan(color)2Dor3Dimage.%Itcalculates:%-Thesumofsquareddifference(SSDBlockMatching),robusttemplate%matching.%-Thenormalizedcrosscorrelation(NCC),independentofillumination,%onlydependentontexture%Theusercancombinethetwoimages,togettemplatematchingwhich%worksrobustwithhisapplication.%BothmeasuresareimplementedusingFFTbasedcorrelation.%%[I_SSD,I_NCC,Idata]=template_matching(T,I,Idata)%%inputs,%T:ImageTemplate,canbegrayscaleorcolor2Dor3D.%I:Colorimage,canbegrayscaleorcolor2Dor3D.%(optional)%Idata:StorageoftemporaryvariablesfromtheimageI,toallow%fastersearchformultipletemplatesinthesameimage.%%outputs,%I_SSD:Thesumofsquareddifference2D/3Dimage.TheSSDsignis%reversedandnormalizedtorange[01]%I_NCC:Thenormalizedcrosscorrelation2D/3Dimage.Thevalues%rangebetween0and1%Idata:StorageoftemporaryvariablesfromtheimageI,toallow%fastersearchformultipletemplatesinthesameimage.%%Example2D,%%Findmaximumresponse%I=im2double(imread('lena.jpg'));%%TemplateofEyeLena%T=I(124:140,124:140,:);%%CalculateSSDandNCCbetweenTemplateandImage%[I_SSD,I_NCC]=template_matching(T,I);%%FindmaximumcorrespondenceinI_SDDimage%[x,y]=find(I_SSD==max(I_SSD(:)));%%Showresult%figure,%subplot(2,2,1),imshow(I);holdon;plot(y,x,'r*');title('Result')%subplot(2,2,2),imshow(T);title('Theeyetemplate');%subplot(2,2,3),imshow(I_SSD);title('SSDMatching');%subplot(2,2,4),imshow(I_NCC);title('Normalized-CC');%%Example3D,%%Makesomerandomdata%I=rand(50,60,50);%%Getasmallvolumeastemplate%T=I(20:30,20:30,20:30);%%CalculateSDDbetweentemplateandimage%I_SSD=template_matching(T,I);%%Findmaximumcorrespondence%[x,y,z]=ind2sub(size(I_SSD),find(I_SSD==max(I_SSD(:))));%disp(x);%disp(y);%disp(z);%%FunctioniswrittenbyD.KroonUniversityofTwente(February2011)if(nargin3),IdataIn=[];end%ConvertimagestodoubleT=double(T);I=double(I);if(size(T,3)==3)%ColorImagedetected[I_SSD,I_NCC,Idata]=template_matching_color(T,I,IdataIn);else%Grayscaleimageor3Dvolume[I_SSD,I_NCC,Idata]=template_matching_gray(T,I,IdataIn);endfunction[I_SSD,I_NCC,Idata]=template_matching_color(T,I,IdataIn)if(isempty(IdataIn)),IdataIn.r=[];IdataIn.g=[];IdataIn.b=[];end%Splitecolorimage,anddotemplatematchingonR,GandBimage[I_SSD_R,I_NCC_R,Idata.r]=template_matching_gray(T(:,:,1),I(:,:,1),IdataIn.r);[I_SSD_G,I_NCC_G,Idata.g]=template_matching_gray(T(:,:,2),I(:,:,2),IdataIn.g);[I_SSD_B,I_NCC_B,Idata.b]=template_matching_gray(T(:,:,3),I(:,:,3),IdataIn.b);%CombinetheresultsI_SSD=(I_SSD_R+I_SSD_G+I_SSD_B)/3;I_NCC=(I_NCC_R+I_NCC_G+I_NCC_B)/3;function[I_SSD,I_NCC,Idata]=template_matching_gray(T,I,IdataIn)%Calculatecorrelationoutputsize=inputsize+paddingtemplateT_size=size(T);I_size=size(I);outsize=I_size+T_size-1;%calculatecorrelationinfrequencydomainif(length(T_size)==2)FT=fft2(rot90(T,2),outsize(1),outsize(2));if(isempty(IdataIn))Idata.FI=fft2(I,outsize(1),outsize(2));elseIdata.FI=IdataIn.FI;endIcorr=real(ifft2(Idata.FI.*FT));elseFT=fftn(rot90_3D(T),outsize);FI=fftn(I,outsize);Icorr=real(ifftn(FI.*FT));end%CalculateLocalQuadraticsumofImageandTemplateif(isempty(IdataIn))Idata.LocalQSumI=local_sum(I.*I,T_size);elseIdata.LocalQSumI=IdataIn.LocalQSumI;endQSumT=sum(T(:).^2);%SSDbetweentemplateandimageI_SSD=Idata.LocalQSumI+QSumT-2*Icorr;%Normalizetorange0..1I_SSD=I_SSD-min(I_SSD(:));I_SSD=1-(I_SSD./max(I_SSD(:)));%RemovepaddingI_SSD=unpadarray(I_SSD,size(I));if(nargout1)%NormalizedcrosscorrelationSTDif(isempty(IdataIn))Idata.LocalSumI=local_sum(I,T_size);elseIdata.LocalSumI=IdataIn.LocalSumI;end%Standarddeviationif(isempty(IdataIn))Idata.stdI=sqrt(max(Idata.LocalQSumI-(Idata.LocalSumI.^2)/numel(T),0));elseIdata.stdI=IdataIn.stdI;endstdT=sqrt(numel(T)-1)*std(T(:));%MeancompensationmeanIT=Idata.LocalSumI*sum(T(:))/numel(T);I_NCC=0.5+(Icorr-meanIT)./(2*stdT*max(Idata.stdI,stdT/1e5));%RemovepaddingI_NCC=unpadarray(I_NCC,size(I));endfunctionT=rot90_3D(T)T=flipdim(flipdim(flipdim(T,1),2),3);functionB=unpadarray(A,Bsize)Bstart=ceil((size(A)-Bsize)/2)+1;Bend=Bstart+Bsize-1;if(ndims(A)==2)B=A(Bstart(1):Bend(1),Bstart(2):Bend(2));elseif(ndims(A)==3)B=A(Bstart(1):Bend(1),Bstart(2):Bend(2),Bstart(3):Bend(3));endfunctionlocal_sum_I=local_sum(I,T_size)%AddpaddingtotheimageB=padarray(I,T_size);%Calculateforeachpixelthesumoftheregionaroundit,%withtheregionthesizeofthetemplate.if(length(T_size)==2)%2Dlocalsums=cumsum(B,1);c=s(1+T_size(1):end-1,:)-s(1:end-T_size(1)-1,:);s=cumsum(c,2);local_sum_I=s(:,1+T_size(2):end-1)-s(:,1:end-T_size(2)-1);else%3DLocalsums=cumsum(B,1);c=s(1+T_size(1):end-1,:,:)-s(1:end-T_size(1)-1,:,:);s=cumsum(c,2);c=s(:,1+T_size(2):end-1,:)-s(:,1:end-T_size(2)-1,:);s=cumsum(c,3);local_sum_I=s(:,:,1+T_size(3):end-1)-s(:,:,1:end-T_size(3)-1);end