# Hermite polynomials for 128x128
# naive implementation 
if( ! exist( "I0") )
I=imread("lena.ppm") ;
I0 = double(I)/256-0.5;
endif
i512 = max(size(I0));

rmax = 128;
#radius=round(rmax/(2*pi));
radius=(rmax/(2*pi));
r128=sqrt(rmax/(2*pi));
[Xm,Ym] = size(I0);
X= ((1:Xm)'-1-round(Xm/2)) / radius;

Y= ((1:rmax)'-1-round(rmax/2)) / r128;

H0X = exp(-0.5*X.^2)/sqrt(radius*sqrt(pi));
#H0X /= norm(H0X);
H1X = 2*X.*H0X;
HX=[H0X,H1X];
H0Y = exp(-0.5*Y.^2); H0Y /= norm(H0Y);
HY = [H0Y,2*Y.*H0Y];
for (k=2:rmax)
  HX(:,k+1) =(X.*HX(:,k) - (k-1)* HX(:,k-1))*2;
  HY(:,k+1) =(Y.*HY(:,k) - (k-1)* HY(:,k-1))*2; endfor;
c = 1;
for (k=1:rmax) c *= sqrt(k*2); HX(:,k+1) /= c; HY(:,k+1) /= c; endfor


H=[];
for( j=1:rmax)
for( k=1:Xm) Hz(k,j) = sum( I0(k,:) .* HX(:,j)' ) ; endfor
for( i=1:(1+rmax-j))
for( k=1:Xm) H(i,j) = sum( Hz(:,j).* HX(:,i) ) ; endfor
endfor
endfor

I1x= zeros(4*Xm,Xm);
I1f = zeros(rmax,rmax);
for( j=1:rmax)
for( i=1:(1+rmax-j))
I1x(Xm*rem(j+i,4)+(1:Xm),:) += H(i,j) * (HX(:,i)*HX(:,j)');
I1f += 1i ^ (j+i-2) * H(i,j) * (HY(:,i)*HY(:,j)');
endfor
endfor

nshow((I1x(Xm*0+(1:Xm),:)+I1x(Xm*1+(1:Xm),:)+I1x(Xm*2+(1:Xm),:)+I1x(Xm*3+(1:Xm),:)))

#nshow(H(2:2:rmax-1,1:2:rmax-1)/200)

#plot(HX(:,1:5))
X2 = (X*ones(1,Xm)).^2 + (ones(Xm,1)*X').^2;
Y2 = (Y*ones(1,rmax)).^2 + (ones(rmax,1)*Y').^2;
#D = I1f .* (Y2<200);
#for i=1:5000; D -= 0.0001i*pi*Ham(D,Y2/2,sqrt(2)*r128).*(Y2<195); endfor
