#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "genUtil.h"
#include "mrcImage.h"
#include "Vector.h"
#include "./lmrcFFTInfo.h"

floatVector*
lmrcFSInfoXAxisMag(mrcImage* fft)
{
    floatVector* fv;
    float iX, iY;
    double mag;
    unsigned long size, i;

    size = fft->HeaderN.x/2+1;

	fv = floatVectorInit(NULL, size);
    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
    }

    for(iX=0.0, i=0; iX>=-fft->HeaderN.x/2.0; iX--, i++) {
        mrcPixelDataGet(fft, iX, 0.0, 0.0, &mag, mrcPixelPow, mrcPixelHowNearest);
        fv->data[i] = mag;
    }
    return fv;
}

floatVector*
lmrcFSInfoXAxisPhase(mrcImage* fft)
{
    floatVector* fv;
    float iX, iY;
    double phase;
    unsigned long size, i;

    size = fft->HeaderN.x/2+1;
	fv = floatVectorInit(NULL, size);

    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
    }

    for(iX=0.0, i=0; iX>=-fft->HeaderN.x/2.0; iX--, i++) {
        mrcPixelDataGet(fft, iX, 0.0, 0.0, &phase, mrcPixelPhase, mrcPixelHowNearest);
        fv->data[i] = phase;
    }
    return fv;
}

floatVector*
lmrcFSInfoYAxisMagAverage(mrcImage* fft)
{
    floatVector* fv;
    float iX, iY;
    double mag;
    unsigned long size, i;

    size = fft->HeaderN.y/2+1;
	fv = floatVectorInit(NULL, size);

    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
    }

    for(iY=0.0, i=0; iY>=-fft->HeaderN.y/2.0; iY--, i++) {
		fv->data[i] = 0.0;
    	for(iX=0.0; iX<fft->HeaderN.x/2.0; iX++) {
       		mrcPixelDataGet(fft, iX, iY, 0.0, &mag, mrcPixelMag, mrcPixelHowNearest);
        		fv->data[i] += mag;
			}
		}
		fv->data[i] /= (fft->HeaderN.x/2.0);
    }
    return fv;
}

floatVector*
lmrcFSInfoYAxisMag(mrcImage* fft)
{
    floatVector* fv;
    float iX, iY;
    double mag;
    unsigned long size, i;

    size = fft->HeaderN.y/2+1;
	fv = floatVectorInit(NULL, size);

    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
    }

    for(iY=0.0, i=0; iY>=-fft->HeaderN.y/2.0; iY--, i++) {
        mrcPixelDataGet(fft, 0.0, iY, 0.0, &mag, mrcPixelMag, mrcPixelHowNearest);
        fv->data[i] = mag;
    }
    return fv;
}

floatVector*
lmrcFSInfoYAxisPhase(mrcImage* fft)
{
    floatVector* fv;
    float iX, iY;
    double phase;
    unsigned long size, i;

    size = fft->HeaderN.y/2+1;
	fv = floatVectorInit(NULL, size);

    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
    }

    for(iY=0.0, i=0; iY>=-fft->HeaderN.y/2.0; iY--, i++) {
        mrcPixelDataGet(fft, 0.0, iY, 0.0, &phase, mrcPixelPhase, mrcPixelHowNearest);
        fv->data[i] = phase;
    }
    return fv;
}

floatVector*
lmrcFSInfoScatteringAngularDistribution(mrcImage* fft)
{
    floatVector* fv;
    float iX, iY;
    double re, im, pow;
    int flag, i, R;
    float mul, mul2;
    float* data;
    unsigned long size;

    if(fft->HeaderN.x<fft->HeaderN.y) {
        flag = 0;
        size = fft->HeaderN.x/2+1;
        mul = ((float)fft->HeaderN.x)/((float)fft->HeaderN.y);
    } else {
        flag = 1;
        size = fft->HeaderN.y/2+1;
        mul = ((float)fft->HeaderN.y)/((float)fft->HeaderN.x);
    }
    mul2 = mul*mul;

	fv = floatVectorInit(NULL, size);
    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
    }


    if(flag==0) {
        for(iX=-fft->HeaderN.x/2.0+1; iX<fft->HeaderN.x/2.0; iX++) {
            for(iY=-fft->HeaderN.y/2.0+1; iY<fft->HeaderN.y/2.0; iY++) {
                mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
                R = (int)(sqrt(iX*iX+iY*iY*mul2)+0.5);
				if(R<size) {
                	fv->data[R] += pow;
				}
            }
        }
    } else {
        for(iX=-fft->HeaderN.x/2.0+1; iX<=fft->HeaderN.x/2.0; iX++) {
            for(iY=-fft->HeaderN.y/2.0+1; iY<fft->HeaderN.y/2.0; iY++) {
                mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
                R = (int)(sqrt(iX*iX*mul2+iY*iY)+0.5);
				if(R<size) {
                	fv->data[R] += pow;
				}
            }
        }
    }
    return fv;
}

