ImageRecognition/ismdim_matching.m

function [params]=ismdim_matching(I,patches,patchClassList,Locations,sigmaLGN,scales,reScale,flip)
% Performs the matching process of the implicit shape model (ISM). The method
% differs from the standard in (1) the way patches are matched to the image
% (this method uses DIM), and (2) the way voting is performed (this method uses
% DIM).
%
% Inputs:
%    I - an image on which object detection is to be performed
%    patches - the codebook image patches
%    patchClassList - the object class associated with each codebook patch
%    Locations - the locations of the votes cast by each codebook patch
%    sigmaLGN - sets of the size of the Gaussian used to preprocess grayscale images
%
% patches, patchClassList, Locations, and sigmaLGN are all output by the ismdim_training function.
%
% Output:
%    params - the parameters of potential object locations 
%
% params is a cell array. The number of cells equals the number of classes. Each
% cell array is a 3xn or 4xn matrix, where n is the number of potential object
% locations found. For each location the first parameter is the amplitude of the
% votes at that location, the following 2 or 3 parameters are the coordinates of
% that location (in 2D or 3D).
%
% NOTE: this implementation confounds scale and class!!!

[numTemplates,patchLen]=size(patches);
patchLen=sqrt(patchLen);
patchHalfLen=(patchLen-1)/2;
ONOFF=1; 
recon=0;

%CONVERT PATCH CODEBOOK TO DIM WEIGHTS
if nargin>6 && reScale==1
  %convert single scale codebook to a multiscale codebook
  [w,v]=define_dictionary(patches,ONOFF,sigmaLGN,scales);
  disp('expanded single-scale coodebook into a multiscale codebook (patches)')
else
  [w,v]=define_dictionary(patches,ONOFF,sigmaLGN);
end
if nargin>7 && flip==1
  [nMasks,nChannels]=size(w);
  for i=1:nMasks
    for j=1:nChannels
      w{i+nMasks,j}=fliplr(w{i,j});
      v{i+nMasks,j}=fliplr(v{i,j});
    end
  end
  disp('added flipped image patches to the codebook');
end

%MATCH CODEBOOK PATCHES TO IMAGE using DIM
if ONOFF
  [~,~,X{1},X{2},trueRange]=imnorm(I,sigmaLGN,[],1);
  y=dim_activation_conv_recurrent(w,X,v,[],[],trueRange);
else
  y=dim_activation_conv_recurrent(w,{I},v,[]);
end
%ymax=max(max(max(cat(3,y{:}),[],3)))

%select only those matches that come from voting elements (non-background
%patches) and exceed a threshold
indPatchesInClass=find(patchClassList>0); %all non-background patches
threshold=1e-3;
k=0;
for j=indPatchesInClass
  k=k+1;
  y{j}=imregionalmax(y{j}).*y{j};
  y{j}(y{j}<threshold)=0;
  Mmaps{k}=y{j};
  if recon
    if ONOFF
      Vr{k}=v{j,1}-v{j,2};
    else
      Vr{k}=v{j,1};
    end
  end
end

%plot matches
figured(1),clf, plot_image(I), hold on, colormap('gray');
set(gca,'FontSize',24);
set(gcf,'PaperSize',[20 15],'PaperPosition',[0 0 20 15]);
MmapAll=max(cat(3,Mmaps{:}),[],3);
step=0.1;
for thres=threshold:step:1
  [ptx,pty]=find(MmapAll>thres);
  c=[0.1,min(1,0.5+thres),0.1];
  plot(pty,ptx,'o','LineWidth',1,'Color',c,'MarkerSize',ceil(2*thres/step),'MarkerFaceColor',c); 
end


%APPLY GENERALISED HOUGH TRANSFORM TO COUNT VOTES using DIM
tic
%calculate the accumulator array
%if prod(size(I))*prod(size(w))<1e8
%method that used the matrix multiplication version of DIM - fast, memory intensive:
%H=ismdim_voting(Locations,Mmaps);
%else
%alternative method using the convolutional version of DIM - slow, but uses much less memory:
Hc=dim_activation_conv_recurrent(Locations,Mmaps); H=cat(3,Hc{:});
%end
toc

