function [p_Igivenv,Vx,Vy] = ComputeVelLikeli(Ix,Iy,It,Measure_sigma,pixellocation_ijt,window_sigma,vellims)
% Ix,Iy,It:  Image derivatives for a movie in format returned by ComputeImageDerivs
% Measure_sigma: Gaussian noise standard deviation assumed for image
%               measurements-- controls width of likelihood in speed
% pixellocation_ijt: pixel and frame around which velocity computation is
%                   performed
% window_sigma:  controls the width of the averaging Gaussian window (e.g.
%               window_sigma= 1 standard deviation produces a weighted average 
%               over a 2-3 pixel region
% vellims= [velmin, velmax]:  optional argument allows you to specify the minimum and maximum velocities the
%               likelihood will be computed over.
%
% The function computes the velocity likelihood centered at pixellocation i,j, and frame t (e.g.
% middle of a 64x64 image, and middle frame=16 ~ pixellocation_ijt = [32 32 16].  
% 
% The function automatically computes the velocity range and returns
% the conditional probability at each velocity sampled over a 64x64 grid
% p_Igivenv=p(vx,vy|Ix,Iy,It)  (however, the normalization constant is
%       neglected)
% 

ix = pixellocation_ijt(1);
jy = pixellocation_ijt(2);
t = pixellocation_ijt(3);
% compute mean velocity
[px,py] = meshgrid(1:size(Ix,1),1:size(Ix,2));
w = exp(-0.5*((px-ix).^2+(py-jy).^2)/window_sigma^2);
sIxx = sum(sum(Ix(:,:,t).^2.*w));
sIxy = sum(sum(Ix(:,:,t).*Iy(:,:,t).*w));
sIyy = sum(sum(Iy(:,:,t).^2.*w));
sIxt = sum(sum(Ix(:,:,t).*It(:,:,t).*w));
sIyt = sum(sum(Iy(:,:,t).*It(:,:,t).*w));
sItt = sum(sum(It(:,:,t).*It(:,:,t).*w));


if nargin<7,
    meanv = inv([sIxx sIxy; sIxy sIyy])*[sIxt;sIyt];
    vlimmin = -2*norm(meanv);
    vlimmax = 2*norm(meanv);
else,
    vlimmin =vellims(1);
    vlimmax = vellims(2);
end
    
vrange = vlimmax-vlimmin;
[Vx,Vy] = meshgrid(vlimmin:vrange/63:vlimmax,vlimmin:vrange/63:vlimmax);

% the likelihood is given by 
% exp(1/(2sig^2) {Sum_xy w(x,y)(Ix vx + Iy vy +It)^2})
%  {Sum_xy w(x,y)(Ix vx + Iy vy +It)^2} =
%   Sum_xy w(x,y)[(Ix vx + Iy vy +It)(Ix vx + Iy vy +It)]
%   Sum_xy w(x,y)[Ix^2 vx^2 +2*IxIy vxvy + 2*IxIt vx +2*IyIt vy +It^2]
%   (Sum_xy w(x,y)Ix^2) vx^2 +  (Sum_xy w(x,y)Iy^2) vy^2 + ...
%   sIxx vx^2 + sIyy vy^2 +...

p_Igivenv = exp(-0.5*(sIxx*Vx.^2+sIyy*Vy.^2+2*sIxy*Vx.*Vy+2*sIxt*Vx+2*sIyt*Vy+sItt)/Measure_sigma^2);