correspwarp.m

% Written by Carl DOersch (cdoersch at cs dot cmu dot edu)
%
% A utility to convert a corresp datastructure into an image warping
%
% corresp: as generated by contextpredict
%
% sbins: HOG sbins parameter.  Output is of size (size(corresp)+2)*sbins.
%
% im: the predictor image to be warped into the coordinates of whatever query image
%     generated corresp.
%
% pyrs: the full set of pyramids used to generate corresp.
%
% pyridx: the index into corresp{y,x}.mu, corresp{y,x}.covar, and pyrs that corresponds
%         to im.
% 
% conf: struct with additional optional parameters.  They may include:
%
%    - nowarp: don't actually warp the image; if 1, imwarp will be returned as [].
%
% imwarp: the warped image
%
% warpx and warpy: matrices the same size as im(:,:,1) saying where in the query image
%                  each pixel in im corresponds to.
%
% sumtotal: the total amount of probability mass that fell on each pixel in warpx and warpy.
%
function [imwarp,warpx,warpy,sumtotal]=correspwarp(corresp,sbins,im,pyrs,pyridx,conf)
if(~exist('conf','var'))
  conf=struct();
end
try
  im=im2double(im);

  % First goal is to compute warpx and warpy.  This is done by computing the coordinates
  % of the center of each cell of corresp in the query image.  Then for each pixel in
  % warpx and warpy, we average all the locations in the query image which correspond
  % to it.  To build up this average for a given location u,v in warpx and warpy,, we 
  % iterate over all cells in corresp and add up
  % the locations where corresp{y,x} has high probability of corresponding to u,v according
  % to the gaussian defined by corresp{y,x}.  
  %
  % Note that this can result in a fisheye effect if the covariance matrices are large,
  % since points in warpx and warpy may fall under some of the gaussian bumps
  % even they are far from the mean of any of these gaussians.
  warpx=zeros(size(im(:,:,1)));
  warpy=warpx;
  sumtotal=warpx;
  for(x=1:size(corresp,2))
    for(y=1:size(corresp,1))
      curpos=[(x+.5)*sbins (y+.5)*sbins]+.5;
      mycorresp=corresp{y,x};
      if(isempty(mycorresp))
        continue;
      end
      currcovar=reshape(mycorresp.covar(:,pyridx),[2,2]);
      % the window in warpx and warpy that's affected by this particular gaussian is
      % computed in the same way as in contextpredict, but there's an extra
      % coordinate transformation because corresp is in HOG cell coordinates, which isn't
      % the same as pixels.
      [V,D]=eig(currcovar);
      dist=sqrt(abs((log(.0001)+log(det(currcovar))/2)/(V(:,1)'*inv(currcovar)*V(:,1)/2)));
      dist2=sqrt(abs((log(.0001)+log(det(currcovar))/2)/(V(:,2)'*inv(currcovar)*V(:,2)/2)));
      V=max(abs(dist*V(:,1)),abs(dist2*V(:,2)));
      pyrlevel=mycorresp.level(pyridx);
      coordscale=sbins*pyrs{pyridx}.scales(pyrlevel);
      coordoffset=(.5)*coordscale+.5;
      windowx=round((((corresp{y,x}.mu(1,pyridx)-V(1))*coordscale):((corresp{y,x}.mu(1,pyridx)+V(1))*coordscale))+coordoffset);
      windowy=round((((corresp{y,x}.mu(2,pyridx)-V(2))*coordscale):((corresp{y,x}.mu(2,pyridx)+V(2))*coordscale))+coordoffset);
      windowxinhog=(windowx-coordoffset)/coordscale;
      windowyinhog=(windowy-coordoffset)/coordscale;
      [gridx,gridy]=meshgrid(windowxinhog,windowyinhog);
      pts=[gridx(:)';gridy(:)'];
      pts=bsxfun(@minus,pts,mycorresp.mu(:,pyridx));
      probs=reshape(1/(2*pi*sqrt(det(currcovar)))*exp(-sum(pts.*(inv(currcovar)*pts),1)/2),size(gridx));
      probs=probs./sum(probs(:));%*confidence(y,x);
      probs(windowy<1 | windowy>size(im,1),:)=[];
      probs(:,windowx<1 | windowx>size(im,2))=[];
      windowx(windowx<1 | windowx>size(im,2))=[];
      windowy(windowy<1 | windowy>size(im,1))=[];
      warpx(windowy,windowx)=warpx(windowy,windowx)+curpos(1)*probs;
      warpy(windowy,windowx)=warpy(windowy,windowx)+curpos(2)*probs;
      sumtotal(windowy,windowx)=sumtotal(windowy,windowx)+probs;
    end
  end
  warpx=warpx./sumtotal;
  warpy=warpy./sumtotal;
  if(dsbool(conf,'nowarp'))
    imwarp=[];
    return
  end
  % now that we have warpx and warpy computed, TriScatteredInterp will do all the work
  % of computing the warped image.  We do, however, need to get rid of points in
  % warpx and warpy that don't correspond to anything.
  [outgridx,outgridy]=meshgrid(1:(size(corresp,2)+2)*sbins,1:(size(corresp,1)+2)*sbins);
  for(i=1:size(im,3))
    imch=im(:,:,i);
    valid=~isnan(warpx(:))&~isinf(warpx(:))&~isnan(warpy(:))&~isinf(warpy(:));
    F=TriScatteredInterp(warpx(valid),warpy(valid),imch(valid));
    imwarp(:,:,i)=F(outgridx,outgridy);
  end
  imwarp=max(min(imwarp,1),0);
catch ex,dsprinterr;end
end