//http://www.nauticom.net/www/jdtaft/JavaFFT.htm
//http://sepwww.stanford.edu/oldsep/hale/FftLab.html fft demo applet

import java.awt.*;
public class fft{

public static int DB=30; // kiek decibelu rodyti (keicia autocorr2)

public static double meanuofautocorr2;

public static complexArr fft1d(complexArr u){
 int n = u.length;
 double[][] U = new double[n][2];
 for (int i = 0; i < n; i++ ) { U[i][0]=u.re[i]; U[i][1]=u.im[i]; }
 return new complexArr(fft_1d(U));
}
public static complexArr ifft1d(complexArr uf){
 int n = uf.length;
 double[][] U = new double[n][2];
 for (int i = 0; i < n; i++ ) { U[i][0]=uf.re[i]; U[i][1]=-uf.im[i]; }
 complexArr u = new complexArr(fft_1d(U));
 for (int i = 0; i < n; i++ ) { u.im[i]=-u.im[i]; }
 return u;
}

public static complexArr2 fft2d(complexArr2 u){
 int h=u.h, w=u.w;
// my.print(h,"h"); my.print(w,"w");
 complexArr2 uf = new complexArr2(u.re,u.im);
 for (int i=0; i<h; i++) uf.putRow(i, fft1d(uf.getRow(i)));
 for (int j=0; j<w; j++) uf.putCol(j, fft1d(uf.getCol(j)));
 return uf;
}	
public static complexArr2 ifft2d(complexArr2 uf){
 int h=uf.h, w=uf.w;
// my.print(h,"h"); my.print(w,"w");
 complexArr2 u = new complexArr2(uf.re,uf.im);
 for (int i=0; i<h; i++) u.putRow(i, ifft1d(u.getRow(i)));
 for (int j=0; j<w; j++) u.putCol(j, ifft1d(u.getCol(j)));
 return u;
}	
public static int[] mexHatflt(byte[] u, int w, int h){
int alfa=(255 << 24), R=(255<<16), G=(255<<8), B=(255), Y=R|G;
int [][] flt=new int [][]{
	{ 0,  0, -1, -1, -1,  0,  0},
	{ 0, -1, -3, -3, -3, -1,  0},
	{-1, -3,  0, +7,  0, -3, -1},
	{-1, -3, +7, 24, +7, -3, -1},
	{-1, -3,  0, +7,  0, -3, -1},
	{ 0, -1, -3, -3, -3, -1,  0},
	{ 0,  0, -1, -1, -1,  0,  0}};
 int wh=w*h,w1=wh-1,w2=2*w,w3=3*w;
 int[] U = new int[wh], v = new int[wh], Sw = new int[wh], Sh = new int[wh];;
 for (int i=0; i<wh; i++) U[i]=127-(int)u[i];
 int[] b = new int[wh];
 int n=flt.length; 
// for (int i=0; i<n; i++) for (int j=0; j<n; j++) q+=Math.abs(flt[i][j]);
 for (int i=0; i<wh; i++){
	Sw[i]=U[Math.max(0,i-1)]+U[i]+U[Math.min(w1,i+1)];
	Sh[i]=U[Math.max(0,i-w)]+U[i]+U[Math.min(w1,i+w)];
 }
 int q=0, Q=0;
 for (int i=0; i<wh; i++){ v[i] = 24*U[i]+
	7*(U[Math.max(0,i-1)]+U[Math.min(w1,i+1)]+
		U[Math.max(0,i-w)]+U[Math.min(w1,i+w)])-
	3*(Sh[Math.max(0,i-2)]+Sh[Math.min(w1,i+2)]+
		Sw[Math.max(0,i-w2)]+Sw[Math.min(w1,i+w2)])-
	(Sh[Math.max(0,i-3)]+Sh[Math.min(w1,i+3)]+
		Sw[Math.max(0,i-w3)]+U[i]+Sw[Math.min(w1,i+w3)])-
	(U[Math.max(0,i-2-w2)]+U[Math.min(w1,i-2+w2)]+
		U[Math.max(0,i+2-w2)]+U[Math.min(w1,i+2+w2)]);
	if (Q<v[i]) Q=v[i]; else if (q>v[i]) q=v[i];
 }
 if (Q<-q) Q=-q;
 for (int i=0; i<wh; i++){ q=(255*v[i])/Q;
  if (q>0)
	b[i] = alfa | R | (255-q)<<8| (255-q);
  else b[i] = alfa | (255+q)<<16|G|(255+q);
 }
 return b;
}

public static byte[] autocorr2(byte[] u, int w, int h, int L, int tozero){
// u puse LxL viduriuko autokoreliacija su savimi
 if (tozero==0) DB*=3;
 int trsh=255/4; // kontrasto minimalus dydis kai dar apsimoka slaiciuoti
 int wh=w*h, L2=L/2, L4=L2/2;
 byte[] b = new byte[wh];
 double[] C = new double[L*L];
 float[] c = new float[wh];
 double[][] Ure = new double[L][L];
 double[][] ure = new double[L][L];
 boolean[][] B = new boolean[L][L];
 int iB = 0, R2=L2*L2;
 for (int i=0; i<L; i++) for (int j=0; j<L; j++)
	if (((i-L2)*(i-L2)+(j-L2)*(j-L2))<=R2) { iB+=1; B[i][j]=true;} else B[i][j]=false;
 double s,S;
 int p,umx,umi;
 complexArr2 uc; 
 int I,J=(h%L)/2,ii,jj,ij,II=0,JJ=0;
 I=(w%L)/2; while (I+L<=w) {II+=1; I+=L;}
 J=(h%L)/2; while (J+L<=h) {JJ+=1; J+=L;}
 double[][] cmax = new double[II][JJ], cmin = new double[II][JJ];
 double Cmax=0;

 J=(h%L)/2; S = 0; jj=0;
 while (J+L<=h){
  I=(w%L)/2; ii=0;
  while (I+L<=w) {
	umx=tozero-(int)u[I+L2+(L2+J)*w]; umi=umx; s=0;
	for (int i=0; i<L; i++) for (int j=0; j<L; j++){	
		Ure[i][j]=tozero-(int)u[I+i+(j+J)*w];
		if (B[i][j]){
			p = (int)Ure[i][j];
			s+=p;
			if (umx<p) umx=p; else if (umi>p) umi=p;
			ure[i][j]=p;
		}
	}
//	my.print(new int[]{umx,umi,umx-umi,trsh},"umx,umi,umx-umi,trsh");
	if (umx-umi>trsh){
		s/=iB;
		if (S<s) S=s;
//		my.print(new int[]{(int)umx,(int)umi,(int)((umx+umi)/2),(int)s},"umx,umi,(umx+umi)/2,s");
		uc = ifft2d(complexMult2(fft2d(new complexArr2(Ure)),
						conj(fft2d(new complexArr2(ure)))));
//		for (int i=0; i<L; i++) for (int j=0; j<L; j++) C[I+(i+L2)%L+((j+L2)%L+J)*w] = uc.re[i][j];
		for (int i=0; i<L; i++) for (int j=0; j<L; j++) {
			C[i+j*L] = uc.re[i][j]; c[I+(i+L2)%L+((j+L2)%L+J)*w]=(float)C[i+j*L];
		}
		cmax[ii][jj] = my.max(C)+1; Cmax=Math.max(Cmax,cmax[ii][jj]);
		cmin[ii][jj] = my.min(C);
	}
/*
	Cmin=Cmax*Math.exp(-DB/10/my.log10);
	double q = Math.log(Cmax/Cmin);
	for (int i=0; i<L; i++) for (int j=0; j<L; j++) b[I+(i+L2)%L+((j+L2)%L+J)*w] = (byte) (Math.log(Math.max(1,C[i+j*L]/Cmin))/q*255-128);
*/

	I+=L; ii+=1;
  }
  J+=L; jj+=1;
 }
 for (ii=0; ii<II; ii++) for (jj=0; jj<JJ; jj++) if (cmax[ii][jj]>0){	
	double cmx=(8*cmax[ii][jj]+Cmax)/9;
	double cmi=cmx*Math.exp(-DB/10/my.log10);
	double q = Math.log(cmx/cmi);
	I=(w%L)/2+ii*L; J=(h%L)/2+jj*L;
	if (cmi+1<cmin[ii][jj]){
		double pp=(cmax[ii][jj]-cmi)/(cmax[ii][jj]-cmin[ii][jj]);
		for (int i=0; i<L; i++) for (int j=0; j<L; j++){
			ij = I+i+(j+J)*w; c[ij] = (float)(pp*(c[ij]-cmin[ii][jj])+cmi);
		}
	}
	for (int i=0; i<L; i++) for (int j=0; j<L; j++){
		ij = I+i+(j+J)*w;
		b[ij] = (byte) (Math.log(Math.max(1,c[ij]/cmi))/q*255-128);
	}
 }
 meanuofautocorr2 = S;
 my.print((int) S,"meanuofautocorr2");
/*
 double Cmax = my.max(C), Cmin=Cmax*Math.exp(-DB/10/my.log10);
 double q = Math.log(Cmax/Cmin);
 for (int i=0; i<wh; i++) b[i] = (byte) (Math.log(Math.max(1,C[i]/Cmin))/q*255-128);
*/
 if (tozero==0) DB/=3;
 return b;
}

public static byte[] autocorr2(int[] u, int w, int h, int L){
 return autocorr2(int2byte(u),w,h,L,0);
}


public static byte[] int2byte(int[] u){
 int n=u.length, r, g, G=(255<<8),B=255,j,q;
 byte[] b = new byte[n];
 for (int i=0; i<n; i++){
	j = u[i]; if ((j&G)==G)	b[i]=(byte)(((j&B)-B)/2);
		else b[i]=(byte)((B-(j&B))/2);
	}
// if (q>0){ b[i] = alfa | R | (255-q)<<8| (255-q);
//	  else b[i] = alfa | (255+q)<<16|G|(255+q);
 return b;
}

public static double[][] fft_1d( double[][] array )
{
	double  u_r,u_i, w_r,w_i, t_r,t_i;
	int     ln, nv2, k, l, le, le1, j, ip, i, n;

	n = array.length;
    	ln = (int)( Math.log( (double)n )/Math.log(2) + 0.5 );
    	nv2 = n / 2;
    	j = 1;
 	for (i = 1; i < n; i++ )
    	{
		if (i < j)
		{
	    		t_r = array[i - 1][0];
	    		t_i = array[i - 1][1];
	    		array[i - 1][0] = array[j - 1][0];
	    		array[i - 1][1] = array[j - 1][1];
	    		array[j - 1][0] = t_r;
	    		array[j - 1][1] = t_i;
		}
		k = nv2;
		while (k < j)
		{
	    		j = j - k;
	    		k = k / 2;
		}
		j = j + k;
    	}

 	for (l = 1; l <= ln; l++) /* loops thru stages */
    	{	
     	 	le = (int)(Math.exp( (double)l * Math.log(2) ) + 0.5 );
	  	le1 = le / 2;
		u_r = 1.0;
		u_i = 0.0;
		w_r =  Math.cos( Math.PI / (double)le1 );
		w_i = -Math.sin( Math.PI / (double)le1 );
		for (j = 1; j <= le1; j++) /* loops thru 1/2 twiddle values per stage */
		{
	    		for (i = j; i <= n; i += le) /* loops thru points per 1/2 twiddle */
	    		{
				ip = i + le1;
				t_r = array[ip - 1][0] * u_r - u_i * array[ip - 1][1];
				t_i = array[ip - 1][1] * u_r + u_i * array[ip - 1][0];

				array[ip - 1][0] = array[i - 1][0] - t_r;
				array[ip - 1][1] = array[i - 1][1] - t_i; 

				array[i - 1][0] =  array[i - 1][0] + t_r;
				array[i - 1][1] =  array[i - 1][1] + t_i;  
	    		} 
	    		t_r = u_r * w_r - w_i * u_i;
	    		u_i = w_r * u_i + w_i * u_r;
	    		u_r = t_r;
		} 
    	}
	double q = ((double) 1)/Math.sqrt((double) n);
	for (i = 0; i < n; i++){ array[i][0] *=q; array[i][1] *=q;}
	return array;
} /* end of FFT_1d */

