import java.awt.*;
import java.util.*;
import java.awt.image.*;
public class ComplexImage //implements ImageConsumer
{
    public int w;
    public int h;
    
    public boolean complete;
    
    int widthp2;
    int heightp2;

    public int[] localImageBuffer;
    public MemoryImageSource localImageSource;    
    //public Image localImage;
    
    public ComplexA imagedata;
    
    ImageConsumer notifyThis = null;
    public ComplexImageMapping map = new ComplexImageMapping();
    public void ComplexImage()
    {
    }
    
    public int closestPower2(int p)
    {
        double l = Math.log(p)/Math.log(2);
        return (int)Math.ceil(Math.exp(Math.log(2)*Math.ceil(l)));        
    }
    
    public ComplexImage(int wi,int he)
    {
        w = wi;
        h = he;
        widthp2 = closestPower2(wi);
        heightp2 = closestPower2(he);  
        imagedata = new ComplexA(w*h);
        localImageBuffer = new int[w*h];
        localImageSource = new MemoryImageSource(w,h,localImageBuffer,0,w);
        localImageSource.setAnimated(true);
        localImageSource.setFullBufferUpdates(false);
        //localImage = createImage(localImageSource);
    }
    public void copyImage(ComplexImage ci)
    { // assumes same size --- getting to close to deadline to throw excpetions
        for(int x=0;x<w;x++)
			{
					for(int y=0;y<h;y++)
					{
							imagedata.r[y*w+x]=ci.imagedata.r[y*w+x];
							imagedata.i[y*w+x]=ci.imagedata.i[y*w+x];							
							localImageBuffer[y*w+x]=ci.localImageBuffer[y*w+x];
					}
			}			
    }
    public void grabImage(Image im,ImageProducer ip)
    {
        PixelGrabber pg = new PixelGrabber(im,0,0,-1,-1,true);
        try{
        pg.grabPixels();
        } catch (InterruptedException ex){ex.printStackTrace();}
        localImageBuffer = (int[])pg.getPixels();
        
        // now essentially do the constructor's job
        w = pg.getWidth();
        h = pg.getHeight();
        widthp2 = closestPower2(w);
        heightp2 = closestPower2(h);
        localImageSource = new MemoryImageSource(w,h,localImageBuffer,0,w);
        localImageSource.setAnimated(true);
        localImageSource.setFullBufferUpdates(false);
        // now reinterpret the localImageBuffer as a complex array
        imagedata = new ComplexA(w*h);
        int x,y;
        for( x=0;x<w;x++)
        {
            for( y=0;y<h;y++)
            {
                //imagedata.r[y*w+x] = localImageBuffer[y*w+x] & 255;
                int i = localImageBuffer[y*w+x];
                imagedata.r[y*w+x] = ((i & 255) + ((i & (255 << 8)) >>> 8) + ((i & (255 << 16)) >>> 16))/3;
                imagedata.i[y*w+x] = 0;             
            }
        }
        System.out.println("Completed loading image");
    }
    public ComplexImage(ComplexImage ci)
    {
        w = ci.w;
        h = ci.h;
        widthp2 = ci.widthp2;
        heightp2 = ci.heightp2;
        imagedata = new ComplexA(w*h);
        localImageBuffer = new int[w*h];	
        localImageSource = new MemoryImageSource(w,h,localImageBuffer,0,w);
        localImageSource.setAnimated(true);
        localImageSource.setFullBufferUpdates(false);
        //localImage = createImage(localImageSource);
        copyImage(ci);			         
    }
    public ComplexImage(){}
    public void setDimensions(int wi,int he)
    {
        w = wi;
        h = he;
        widthp2 = closestPower2(wi);   /// argh BUG
        heightp2 = closestPower2(he);
        imagedata = new ComplexA(w*h);
        localImageBuffer = new int[w*h];	
        localImageSource = new MemoryImageSource(w,h,localImageBuffer,0,w);
        localImageSource.setAnimated(true);
        localImageSource.setFullBufferUpdates(false);
        complete = false;
        //localImage = createImage(localImageSource);
    }
    
