title: Horn-Schunck Optical Flow (Averaging Velocity)
body: <p>I'm currently trying to implement the HS-method for optical flow but my <code>u</code> and <code>v</code> always seem to have only zeros in them. I can't seem to figure out my error in here:</p>
  
  <pre><code>vid=VideoReader('outback.AVI');
  vid.CurrentTime = 1.5;
  
  alpha=1;
  iterations=10;    
  
  frame_one = readFrame(vid);
  vid.CurrentTime = 1.6;
  frame_two = readFrame(vid);
  
  % convert to grayscale
  
  fone_gr = rgb2gray(frame_one);
  ftwo_gr = rgb2gray(frame_two);
  
  % construct  for each image
  sobelx=[-1 -2 -1; 0 0 0; 1 2 1];
  sobely=sobelx';
  time=[-1 1];
  
  fx_fone=imfilter(fone_gr, sobelx);
  fy_fone=imfilter(fone_gr, sobely);
  ft_fone=imfilter(fone_gr, time);
  
  fx_ftwo=imfilter(ftwo_gr, sobelx);
  fy_ftwo=imfilter(ftwo_gr, sobely);
  ft_ftwo=imfilter(ftwo_gr, time);
  
  Ix=double(fx_fone+fx_ftwo);
  Iy=double(fy_fone+fy_ftwo);
  It=double(ft_fone+ft_ftwo);
  
  % flow-variables (velocity = 0 assumption)
  velocity_kernel=[0 1 0; 1 0 1; 0 1 0];
  
  u = double(0);
  v = double(0);
  
  % iteratively solve for u and v
  for i=1:iterations
      neighborhood_average_u=conv2(u, velocity_kernel, 'same');
      neighborhood_average_v=conv2(v, velocity_kernel, 'same');
  
      data_term = (Ix .* neighborhood_average_u + Iy .* neighborhood_average_v + It);
      smoothness_term = alpha^2 + (Ix).^2 + (Iy).^2;
  
      numerator_u = Ix .* data_term;
      numerator_v = Iy .* data_term;
  
      u = neighborhood_average_u - ( numerator_u ./ smoothness_term );
      v = neighborhood_average_v - ( numerator_v ./ smoothness_term );
  end
  
  u(isnan(u))=0;
  v(isnan(v))=0;
  
  figure
  imshow(frame_one); hold on;
  quiver(u, v, 5, 'color', 'b', 'linewidth', 2);
  set(gca, 'YDir', 'reverse');
  </code></pre>
  
  <p>The only thing I'm not really confident about is the computation of the neighborhood average:</p>
  
  <pre><code>velocity_kernel=[0 1 0; 1 0 1; 0 1 0];
  
  u = double(0);
  v = double(0);
  
  [..]
  neighborhood_average_u=conv2(u, velocity_kernel, 'same');
  neighborhood_average_v=conv2(v, velocity_kernel, 'same');
  </code></pre>
  
  <p>Wouldn't that always result in a convolution matrix with only zeros?</p>
  
  <p>I thought about changing it to the following, since I need to compute the average velocity using the velocity kernel on each pixel of my images:</p>
  
  <pre><code>velocity_kernel=[0 1 0; 1 0 1; 0 1 0];
  u = double(0);
  v = double(0);
  
  [..]
  neighborhood_average_u=conv2(Ix, velocity_kernel, 'same');
  neighborhood_average_v=conv2(Iy, velocity_kernel, 'same');
  </code></pre>
  
  <p>But I still don't know if that would be the correct way. I followed the instructions on the bottom of this MATLAB page: 
  <a href="http://de.mathworks.com/help/vision/ref/opticalflowhs-class.html" rel="nofollow">http://de.mathworks.com/help/vision/ref/opticalflowhs-class.html</a></p>