 public static complexArr2 complexMult2(complexArr2 u, complexArr2 v){
	int h=u.h, w=u.w;
	complexArr2 uv = new complexArr2(h,w);
	for (int i=0; i<h; i++) for (int j=0; j<w; j++){
		uv.re[i][j]=u.re[i][j]*v.re[i][j]-u.im[i][j]*v.im[i][j];
		uv.im[i][j]=u.re[i][j]*v.im[i][j]+u.im[i][j]*v.re[i][j];
	 	}
	return uv;
 }
 public static complexArr2 conj(complexArr2 u){
	int h=u.h, w=u.w;
	complexArr2 v = new complexArr2(h,w);
	for (int i=0; i<h; i++) for (int j=0; j<w; j++){
		v.re[i][j]=u.re[i][j]; v.im[i][j]=-u.im[i][j];
	 	}
	return v;
 }
}
class complexArr{
 public double[] re;
 public double[] im;
 int length;
 public complexArr(double[] Re, double[] Im){
  length=Re.length;
  re= new double[length];
  im= new double[length];
  for (int i=0; i<length; i++) {re[i]=Re[i]; im[i]=Im[i];}
 }
 public complexArr(double[][] u){
  length=u.length;
  re= new double[length];
  im= new double[length];
  for (int i=0; i<length; i++) {re[i]=u[i][0]; im[i]=u[i][1];}
 }
 public double getE(){
	double e=0; for (int i=0; i<length; i++) e+=re[i]*re[i]+im[i]*im[i];
	return e;
 }
} //complexArr

