function [g,A] = estimateG(laserhdrimg, tau, zd, min_radius, max_radius)

% [g,A] = estimateG(laserhdrimg, tau, zd, zw, min_radius, max_radius)
%
% Estimate the g parameter of the phase function (Heyney-Greenstein). g
% should be in [-1,1]. And also estimate the scale factor. From the HDR
% image, hdrlaser, we will use the data in the ring region from
% min_radius to max_radius. These two are in [0,1]. 1 means half size of
% the width of hdrlaser. For a radius in this range, we compute the average
% along the circle and use it as a point for fitting the HG function.
%
% Input:
% 
% laserhdrimg   HDR image of the laser scattering pattern (see getHDRlaser.m)
% tau           optical thickness of the dirty pattern
% zd            distance from the window to the back board ( in pixel )
% min_radius    the minimum radius of the data used for fitting, in [0,1]
% max_radius    the maximum radius of the data used for fitting, in [0,1]
%
% Output:
%
% g           parameter for the HG function
% A           scale factor
%
%
% Jinwei Gu
% jwgu@cs.columbia.edu
% 2006/11/2

% parsing the input parameters
if nargin < 4
    min_radius = 0.3;
end

if nargin < 5
    max_radius = 0.7;
end

%--------------------- main ----------------------

a = laserhdrimg(:,:,2); % green laser
wsize = size(laserhdrimg, 1);
xc = wsize/2;
yc = wsize/2;

rmin = min_radius*wsize/2;
rmax = max_radius*wsize/2;

I = [];
for r = rmin : 2 : rmax
    
    n = round(2*pi*r);
    tt = linspace(0,2*pi,n);
    
    j = round(xc + r*cos(tt));
    i = round(yc + r*sin(tt));
    
    f = 0;
    for k=1:length(tt)
        ta = a(i(k)-1:i(k)+1,j(k)-1:j(k)+1);   
        f = f+mean(ta(:));
    end
    f = f/length(tt);
    theta = atan2(r,zd);
    
    mu = cos(theta);
    f = (1-mu)*f/(exp(-tau)-exp(-tau/mu))/(mu^4);

    I = [I; theta, f];
end

% fitting
xdata=I(:,1);
ydata=I(:,2);
[estimates, model]=fit_g_curve(xdata,ydata);
g = estimates(1);
A = estimates(2);

% plot the result
figure;
hold on;
plot(I(:,1),I(:,2),'ro');
x=sort(xdata);
y=A*(1-g*g)./((1+g*g-2*g*cos(x)).^(1.5));
plot(x,y);
title(sprintf('g=%f',g));
xlabel('\theta');
ylabel('Radiance');
hold off;