www.pudn.com > snippets(1).rar > flowmap_withties.m, change:2009-10-12,size:2952b


function [Fy,Fx] = flowmap(I,winsz,maxv,ssd_flag,smooth_flag)
% function [Fy,Fx] = flowmap(I,winsz,maxv,ssd_flag)
%
% get optic flow in x- and y-direction, by direct template matching
%
% I ... input sequence
% winsz ... size of matching window
% maxv ... maximum possible flow between neighbouring frames
% ssd_flag ... use sum of ABSOLUTE or SQUARED differences ?
% smooth_flag ... smooth integer flow-field after computations ?

% show estimated flow-field ?
DISPLAY = 0;

% preparations
[ht,wt,fr] = size(I);

I = single(I);

Fy = zeros(ht,wt,fr-1,'single');
Fx = zeros(ht,wt,fr-1,'single');
hy = zeros(ht,wt,'single');
hx = zeros(ht,wt,'single');

% window size should be odd
if mod(winsz,2)==0
  winsz = winsz+1;
end

woff = floor(winsz/2);
numpos = (2*maxv+1)^2;

% filter definitions
fsum = fspecial('gaussian',winsz,winsz/2);
fvsum = fsum(ceil(winsz/2),:); % convolution is separable -> much faster
fsmooth = fspecial('gaussian',winsz,winsz/2);

% loop through frames
for u = 2:fr

  prev = I(:,:,u-1);
  curr = I(:,:,u);

  flowerr = zeros(ht,wt,numpos);
  offset = zeros(numpos,2);
  fct = 1;

  % shift previous frames to all possible positions within max-velocity
  for v = -maxv:maxv
    for w = -maxv:maxv

      % compute difference image
      if ssd_flag % ...with sum of squared differences
	adiff = (curr-imgshift(prev,v,w)).^2;
      else % ...or with sum of absolute differences
	adiff = abs(curr-imgshift(prev,v,w));
      end
      % sum locally in each receptive field
      % convolution is separable -> much faster
      % flowerr(:,:,fct) = conv2(adiff,fsum,'same');
      flowerr(:,:,fct) = conv2(fvsum,fvsum',adiff,'same');
      offset(fct,:) = [w v];
      fct = fct+1;

    end % vor w
  end % vor v

  % find flow vector with smallest local difference
  [smallest,best] = min(flowerr,[],3);
  % resolve ties: prefer small flow
  bestwithties = zeros(ht,wt);
  for v = 1:ht
    for w = 1:wt
      list = find(flowerr(v,w,:)==smallest(v,w)); % find all with minimum value
      dist = sum(offset(list).^2,2);              % get their flow magnitude
      [dummy,pos] = min(dist);                    % get the one with smallest flow
      best(v,w) = list(pos);                      % replace best value
    end
  end

  for v = 1:numpos
    idx = (best==v);
    hy(idx) = offset(v,1);
    hx(idx) = offset(v,2);
  end
  % smooth flow-field in receptive field (lateral interaction!)
  if smooth_flag
    hy = conv2(hy,fsmooth,'same');
    hx = conv2(hx,fsmooth,'same');
  end

  Fy(:,:,u-1) = hy;
  Fx(:,:,u-1) = hx;

  % show flow field overlayed on image
  if DISPLAY

    rst = maxv+3;
    figure;
    imshow(I(:,:,u)/255);
    hold on;
    for a = maxv:rst:ht
      for b = maxv:rst:wt
        if Fx(a,b,u-1)<0.01 && Fy(a,b,u-1)<0.01, continue; end
	plot([b,b+Fx(a,b,u-1)],[a,a+Fy(a,b,u-1)],'r-','linewidth',2);
	plot(b,a,'go','markersize',3,'markerfacecolor','g');
      end
    end

  end % if DISPLAY

end % for u