class complexArr2{
 public double[][] re;
 public double[][] im;
 int[] length;
 int h, w;
 public complexArr2(int h, int w){
  this.h = h; this.w = w; length = new int[]{h, w};
  re= new double[h][w]; im= new double[h][w];
 }
 public complexArr2(double[][] Re, double[][] Im){
  h=Re.length; w=Re[0].length; length = new int[]{h, w};
  re= new double[h][w]; im= new double[h][w];
  for (int i=0; i<h; i++) for (int j=0; j<w; j++) 
	{re[i][j]=Re[i][j]; im[i][j]=Im[i][j];}
 }
 public complexArr2(double[][] Re){
  h=Re.length; w=Re[0].length; length = new int[]{h, w};
  re= new double[h][w]; im= new double[h][w];
  for (int i=0; i<h; i++) for (int j=0; j<w; j++) re[i][j]=Re[i][j];
 }

 public complexArr getCol(int J){
	double[] ur = new double[h], ui = new double[h];
	for (int i=0; i<h; i++) { ur[i] = re[i][J]; ui[i] = im[i][J]; }
	return new complexArr(ur,ui);
 }
 public complexArr getRow(int I){
	double[] ur = new double[w], ui = new double[w];
	for (int j=0; j<w; j++) { ur[j] = re[I][j]; ui[j] = im[I][j]; }
	return new complexArr(ur,ui);
 }
 public void putCol(int J, complexArr u){
	for (int i=0; i<h; i++) { re[i][J] = u.re[i]; im[i][J] = u.im[i]; }
 }
 public void putRow(int I, complexArr u){
	for (int j=0; j<w; j++) { re[I][j] = u.re[j]; im[I][j] = u.im[j]; }
 }
 public double getE(){
	double e=0;
	for (int i=0; i<h; i++) for (int j=0; j<w; j++)
		e+=re[i][j]*re[i][j]+im[i][j]*im[i][j];
	return e;
 }
} //complexArr2