    public void setColorModel(ColorModel cm)
    {
            // ignore this?
            System.out.println("ComplexImage as ImageConsumer: cm set as "+cm);
    }
    public void setHints(int i)
    {
            // ignor this
            System.out.println("ComplexImage as ImageConsumer: hints set as "+i);
    }
    public void setProperties(Hashtable h)
    {       // ignor this
            System.out.println("ComplexImage as ImageConsumer: hash set as "+h);    
    }
    public void setPixels(int x,int y,int xw, int yh, ColorModel cm, int pixels[],int off,int scan)
    {
        // need to get the magnitude of all the pixels and jam them in there
        System.out.println("ComplexImage as ImageConsumer: setPixels(int)");
        for(int sx=x;sx<xw;sx++)
        {
                for(int sy=y;sy<yh;sy++)
                {
                        int pixel = pixels[sy*scan+sx+off];
                        float m = cm.getRed(pixel) + cm.getBlue(pixel) + cm.getGreen(pixel);
                        m/=3;
                        imagedata.r[sx*w+sy] = m;
                        imagedata.i[sx*w+sy] = 0;
                }
        }
    }
    public void setPixels(int x,int y,int xw, int yh, ColorModel cm, byte pixels[],int off,int scan)
    {
        // need to get the magnitude of all the pixels and jam them in there
        for(int sx=x;sx<xw;sx++)
        {
                for(int sy=y;sy<yh;sy++)
                {
                        int pixel = (int)pixels[sy*scan+sx+off];
                        float m = cm.getRed(pixel) + cm.getBlue(pixel) + cm.getGreen(pixel);
                        m/=3;
                        imagedata.r[sx*w+sy] = m;
                        imagedata.i[sx*w+sy] = 0;
                }
        }
    }
    public void imageComplete(int i)
    {
        complete = true;               
    }
    
    public void setPixel(int x,int y, Complex c)
    {
        imagedata.setComplex(y*w+x,c);
    }
    
    public Complex getPixel(int x,int y)
    {
        return(imagedata.getComplex(y*w+x));
    }
    
    public void setPixelMagnitude(int x, int y, float m)
    {
        /*imagedata[y*w+x].setMag(m);
        
        // change this pixel (?)
        if (notifyThis!=null)
        {
            int[] i = new int[1];
            Complex c;
            map.sample(c=imagedata[y*w+x]);
            i[0] = map.perform(c); // problem with screen mapping here!!!!!!!!1
            notifyThis.setPixels(x,y,1,1,ColorModel.getRGBdefault(),i,0,1);
        }*/
    }
    
    public void rotatePixelPhase(int x, int y, float r)
    {
        //imagedata[y*w+x].rotPhase(r);
    }
    
    public int[] allocateIntArray()
    {
    	return new int[w*h];   	
    }