%find peaks in accumulator array
numClasses=size(Locations,1);
if nargin>5 && ~isempty(scales)
  %find peaks in a 3-dimensional accumulator array where dimentions represent x-
  %and y-coordinates and scale
  [params{1},Hone_peak{1}]=find_accumulator_peaks(H,[],[],scales);
  for c=2:numClasses, params{c}=[]; Hone_peak{c}=[]; end, 
else
  for c=1:numClasses
    %for each class separately find peaks in a 2-dimensional accumulator array where
    %dimentions represent x- and y-coordinates
    Htmp=H(:,:,c);
    [params{c},Hone_peak{c}]=find_accumulator_peaks(Htmp,[],[],[]);
  end
end

%plot accumulator array
figured(2),clf,
[plotRows,plotCols]=best_subplot_rows_cols(numClasses);
for c=1:numClasses
  maxsubplot(plotRows,plotCols,c);
  plot_accumulator(H(:,:,c),0.95*max(H(:)),params{c});
end

if recon
  %plot back-projections of best matches
  if nargin>5 && ~isempty(scales), 
    for match=1:3
      figured(3+match),clf,
      if length(Hone_peak{1})>=match 
        for c=1:numClasses
          Hc{c}=Hone_peak{1}{match}(:,:,c);
        end
        [~,~,R1]=dim_activation_conv_recurrent(Locations,Mmaps,[],Hc,1);
        %R1All=max(cat(3,R1{:}),[],3); subplot(1,2,1), plot_image(R1All);
        R0=zeros(size(I));
        for j=1:length(indPatchesInClass)
          R0=R0+conv2(R1{j}.*Mmaps{j},Vr{j},'same');
        end
        plot_image(R0);
      end
    end
  else
    for match=1:3
      figured(3+match),clf,
      for c=1:numClasses
        if length(Hone_peak{c})>=match
          [~,~,R1]=dim_activation_conv_recurrent(Locations,Mmaps,[],{Hone_peak{c}{match}},1);
          %R1All=max(cat(3,R1{:}),[],3); subplot(1,2,1), plot_image(R1All);
          R0=zeros(size(I));
          for j=1:length(indPatchesInClass)
            R0=R0+conv2(R1{j}.*Mmaps{j},Vr{j},'same');
          end
          subplot(plotRows,plotCols,c);
          plot_image(R0);
        end
      end
    end
  end
  drawnow;
end


function plot_accumulator(H,scale,params);
if nargin<2 || isempty(scale)
  scale=max(1,0.95*max(max(H)));
end
plot_image(H,[0,max(1e-6,scale)]), hold on, 
set(gca,'FontSize',24);
set(gcf,'PaperSize',[20 15],'PaperPosition',[0 0 20 15]);
cmap=colormap('gray');cmap=1-cmap;colormap(cmap);

if nargin>=3 && ~isempty(params)
  hold on
  amplitude=params(1,:);
  coords=params(2:3,:);
  for i=1:min(3,length(amplitude))
    text(coords(1,i),coords(2,i),num2str(amplitude(i),'% 1.2f'),'Color',[0.1,0.8,0.1],'VerticalAlignment','bottom');
  end
end
drawnow;



function [param,Hone_peak_ordered]=find_accumulator_peaks(H,T,R,S)
%binarize accumulator array
Hbin=H;
Hbin(Hbin<1e-3)=0;
Hbin(Hbin>=1e-3)=1; 
%divide accumularor array into separate connected components
Hconnected=bwconncomp(Hbin); 

%find parameters associated with each component
L=labelmatrix(Hconnected);
centre=zeros(ndims(H),0);
amplitude=zeros(1,0);
for k=1:Hconnected.NumObjects
  Hone_peak{k}=zeros(size(H));
  Hone_peak{k}(L==k)=1;
  Hone_peak{k}=Hone_peak{k}.*H;

  centre(:,k)=centroid(Hone_peak{k},T,R,S)';
  amplitude(k)=sum(Hone_peak{k}(:));
  %amplitude(k)=sum(Hone_peak{k}(:).^2).^0.5;
end
%return parameters in order of greatest to weakest amplitude (i.e. most to fewest votes)
[~,order]=sort(amplitude,'descend');
param=[amplitude(order);centre(:,order)];
Hone_peak_ordered=[];
for i=1:length(order)
  Hone_peak_ordered{i}=Hone_peak{order(i)};
end