/*
# %M% %Y% %I%
# The latest update : %G% at %U%
#
#%Z% lmrcImageCorrelation ver %I%
#%Z% Created by 
#%Z%
#%Z% Usage : lmrcImageCorrelation
#%Z% Attention
#%Z%
*/
static char __sccs_id[] = "%Z%lmrcImageCorrelation ver%I%; Date:%D% %Z%";
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#undef DEBUG
#include "genUtil.h"
#include "../inc/mrcImage.h"       
#include "./lmrcImageCorrelation.h"
#include "lDescartesIntoPolar.h"
#include "./lmrcImageTransformDescartesIntoPolar.h"

/*
	lmrcImageEuclidDistanceCalc: 
		mode: 0 : Sum of SQR
		      1 : Distance 
*/

double 
lmrcImageEuclidDistanceCalc(mrcImage* in1, mrcImage* in2, int mode)
{
	double data1, data2;
	double data;
	mrcImageParaTypeReal x, y, z;

	if(in1->HeaderN.x!=in2->HeaderN.x
	 ||in1->HeaderN.y!=in2->HeaderN.y
	 ||in1->HeaderN.z!=in2->HeaderN.z) {
		fprintf(stderr, "Attention in lmrcImageEuclidDistanceCalc\n");
		fprintf(stderr, "size is different (%f, %f, %f) != (%f, %f, %f)\n", 
			in1->HeaderN.x, in1->HeaderN.y, in1->HeaderN.z,
			in2->HeaderN.x, in2->HeaderN.y, in2->HeaderN.z);
	}
	data = 0.0;
	for(z=0; z<in1->HeaderN.z; z++) {
	for(y=0; y<in1->HeaderN.y; y++) {
	for(x=0; x<in1->HeaderN.x; x++) {
		 mrcPixelDataGet(in1 ,x ,y ,z ,&data1 ,mrcPixelRePart ,mrcPixelHowNearest);
		 mrcPixelDataGet(in2 ,x ,y ,z ,&data2 ,mrcPixelRePart ,mrcPixelHowNearest);
		 data += SQR(data2-data1);
	}
	}
	}
	if(data!=0 && mode==1) {
		data = sqrt(data);
	} 
	return data;
}

/*

*/
void
lmrcImageAutoRotationCorrelationRotationalCrossCorrelationFunction(
	mrcImage* out,  mrcImage* cor, 
    mrcImage* in, mrcImage* ref, 
	lmrcImageAutoRotationCorrelationInfo* linfo, int mode)
{
	mrcImage inPolar;
	mrcImage refPolar;
	mrcImage Cor;
	lmrcImageTransformDescartesIntoPolarInfo llinfoIn;
	lmrcImageTransformDescartesIntoPolarInfo llinfoRef;
	lDescartesIntoPolarInfo                  llinfo2In;
	lDescartesIntoPolarInfo                  llinfo2Ref;

	llinfoIn.dr     = ref->HeaderLength.x; 
	llinfoIn.dphi   = 2*M_PI/linfo->nRot;
	llinfoIn.dtheta = 2*M_PI;
	llinfoIn.flagImageCentreIsGravityCentre = 1;
	llinfoIn.flagDescartesIntoPolarInfo    = 0;
	llinfoIn.flagrWeight                    = 1;
	llinfoRef = llinfoIn;
	lmrcImageTransformDescartesIntoPolar(&inPolar,  in,  &llinfoIn,  &llinfo2In, 0);  
	lmrcImageTransformDescartesIntoPolar(&refPolar, ref, &llinfoRef, &llinfo2Ref, 0);  
	lmrcImageCorrelation(&Cor, &inPolar, &refPolar, mode);
	linfo->corInfo.mode = meanOfAll;
	lmrcImageInformation(&(linfo->corInfo), &Cor);
	linfo->max      = linfo->corInfo.max;
	linfo->maxP     = linfo->corInfo.maxCoord;
	linfo->maxTheta = linfo->maxP.x*llinfoIn.dphi;
	mrcImageFree(&inPolar, 0);
	mrcImageFree(&refPolar, 0);
	mrcImageFree(&Cor, 0);
}