floatVector*
lmrcFSInfoScatteringAngularDistributionAverage(mrcImage* fft)
{
    floatVector* fv;
    floatVector* count;
    float iX, iY;
    double re, im, pow;
    int flag, i, R;
    float mul, mul2;
    float* data;
    unsigned long size;

    if(fft->HeaderN.x<fft->HeaderN.y) {
        flag = 0;
        size = fft->HeaderN.x/2+1;
        mul = ((float)fft->HeaderN.x)/((float)fft->HeaderN.y);
    } else {
        flag = 1;
        size = fft->HeaderN.y/2+1;
        mul = ((float)fft->HeaderN.y)/((float)fft->HeaderN.x);
    }
    mul2 = mul*mul;
	fv    = floatVectorInit(NULL, size);
	count = floatVectorInit(NULL, size);
	DEBUGPRINT3("%d %d %d\n", count->size, fv->size, size);
    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
        count->data[i] = 0.0;
    }

    if(flag==0) {
        for(iX=-fft->HeaderN.x/2.0+1; iX<fft->HeaderN.x/2.0; iX++) {
            for(iY=-fft->HeaderN.y/2.0+1; iY<fft->HeaderN.y/2.0; iY++) {
                mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
                R = (int)(sqrt(iX*iX+iY*iY*mul2)+0.5);
				if(R<size) {
                	fv->data[R] += pow;
                	count->data[R] ++;
				}
            }
        }
    } else {
        for(iX=-fft->HeaderN.x/2.0+1; iX<fft->HeaderN.x/2.0; iX++) {
            for(iY=-fft->HeaderN.y/2.0+1; iY<fft->HeaderN.y/2.0; iY++) {
                mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
                R = (int)(sqrt(iX*iX*mul2+iY*iY)+0.5);
				if(R<size) {
                	fv->data[R] += pow;
                	count->data[R] ++;
				}
            }
        }
    }
    for(i=0; i<size; i++) {
		if(0!=count->data[i]) {
        	fv->data[i] /= count->data[i] ;
		} else {
			fv->data[i] = 0;
		}
    }
    return fv;
}

floatVector*
lmrcFSInfoScatteringAngularDistributionSD(mrcImage* fft)
{
    floatVector* fv;
    floatVector* ave;
    floatVector* ave2;
    floatVector* count;
    float iX, iY;
    double re, im, pow;
    int flag, i, R;
    float mul, mul2;
    float* data;
    unsigned long size;

    if(fft->HeaderN.x<fft->HeaderN.y) {
        flag = 0;
        size = fft->HeaderN.x/2+1;
        mul = ((float)fft->HeaderN.x)/((float)fft->HeaderN.y);
    } else {
        flag = 1;
        size = fft->HeaderN.y/2+1;
        mul = ((float)fft->HeaderN.y)/((float)fft->HeaderN.x);
    }
    mul2 = mul*mul;
	fv    = floatVectorInit(NULL, size);
	ave   = floatVectorInit(NULL, size);
	ave2  = floatVectorInit(NULL, size);
	count = floatVectorInit(NULL, size);
	DEBUGPRINT3("%d %d %d\n", count->size, fv->size, size);
    /* Clear */
    for(i=0; i<size; i++) {
        fv->data[i] = 0.0;
        ave->data[i] = 0.0;
        ave2->data[i] = 0.0;
        count->data[i] = 0.0;
    }

    if(flag==0) {
        for(iX=-fft->HeaderN.x/2.0+1; iX<fft->HeaderN.x/2.0; iX++) {
            for(iY=-fft->HeaderN.y/2.0+1; iY<fft->HeaderN.y/2.0; iY++) {
                mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
                R = (int)(sqrt(iX*iX+iY*iY*mul2)+0.5);
				if(R<size) {
                	ave->data[R] += pow;
                	ave2->data[R] += pow*pow;
                	count->data[R] ++;
				}
            }
        }
    } else {
        for(iX=-fft->HeaderN.x/2.0+1; iX<fft->HeaderN.x/2.0; iX++) {
            for(iY=-fft->HeaderN.y/2.0+1; iY<fft->HeaderN.y/2.0; iY++) {
                mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
                R = (int)(sqrt(iX*iX*mul2+iY*iY)+0.5);
				if(R<size) {
                	ave->data[R] += pow;
                	ave2->data[R] += pow*pow;
                	count->data[R] ++;
				}
            }
        }
    }
    for(i=0; i<size; i++) {
		if(0!=count->data[i]) {
        	ave->data[i] /= count->data[i] ;
        	ave2->data[i] /= count->data[i] ;
			fv->data[i] = ave2->data[i] - (ave->data[i]*ave->data[i]);
		} else {
			ave->data[i] = 0;
			ave2->data[i] = 0;
			fv->data[i] = 1.0;
		}
    }
    return fv;
}

floatVector*
lmrcFSInfoSpacing(mrcImage* fft)
{
    floatVector* fv;
    int flag, i, R; 
	unsigned long size;

    if(fft->HeaderN.x<fft->HeaderN.y) {
        flag = 0;
        size = fft->HeaderN.x/2+1;
    } else {
        flag = 1;
        size = fft->HeaderN.y/2+1;
    }
	fv    = floatVectorInit(NULL, size);
	DEBUGPRINT3("%d %d %d\n", count->size, fv->size, size);
	if(fft->HeaderN.x<fft->HeaderN.y) {
    	for(i=0; i<size; i++) {
       		fv->data[i] = i/(fft->HeaderN.x*fft->HeaderLength.x);
    	}
	} else {
    	for(i=0; i<size; i++) {
       		fv->data[i] = i/(fft->HeaderN.y*fft->HeaderLength.y);
    	}
	}
    return fv;
}
