clear;

idx = ceil(84*rand);
%idx = 27;
fn  = sprintf('/disk2/people/anuj/research/Operators/Pictures/P%03d',idx);
I = double(imread(fn,'jpeg'));
n = size(I);
if length(n) > 2
	J = I(:,:,1);
	clear I;
	I = J;
	clear J;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Processing
%Id = diff(I);
Id = makeFilt(20,8,I,rand*180);
X = Id(:);
X = X - mean(X);
ku = kurtosis(X);
if ku < 3.05
	ku = 3.05;
end
p = 3/(ku-3);
c = var(X)/p;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Histogram Calculation
[N,t] = hist(X,30);
x1 = min(t);
x2 = max(t);
del = t(2) - t(1);
Nsum = sum(N);
N = N/(Nsum*del);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the bessel stuff for these parameters
sp = sqrt(pi);
gp = gamma(p);

u = (2^(-p + 0.5))*((2/c)^(p/2 + 0.25))/(sp*gp);

num = size(t);
for i = 1:num(2)
   y = abs(t(i));
   f(i) = u*(y^(p - 0.5))*besselk(p-0.5,y*sqrt(2/c));
end
nr = sum(f*del);
f = f/nr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Plotting
subplot(2,2,1);
imagesc(I);
colormap(gray);
axis('image');
title(fn(48:length(fn)));

subplot(2,2,4);
plot(t,N,'g--',t,f,'r');
str = sprintf('p=%0.3f, c=%0.2f',p,c);
title(str);
%legend('Observed','Estimated');
grid;
xlabel('Intensity in Filtered Image');
ylabel('Relative Frequency');