/*
	lmrcImageAutoRotationCorrelation
*/

void
lmrcImageAutoRotationCorrelation(mrcImage* out,  mrcImage* cor, 
                                 mrcImage* in, mrcImage* ref, 
								 lmrcImageAutoRotationCorrelationInfo* linfo, int mode)
{
	int i;
	int ntheta;
	float dtheta;
	float theta;
	float minTheta;
	float maxTheta;
	float maxOld;
	float current;
	mrcImageParaTypeRealCoord currentp;
	mrcImage rotproj;
	mrcImage inFFT;
	mrcImage inImage;
	double data;
	mrcImageParaTypeReal x, y, z;

	DEBUGPRINT("Start lmrcImageAutoRotationCorrelation---\n");
	/* 
		Real Space to Fourier Space 
			Set inImage and inFFT 
	*/
	if(IsImage(in, "IsImage", 0)) {
		lmrcImageFFT(&inFFT, in, 0);
		inImage = *in;
	} else if(IsFT(in, "IsFT", 0)) {
		lmrcImageFFT(&inImage, in, 0);
		inFFT = *in;
	} else {
		fprintf(stderr, "Not supported mode: %ld\n", in->HeaderMode);
	}
	
	/* 
		Calculate corralation between in and rotated ref 
			Set search area
	*/
	if(linfo->flagRestrictionArea == 1) {
		minTheta = linfo->thetaMin;
		maxTheta = linfo->thetaMax;
	} else {
		minTheta = 0;                       /* minimum theta */  
		maxTheta = 2.0*M_PI;       			/* maximum theta */
	}
	ntheta = linfo->nRot;         		  /* Number of Rotation */
	dtheta = (maxTheta-minTheta)/ntheta;  /* Delta of Rotation */  
	linfo->corInfo.mode = meanOfAll;

	switch(linfo->Method) {
		case lmrcImageAutoRotationCorrelationMethodRotationAndCorrelation:
		case lmrcImageAutoRotationCorrelationMethodRotationAndCorrelationSSDA: {
			current  = -1e30; linfo->max = maxOld = -1e30; /* Initial */ 
			for(i=0; i<linfo->iter; i++) {
				for(theta=minTheta; theta<maxTheta; theta+=dtheta) {
					/* Rotate ref to rotproj */ 
					lmrcImageRotation2DPeriodicBoundary(&rotproj, ref, theta, mrcPixelHowLinear); 
					switch(linfo->Method) {
						case lmrcImageAutoRotationCorrelationMethodRotationAndCorrelation: {
							lmrcImageCorrelation(cor, &inFFT, &rotproj, mode); 
							if(linfo->flagXshiftOnly) {
								current = -1e30;
								for(x=0; x<cor->HeaderN.x; x++) {
									mrcPixelDataGet(cor, x, 0, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
									if(current<data) {
										current    = data;
										currentp.x = x; 
										currentp.y = 0.0;
										currentp.z = 0.0;
									}
								}
							} else if(linfo->flagNoShift) {
								currentp.x = currentp.y = currentp.z = 0.0;
								mrcPixelDataGet(cor, x, 0, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
								current = data;   
							} else {
								lmrcImageInformation(&(linfo->corInfo), cor);
								current  = linfo->corInfo.max;
								currentp = linfo->corInfo.maxCoord;
							}
							break;
						}
						case lmrcImageAutoRotationCorrelationMethodRotationAndCorrelationSSDA: {
							lmrcImageSSDAAutomaticThresholdAndFineSearch(&inImage, &rotproj, &linfo->SSDAInfo, mode);
							current = 1/(1+linfo->SSDAInfo.Error);
							currentp = linfo->SSDAInfo.optimump;
							DEBUGPRINT4("Angle: %15.6f Correlation: %15.6g at %15.6f %15.6f :: CurrentMax :: Linear \n", 
								theta*DEGREE, current, currentp.x, currentp.y);
							break;
						}
						default: {
							fprintf(stderr, "Not supported Method in lmrcImageAutoRotationCorrelation: %d\n", 
								linfo->Method);
							exit(EXIT_FAILURE);
							break;
						}
					}
					if(linfo->max<current) { /* When bettern correlation */ 
						linfo->max      = current;
						linfo->maxTheta = theta;
						linfo->maxP     = currentp;
						DEBUGPRINT4("Angle: %15.6f Correlation: %15.6g at %15.6f %15.6f :: CurrentMax :: Linear \n", 
							theta*DEGREE, current, currentp.x, currentp.y);
					}
					mrcImageFree(&rotproj, 0);
					mrcImageFree(cor, 0);
				}
				if(fabs(maxOld+linfo->max)!=0) {
					if(fabs(maxOld-linfo->max)/fabs(maxOld+linfo->max)<1e-6) {
						break;
					} else {
						maxOld = linfo->max;
						minTheta = linfo->maxTheta - dtheta; 
						maxTheta = linfo->maxTheta + dtheta; 
						dtheta  /= 10.0;
					}
				} else {
					maxOld = linfo->max;
					minTheta = linfo->maxTheta - dtheta; 
					maxTheta = linfo->maxTheta + dtheta; 
					dtheta  /= 10.0;
				}
			}
			break;
		}
		case lmrcImageAutoRotationCorrelationMethodRotationalFunction: { /* Rotational Correlation Function */
			lmrcImageAutoRotationCorrelationRotationalCrossCorrelationFunction(out, cor, in, ref, linfo, mode);
			break;
		}
		default: {
			fprintf(stderr, "Not supported Method in lmrcImageAutoRotationCorrelation: %d\n", linfo->Method);
			exit(EXIT_FAILURE);
			break;
		}
	}
	/* Calculate final correlation map and correlation peak */
	lmrcImageRotation2DPeriodicBoundary(&rotproj, ref, linfo->maxTheta, mrcPixelHowCubicConv); 
	lmrcImageCorrelation(cor, &inFFT, &rotproj, mode);
	
	if(linfo->flagXshiftOnly) {
		linfo->max = -1e30;
		for(x=0; x<cor->HeaderN.x; x++) {
			mrcPixelDataGet(cor, x, 0, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
			if(linfo->max<data) {
				linfo->max    = data;
				linfo->maxP.x = x;
				linfo->maxP.y = 0.0;
				linfo->maxP.z = 0.0;
			}
		}
	} else if(linfo->flagNoShift) {
		mrcPixelDataGet(cor, x, 0, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
		linfo->max = data;
		linfo->maxP.x = linfo->maxP.y = linfo->maxP.z = 0.0;
	} else {
		linfo->corInfo.mode = meanOfAll;
		lmrcImageInformation(&(linfo->corInfo), cor);
		DEBUGPRINT4("Angle: %15.6f Correlation: %15.6f at %15.6f %15.6f :: MaxPoint:: Cubic\n", \
			linfo->maxTheta*DEGREE, linfo->corInfo.max, linfo->corInfo.maxCoord.x, linfo->corInfo.maxCoord.y); 
		linfo->max  = linfo->corInfo.max;
		linfo->maxP = linfo->corInfo.maxCoord;
	}
	mrcImageFree(&rotproj, 0);

	/* Calculate outImage from inImage */
	currentp.x = -linfo->maxP.x;
	currentp.y = -linfo->maxP.y;
	currentp.z = -linfo->maxP.z;
	lmrcImageShift(&rotproj, &inImage, currentp, mrcPixelHowCubicConv); 
	lmrcImageRotation2DPeriodicBoundary(out, &rotproj, -linfo->maxTheta, mrcPixelHowCubicConv); 
	linfo->EuclidDistance    = lmrcImageEuclidDistanceCalc(out, ref, 1);  
	linfo->LinearCorrelation = lmrcImageLinearCorrelation (out, ref, 1);  
	DEBUGPRINT1("EuclidDistance:    %f \n", linfo->EuclidDistance); 
	DEBUGPRINT1("LinearCorrelation: %f \n", linfo->LinearCorrelation); 

	mrcImageFree(&rotproj, 0);
	if(IsImage(in, "IsImage", 0)) {
		mrcImageFree(&inFFT, 0);
	} else if(IsFT(in, "IsFT", 0)) {
		mrcImageFree(&inImage, 0);
	} else {
		fprintf(stderr, "Not supported mode: %ld\n", in->HeaderMode);
	}
}

void
lmrcImageCorrelationRefinement(mrcImage* out, mrcImage* shift, mrcImage* in, mrcImage* ref, lmrcImageCorrelationRefinementInfo* info, long mode)
{
    mrcImage inFFT;
    mrcImage refFFT;
    mrcImage shiftFFT;
	long     flagIn  = 0;
	long     flagRef = 0;
	mrcImageParaTypeReal  x,  y,  z;
	mrcImageParaTypeReal ix, iy, iz;
	mrcImageParaTypeReal nx, ny, nz;
	mrcImageParaTypeReal sx, sy, sz;
	mrcImageParaTypeReal dx, dy, dz;
	mrcImageParaTypeReal X, Y, Z;
	double mag, phase;
	double rmag, rphase;
	double eval;
	double sum;
	double deltaPhase;
	double dp;

    if(IsImage(in, "", 0)&&IsImage(ref, "", 0)) {
        lmrcImageFFT(&inFFT,  in, 0);
        lmrcImageFFT(&refFFT, ref, 0);
		flagIn = flagRef = 1;
    } else if(IsFT(in, "", 0)&&IsImage(ref, "", 0)) {
        inFFT = *in;
        lmrcImageFFT(&refFFT, ref, 0);
		flagRef = 1;
    } else if(IsImage(in, "", 0)&&IsFT(ref, "", 0)) {
        lmrcImageFFT(&inFFT,  in, 0);
        refFFT = *ref;
		flagIn = 1;
    } else if(IsFT(in, "", 0)&&IsFT(ref, "", 0)) {
        inFFT = *in;
        refFFT = *ref;
    } else {
        fprintf(stderr, "Nut supported Mode (in %ld, ref %ld)\n", in->HeaderMode, ref->HeaderMode);
        exit(EXIT_FAILURE);
   	}

	if(1<in->HeaderN.x) {
		sx = info->shift.x - info->range.x; 
		nx = 2*info->range.x/info->step.x + 1;
	} else {
		sx = 0;
		nx = 1;
	}
	if(1<in->HeaderN.y) {
		sy = info->shift.y - info->range.y; 
		ny = 2*info->range.y/info->step.y + 1;
	} else {
		sy = 0;
		ny = 1;
	}
	if(1<in->HeaderN.z) {
		sz = info->shift.z - info->range.z; 
		nz = 2*info->range.z/info->step.z + 1;
	} else {
		sz = 0;
		nz = 1;
	}
	dx = info->step.z;
	dy = info->step.z;
	dz = info->step.z;

	mrcDefaultHeaderValueSet(out);
	out->HeaderN.x = nx;
	out->HeaderN.y = ny;
	out->HeaderN.z = nz;
	out->HeaderMode = mrcFloatImage;
	out->HeaderLength.x = dx*in->HeaderLength.x;
	out->HeaderLength.y = dy*in->HeaderLength.y;
	out->HeaderLength.z = dz*in->HeaderLength.z;
	mrcInit(out, NULL);

	for(ix=0, x=sx; ix<nx; ix++, x+=dx) {
	for(iy=0, y=sy; iy<ny; iy++, y+=dy) {
	for(iz=0, z=sz; iz<nz; iz++, z+=dz) {
		eval = sum = 0.0;
		for(X=0;                  X<=in->HeaderN.x/2.0; X++) {
		for(Y=-in->HeaderN.y/2.0; Y< in->HeaderN.y/2.0; Y++) {
		for(Z=-in->HeaderN.z/2.0; Z< in->HeaderN.z/2.0; Z++) {
			deltaPhase = -2*M_PI*x*X/in->HeaderN.x 
						 -2*M_PI*y*Y/in->HeaderN.y
						 -2*M_PI*z*Z/in->HeaderN.z;
			mrcPixelDataGet(&inFFT, X, Y, Z, &mag,    mrcPixelMag,   mrcPixelHowNearest);
			mrcPixelDataGet(&inFFT, X, Y, Z, &phase,  mrcPixelPhase, mrcPixelHowNearest);
			mrcPixelDataGet(&refFFT, X, Y, Z, &rmag,   mrcPixelMag,   mrcPixelHowNearest);
			mrcPixelDataGet(&refFFT, X, Y, Z, &rphase, mrcPixelPhase, mrcPixelHowNearest);
			switch(info->mode&0x0000000f) {
				case 0: { /* Normal Correlation */
					eval += (SQR(mag*cos(phase+deltaPhase) - rmag*cos(rphase)) 
							+SQR(mag*sin(phase+deltaPhase) - rmag*sin(rphase)));
					sum  += SQR(rmag);
					break;
				}
				case 1: { /* Phase Correlation */
					dp = NORMAL_PHASE(phase+deltaPhase - rphase);
					dp = MIN(dp, 2*M_PI - dp);
					eval += rmag*dp;
					sum  += rmag; 
					break;
				}
				case 2: { /* Weighted Phase Correlation */
					dp = NORMAL_PHASE(phase+deltaPhase - rphase);
					dp = MIN(dp, 2*M_PI - dp);
					eval += rmag*SQR(dp);
					sum  += rmag; 
					break;
				}
				default: {
					fprintf(stderr, "Not supported mode: %ld in lmrcImageCorrelation\n", info->mode);
					exit(EXIT_FAILURE);
					break;
				}
			}
		}
		}
		}
		switch(info->mode&0x0000000f) {
			case 0: { /* Normal Correlation */
				eval = sqrt(eval/sum); 
				break;
			}
			case 1: { /* Phase Correlation */
				eval =      eval/sum;
				break;
			}
			case 2: { /* Weighted Phase Correlation */
				eval = sqrt(eval/sum);
				break;
			}
		}
		DEBUGPRINT4("(%g %g %g): %g\n", x, y, z, eval);
		mrcPixelDataSet(out, ix, iy, iz, eval, mrcPixelRePart);
	}
	}
	}
	info->corInfo.mode = meanOfAll;
	lmrcImageInformation(&(info->corInfo), out);
	info->corInfo.minCoord.x = sx + info->corInfo.minCoord.x*dx;   
	info->corInfo.minCoord.y = sy + info->corInfo.minCoord.y*dy;   
	info->corInfo.minCoord.z = sz + info->corInfo.minCoord.z*dz;   

	shiftFFT = inFFT;
	mrcInit(&shiftFFT, NULL);

	x = info->corInfo.minCoord.x;
	y = info->corInfo.minCoord.y;
	z = info->corInfo.minCoord.z;
	DEBUGPRINT4("min (%g %g %g) %g", x, y, z, info->corInfo.min);
	for(X=0;                    X<=inFFT.HeaderN.x/2.0; X++) {
	for(Y=-inFFT.HeaderN.y/2.0; Y< inFFT.HeaderN.y/2.0; Y++) {
	for(Z=-inFFT.HeaderN.z/2.0; Z< inFFT.HeaderN.z/2.0; Z++) {
		deltaPhase = -2*M_PI*x*X/inFFT.HeaderN.x
					 -2*M_PI*y*Y/inFFT.HeaderN.y
					 -2*M_PI*z*Z/inFFT.HeaderN.z;
		mrcPixelDataGet(&inFFT,       X, Y, Z, &mag,    mrcPixelMag,   mrcPixelHowNearest);
		mrcPixelDataGet(&inFFT,       X, Y, Z, &phase,  mrcPixelPhase, mrcPixelHowNearest);
		mrcPixelDataSet(&shiftFFT, X, Y, Z, mag*cos(phase+deltaPhase),  mrcPixelRePart);
		mrcPixelDataSet(&shiftFFT, X, Y, Z, mag*sin(phase+deltaPhase),  mrcPixelImPart);
	}
	}
	}
	lmrcImageFFT(shift, &shiftFFT, 0);
	mrcImageFree(&shiftFFT, "in lmrcImageCorrelation");

	if(flagIn) {
		mrcImageFree(&inFFT, "in lmrcImageCorrelation");
	}
	if(flagRef) {
		mrcImageFree(&refFFT, "in lmrcImageCorrelation");
	}
}

double
lmrcImageLinearCorrelation(mrcImage* in1, mrcImage* in2, long mode)
{
	mrcImageParaTypeReal x, y, z;
	double data1, data2;
	double sum1, sum2;
	double avg1, avg2;
	double xy, xx, yy;
	double count;
	double r;

	sum1 = sum2 = count = 0.0;
	for(x=0; x<in1->HeaderN.x; x++) {
	for(y=0; y<in1->HeaderN.y; y++) {
	for(z=0; z<in1->HeaderN.z; z++) {
		mrcPixelDataGet(in1, x, y, z, &data1, mrcPixelRePart, mrcPixelHowNearest);
		sum1 += data1;
		mrcPixelDataGet(in2, x, y, z, &data2, mrcPixelRePart, mrcPixelHowNearest);
		sum2 += data2;
		count += 1;
	}
	}
	}

	avg1 = sum1/count;
	avg2 = sum2/count;

	xx = yy = xy = 0.0;
	for(x=0; x<in1->HeaderN.x; x++) {
	for(y=0; y<in1->HeaderN.y; y++) {
	for(z=0; z<in1->HeaderN.z; z++) {
		mrcPixelDataGet(in1, x, y, z, &data1, mrcPixelRePart, mrcPixelHowNearest);
		mrcPixelDataGet(in2, x, y, z, &data2, mrcPixelRePart, mrcPixelHowNearest);
		xx += SQR(data1-avg1);
		yy += SQR(data2-avg2);	
		xy += (data1-avg1)*(data2-avg2);
	}
	}
	}

	r = xy/sqrt(xx)/sqrt(yy);

	return r;
}

void 
lmrcImageCorrelation(mrcImage* out, mrcImage* in, mrcImage* ref, long mode)
{
    mrcImage  inFFT;
    mrcImage refFFT;
    mrcImage outFFT;
	long     flagIn, flagRef;
	mrcImageParaTypeReal X, Y, Z;
	double mag, phase;

	flagIn = flagRef = 0;
    if(!(   in->HeaderN.x==ref->HeaderN.x
          &&in->HeaderN.y==ref->HeaderN.y
          &&in->HeaderN.z==ref->HeaderN.z)) {
        fprintf(stderr, "Different size between in(%ld,%ld,%ld) and ref(%ld,%ld,%ld)\n"
                        ,in->HeaderN.x, ref->HeaderN.x
                        ,in->HeaderN.y, ref->HeaderN.y
                        ,in->HeaderN.z, ref->HeaderN.z);
        exit(EXIT_FAILURE);
    }

    if(IsImage(in, "", 0)&&IsImage(ref, "", 0)) {
        lmrcImageFFT(&inFFT,  in, 0);
        lmrcImageFFT(&refFFT, ref, 0);
		flagIn = flagRef = 1;
    } else if(IsFT(in, "", 0)&&IsImage(ref, "", 0)) {
        inFFT = *in;
        lmrcImageFFT(&refFFT, ref, 0);
		flagRef = 1;
    } else if(IsImage(in, "", 0)&&IsFT(ref, "", 0)) {
        lmrcImageFFT(&inFFT,  in, 0);
        refFFT = *ref;
		flagIn = 1;
    } else if(IsFT(in, "", 0)&&IsFT(ref, "", 0)) {
        inFFT = *in;
        refFFT = *ref;
    } else {
        fprintf(stderr, "Nut supported Mode (in %ld, ref %ld)\n", in->HeaderMode, ref->HeaderMode);
        exit(EXIT_FAILURE);
    }
	switch(mode&0x0000000f) {
		case 0: { /* Normal Correlation */
    		lmrcFFTFGconj(&outFFT, &inFFT, &refFFT);
			break;
		}
		case 1: { /* Phase Correlation */
    		lmrcFFTFGconj(&outFFT, &inFFT, &refFFT);
			Z = 0;
			for(X=0; X<=outFFT.HeaderN.x/2.0; X++) {
				for(Y=-outFFT.HeaderN.y/2.0; Y<outFFT.HeaderN.y/2.0; Y++) {
					mrcPixelDataSet(&outFFT, X, Y, Z, 1.0, mrcPixelMag);
				}
			}
			break;
		}
		case 2: { /* Weighted Phase Correlation */
    		lmrcFFTFGconj(&outFFT, &inFFT, &refFFT);
			Z = 0;
			for(X=0; X<=outFFT.HeaderN.x/2.0; X++) {
				for(Y=-outFFT.HeaderN.y/2.0; Y<outFFT.HeaderN.y/2.0; Y++) {
					mrcPixelDataGet(&outFFT, X, Y, Z, &mag, mrcPixelMag, mrcPixelHowNearest);
					mrcPixelDataSet(&outFFT, X, Y, Z, sqrt(mag), mrcPixelMag);
				}
			}
			break;
		}
		case 3: { /* Normalized Correlation */
			double F, G;
			double W;

			mrcPixelDataGet(&inFFT,  0, 0, 0, &F, mrcPixelMag, mrcPixelHowNearest);
			mrcPixelDataGet(&refFFT, 0, 0, 0, &G, mrcPixelMag, mrcPixelHowNearest);
			W = F*G;
    		lmrcFFTFGconj(&outFFT, &inFFT, &refFFT);
			Z = 0;
			for(X=0; X<=outFFT.HeaderN.x/2.0; X++) {
				for(Y=-outFFT.HeaderN.y/2.0; Y<outFFT.HeaderN.y/2.0; Y++) {
					mrcPixelDataGet(&outFFT, X, Y, Z, &mag, mrcPixelMag, mrcPixelHowNearest);
					mrcPixelDataSet(&outFFT, X, Y, Z, mag/W, mrcPixelMag);
				}
			}
			break;
		}
		default: {
			fprintf(stderr, "Not supported mode: %ld in lmrcImageCorrelation\n", mode);
			exit(EXIT_FAILURE);
			break;
		}
	}
	if(flagIn) {
		mrcImageFree(&inFFT, "in lmrcImageCorrelation");
	}
	if(flagRef) {
		mrcImageFree(&refFFT, "in lmrcImageCorrelation");
	}
	if(0 != (mode&0x00000010)) {
		mrcPixelDataSet(&outFFT, 0.0, 0.0, 0.0, 0.0, mrcPixelRePart);
		mrcPixelDataSet(&outFFT, 0.0, 0.0, 0.0, 0.0, mrcPixelImPart);
	}
    lmrcImageFFT(out, &outFFT, 0);
	mrcImageFree(&outFFT, "in lmrcImageCorrelation");
}

void
lmrcImageCorrelationModePrint(FILE* fpt)
{
	fprintf(fpt, "0: Normal Correlation FxG*\n");
	fprintf(fpt, "1: Phase  Correlation FxG*/     |FxG*|\n");
	fprintf(fpt, "2: Phase  Correlation FxG*/sqrt(|FxG*|)\n");
	fprintf(fpt, "3: Normalized Normal Correlation FxG*/|F||G|)\n");
	fprintf(fpt, "16: (0,0) = 0 \n");
	fprintf(fpt, "Refinement Correlation Map\n");
	fprintf(fpt, "0: sqrt(Sum |Fi-Gi|^2                 / Sum |Gi|^2)\n");
	fprintf(fpt, "1:      Sum |Gi||phaseFi - phaseGi|   / Sum |Gi|  \n");
	fprintf(fpt, "2: sqrt(Sum |Gi||phaseFi - phaseGi|^2 / Sum |Gi|  )\n");
}