	public void stuffIntArray(int[] ar)
	{
		if (map.needsPreRun())
		{
			map.initialise();
			for(int x=0;x<w;x++)
			{
					for(int y=0;y<h;y++)
					{
							map.sample(imagedata.getComplex(y*w+x));
					}
			}			
			map.complete();
		}		
		for(int x=0;x<w;x++)
		{
				for(int y=0;y<h;y++)
				{
						ar[y*w+x]= map.perform(imagedata.getComplex(y*w+x));
				}
		}					
	}
	public void stuffIntArray(int[] ar,int lx,int ly,int mx,int my)
	{
		if (map.needsPreRun())
		{
//			map.initialise(); // since this is a partial update
			for(int x=lx;x<mx;x++)
			{
					for(int y=ly;y<my;y++)
					{
							map.sample(imagedata.getComplex(y*w+x));
					}
			}			
//			map.complete();
		}		
		for(int x=lx;x<mx;x++)
		{
				for(int y=ly;y<my;y++)
				{
						ar[y*w+x]= map.perform(imagedata.getComplex(y*w+x));
				}
		}					
	}
	
	
	public void fFFT(ComplexImage output, TwoDimMappingFunction mapper)
	{
		// the big one. - space to frequceny
		
		// this is going to be inplace
		// so copy our array over into a nicely formatted 
		// large real array --- just like numerical recipes wants us
		// and then copy it back into output
		
		// maybe we need to do this in a thread and provide feedback to
		// the user...
		
		float[] data = new float[widthp2*heightp2*2];
		
		for(int x=0;x<widthp2;x++)
		{
				for(int y=0;y<heightp2;y++)
				{
					//if outside window then just zero
					if ( (y<(heightp2-h)/2) || (y>=(heightp2+h)/2) || 
					      (x<(widthp2-w)/2) || (x>=(widthp2+w)/2) )
					{
						data[y*widthp2*2+x*2]   = 0;
						data[y*widthp2*2+x*2+1] = 0;
					}
					else
					{					
						float z = mapper.window( (x-widthp2/2.0)/(widthp2/2.0) , (y-heightp2/2.0)/(heightp2/2.0));
						int index = (x - (widthp2-w)/2) + (y - (heightp2-h)/2) * w;
						//System.out.println(index);
						data[y*widthp2*2+x*2]   = z*imagedata.r[ index ];
						data[y*widthp2*2+x*2+1] = z*imagedata.i[ index ];						
					}
				}
		}		
		
		// here comes the pure numrep
		
		//void fourn(float data[], unsigned long nn[], int ndim, int isign)
		//{
		int idim;
		int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
		int ibit,k1,k2,n,nprev,nrem,ntot;
		float tempi,tempr;
		double theta,wi,wpi,wpr,wr,wtemp; //Double precision for trigonometric recurences
		int isign = 1;
		
		int[] nn = new int[2];
		nn[1] = heightp2;
		nn[0] = widthp2;
		
		for (ntot=1,idim=1;idim<=2;idim++) //Compute total number of complex val-
												//ues. 
		ntot *= nn[idim-1];
		nprev=1;
		for (idim=2;idim>=1;idim--) { //Main loop over the dimensions.
			n=nn[idim-1];
			nrem=ntot/(n*nprev);
			ip1=nprev << 1;
			ip2=ip1*n;
			ip3=ip2*nrem;
			i2rev=1;
			for (i2=1;i2<=ip2;i2+=ip1) { //This is the bit-reversal section of the
									//routine. 
				if (i2 < i2rev) {
					for (i1=i2;i1<=i2+ip1-2;i1+=2) {
						for (i3=i1;i3<=ip3;i3+=ip2) {
							i3rev=i2rev+i3-i2;
							//System.out.println(" Before swap :"+data[i3-1]+" "+data[i3rev-1]);
							//swap(data[i3-1],data[i3rev-1]);
							float t =  data[i3-1]; data[i3-1]=data[i3rev-1]; data[i3rev-1]=t;
							//System.out.println(" After "+data[i3-1]+" "+data[i3rev-1]);
							t =  data[i3]; data[i3]=data[i3rev]; data[i3rev]=t;
//							swap(data[i3],data[i3rev]);
						}
					}
				}
				ibit=ip2 >>> 1;
				while ((ibit >= ip1) && (i2rev > ibit)) {	
				i2rev -= ibit;
				ibit >>>= 1;
			}
			i2rev += ibit;
		}
		ifp1=ip1; //Here begins the Danielson-Lanczos sec-
				  //tion of the routine. 
		while (ifp1 < ip2) {
			ifp2=ifp1 << 1;
			theta=isign*6.28318530717959/(ifp2/ip1); //Initialize for the trig. recur-
												//rence. 
			wtemp=Math.sin(0.5*theta);
			wpr = -2.0*wtemp*wtemp;
			wpi=Math.sin(theta);
			wr=1.0;
			wi=0.0;
			for (i3=1;i3<=ifp1;i3+=ip1) {
				for (i1=i3;i1<=i3+ip1-2;i1+=2) {
					for (i2=i1;i2<=ip3;i2+=ifp2) {
						k1=i2; //Danielson-Lanczos formula:
						k2=k1+ifp1;
						tempr=(float)wr*data[k2-1]-(float)wi*data[k2];
						tempi=(float)wr*data[k2]+(float)wi*data[k2-1];
						data[k2-1]=data[k1-1]-tempr;
						data[k2]=data[k1]-tempi;
						data[k1-1] += tempr;
						data[k1] += tempi;
					}	
				}
				wr=(wtemp=wr)*wpr-wi*wpi+wr; //Trigonometric recurrence.
				wi=wi*wpr+wtemp*wpi+wi;
			}
			ifp1=ifp2;
		}
		nprev *= n;
		
		}
		
		// output must be widthp2 by heightp2
		
		for(int x=0;x<widthp2;x++)
		{
				for(int y=0;y<heightp2;y++)
				{
					//if outside window then just zero
					/*if ( (y<(heightp2-h)/2) || (y>=(heightp2+h)/2) || 
					      (x<(widthp2-w)/2) || (x>=(widthp2+w)/2) )
					{
						data[y*w*2+x*2]   = 0;
						data[y*w*2+x*2+1] = 0;
					}
					else
					{					
						float z = mapper.window( (x-widthp2/2)/(widthp2/2) , (y-heightp2/2)/(heightp2/2));
						int index = (x - (widthp2-w)/2) + (y - (heightp2-h)/2) * w;
						data[y*w*2+x*2]   = z*imagedata.r[ index ];
						data[y*w*2+x*2+1] = z*imagedata.i[ index ];						
					}*/
					
					output.imagedata.r[y*widthp2+x] = (float)(data[((y + heightp2/2) % (heightp2))*widthp2*2+((x+widthp2/2) % widthp2)*2]/Math.sqrt(widthp2*heightp2));
					output.imagedata.i[y*widthp2+x] = (float)(data[((y + heightp2/2) % (heightp2))*widthp2*2+((x+widthp2/2) % widthp2)*2 + 1]/Math.sqrt(widthp2*heightp2));					
					//output.imagedata.r[y*widthp2+x] = data[y*widthp2*2 + x*2];
					//output.imagedata.i[y*widthp2+x] = data[y*widthp2*2 + x*2+1];
					
    				//System.out.println("Sig : "+data[y*widthp2*2+x*2]+" "+data[y*widthp2*2+x*2+1]);	
				}
				
		}		
				
	}
	public void iFFT(ComplexImage output, TwoDimMappingFunction mapper)
	{
		// the big one. - freq to space
		
		// this is going to be inplace
		// so copy our array over into a nicely formatted 
		// large real array --- just like numerical recipes wants us
		// and then copy it back into output
		
		// maybe we need to do this in a thread and provide feedback to
		// the user...
		
		float[] data = new float[widthp2*heightp2*2];
		
		for(int x=0;x<widthp2;x++)
		{
				for(int y=0;y<heightp2;y++)
				{
				    data[y*widthp2*2+x*2]   = imagedata.r[ ((y + heightp2/2) % (heightp2))*widthp2 + ((x+widthp2/2) % widthp2)];
					data[y*widthp2*2+x*2+1] = imagedata.i[ ((y + heightp2/2) % (heightp2))*widthp2 + ((x+widthp2/2) % widthp2)];		
    				//System.out.println("Sig : "+data[y*widthp2*2+x*2]+" "+data[y*widthp2*2+x*2+1]);
				}
		}		
		
		// here comes the pure numrep
		
		//void fourn(float data[], unsigned long nn[], int ndim, int isign)
		//{
		int idim;
		int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
		int ibit,k1,k2,n,nprev,nrem,ntot;
		float tempi,tempr;
		double theta,wi,wpi,wpr,wr,wtemp; //Double precision for trigonometric recurences
		int isign = -1;
		
		int[] nn = new int[2];
		nn[1] = heightp2; // is this the right way round?
		nn[0] = widthp2;
		
		for (ntot=1,idim=1;idim<=2;idim++)//Compute total number of complex val-
			ntot *= nn[idim-1];
		nprev=1;
		for (idim=2;idim>=1;idim--) { //Main  loop over the dimensions.
			n=nn[idim-1];
			nrem=ntot/(n*nprev);
			ip1=nprev << 1;
			ip2=ip1*n;
			ip3=ip2*nrem;
			i2rev=1;
			for (i2=1;i2<=ip2;i2+=ip1) { //This is the bit-reversal section of the
									//routine. 
				if (i2 < i2rev) {
					for (i1=i2;i1<=i2+ip1-2;i1+=2) {
						for (i3=i1;i3<=ip3;i3+=ip2) {
							i3rev=i2rev+i3-i2;
							float t =  data[i3-1]; data[i3-1]=data[i3rev-1]; data[i3rev-1]=t;
							//System.out.println(" After "+data[i3-1]+" "+data[i3rev-1]);
							t =  data[i3]; data[i3]=data[i3rev]; data[i3rev]=t;

//							swap(data[i3-1],data[i3rev-1]);
	//						swap(data[i3],data[i3rev]);
						}
					}
				}
				ibit=ip2 >>> 1;
				while ((ibit >= ip1) && (i2rev > ibit)) {	
    				i2rev -= ibit;
	    			ibit >>>= 1;
		    	}
			i2rev += ibit;
		}
		ifp1=ip1; //Here begins the Danielson-Lanczos sec-
				  //tion of the routine. 
		while (ifp1 < ip2) {
			ifp2=ifp1 << 1;
			theta=isign*6.28318530717959/(ifp2/ip1); //Initialize for the trig. recur-
												//rence. 
			wtemp=Math.sin(0.5*theta);
			wpr = -2.0*wtemp*wtemp;
			wpi=Math.sin(theta);
			wr=1.0;
			wi=0.0;
			for (i3=1;i3<=ifp1;i3+=ip1) {
				for (i1=i3;i1<=i3+ip1-2;i1+=2) {
					for (i2=i1;i2<=ip3;i2+=ifp2) {
						k1=i2; //Danielson-Lanczos formula:
						k2=k1+ifp1;
						tempr=(float)wr*data[k2-1]-(float)wi*data[k2];
						tempi=(float)wr*data[k2]+(float)wi*data[k2-1];
						data[k2-1]=data[k1-1]-tempr;
						data[k2]=data[k1]-tempi;
						data[k1-1] += tempr; //System.out.println(k1-1);
						data[k1] += tempi;  //System.out.println(k1+" "+data[k1]);
					}	
				}
				wr=(wtemp=wr)*wpr-wi*wpi+wr; //Trigonometric recurrence.
				wi=wi*wpr+wtemp*wpi+wi;
			}
			ifp1=ifp2;
		}
		nprev *= n;
		
		}
		
		// output must be widthp2 by heightp2
		
		for(int x=0;x<output.w;x++)
		{
				for(int y=0;y<output.h;y++)
				{
					//if outside window then just zero
					/*if ( (y<(heightp2-h)/2) || (y>=(heightp2+h)/2) || 
					      (x<(widthp2-w)/2) || (x>=(widthp2+w)/2) )
					{
						data[y*w*2+x*2]   = 0;
						data[y*w*2+x*2+1] = 0;
					}
					else
					{					
						float z = mapper.window( (x-widthp2/2)/(widthp2/2) , (y-heightp2/2)/(heightp2/2));
						int index = (x - (widthp2-w)/2) + (y - (heightp2-h)/2) * w;
						data[y*w*2+x*2]   = z*imagedata.r[ index ];
						data[y*w*2+x*2+1] = z*imagedata.i[ index ];						
					}*/
					
					int index = (y+(heightp2-output.h)/2)*widthp2*2 + (x+(widthp2-output.w)/2)*2;
                   // System.out.println(index+" "+data[index]+" "+data[index+1]);
					output.imagedata.r[y*output.w+x] = (float)(data[index]/Math.sqrt(widthp2*heightp2));
					output.imagedata.i[y*output.w+x] = (float)(data[index+1]/Math.sqrt(widthp2*heightp2));					
					
				}
		}		
				
	}
	
}