/*
# %M% %Y% %I%
# The latest update : %G% at %U%
#
#%Z% mrcImageDensityInfo ver %I%
#%Z% Created by 
#%Z%
#%Z% Usage : mrcImageDensityInfo
#%Z% Attention
#%Z%
*/
static char __sccs_id[] = "%Z%mrcImageDensityInfo ver%I%; Date:%D% %Z%";
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>                  
#define GLOBAL_DECLARATION
#include "../inc/config.h"

#define DEBUG
#include "genUtil.h"
#include "Memory.h"
#include "mrcImage.h"
#include "lmrcImageModeChange.h"
#include "Vector.h"

typedef struct lmrcImageDensityInfoInfo {
	floatVector percentageVolume;
	mrcImageParaTypeRealCoord centreOfGravity;
	mrcImageParaTypeRealCoord centreOfGravitySD;
	mrcImageParaTypeRealCoord max;
	mrcImageParaTypeRealCoord maxSD;
	double                    deltaVolume; 
	mrcImage in;
} lmrcImageDensityInfoInfo;

typedef enum lmrcImageDensityInfoMode {
	mrcImageDensityInfoPercentageVolume  =1,
	mrcImageDensityInfoCentreOfGravity   =2,
	mrcImageDensityInfoMax               =4,
	mrcImageDensityInfoCentreOfGravitySD =8,
	mrcImageDensityInfoMaxSD             =16
} lmrcImageDensityInfoMode;

extern void lmrcImageDensityInfo(lmrcImageDensityInfoInfo* linfo, mrcImage* in, long mode);
extern void lmrcImageDensityInfoWrite(FILE* fpt, lmrcImageDensityInfoInfo* linfo, long mode);
extern floatVector lmrcImageDensityInfoPercentageVolume(mrcImage* in, double delta, long mode);
extern mrcImageParaTypeRealCoord lmrcImageDensityInfoCentreOfGravity(mrcImage* in, long mode);
extern mrcImageParaTypeRealCoord lmrcImageDensityInfoSD(mrcImage* in, mrcImageParaTypeRealCoord centre, long mode);
extern mrcImageParaTypeRealCoord lmrcImageDensityInfoMax(mrcImage* in, long mode);

int
main(int argc, char* argv[]) 
{
	mrcImageDensityInfoInfo  info;
	lmrcImageDensityInfoInfo linfo;
	mrcImage in;

	init0(&info);
    argCheck(&info, argc, argv);
    init1(&info);

	DEBUGPRINT("Program Start\n");
	mrcFileRead(&in, info.In, "main", 0);

	linfo.deltaVolume = info.Contour;
	lmrcImageDensityInfo(&linfo, &in, info.mode);

	linfo.in = in;
	lmrcImageDensityInfoWrite(info.fptOut, &linfo, info.mode);
	exit(EXIT_SUCCESS);
}

void
additionalUsage()
{
	fprintf(stderr, "----- Attention -----\n");
	fprintf(stderr, "In this program, positive values are regarded as density. \n");
	fprintf(stderr, "If a pixel has a negative value, the value of the pixel is regarded as zero.");
	fprintf(stderr, "----- Additional Usage -----\n");
	fprintf(stderr, "Mode: \n");
	fprintf(stderr, "%d: PercentageVolume\n", mrcImageDensityInfoPercentageVolume);
	fprintf(stderr, "%d: CentreOfGravity \n", mrcImageDensityInfoCentreOfGravity);
	fprintf(stderr, "%d: Max             \n", mrcImageDensityInfoMax);
}

void
lmrcImageDensityInfo(lmrcImageDensityInfoInfo* linfo, mrcImage* in, long mode)
{
	if(mode&mrcImageDensityInfoPercentageVolume) {
		linfo->percentageVolume = lmrcImageDensityInfoPercentageVolume(in, linfo->deltaVolume, mode);
	}
	if(mode&mrcImageDensityInfoCentreOfGravity) {
		linfo->centreOfGravity = lmrcImageDensityInfoCentreOfGravity(in, mode);
	}
	if(mode&mrcImageDensityInfoCentreOfGravitySD) {
		if(!(mode&mrcImageDensityInfoCentreOfGravity)) {
			linfo->centreOfGravity = lmrcImageDensityInfoCentreOfGravity(in, mode);
		}
		linfo->centreOfGravitySD = lmrcImageDensityInfoSD(in, linfo->centreOfGravity, mode);
	}
	if(mode&mrcImageDensityInfoMax) {
		linfo->max = lmrcImageDensityInfoMax(in, mode);
	}
	if(mode&mrcImageDensityInfoMaxSD) {
		if(!(mode&mrcImageDensityInfoMax)) {
			linfo->max = lmrcImageDensityInfoMax(in, mode);
		}
		linfo->maxSD = lmrcImageDensityInfoSD(in, linfo->max, mode);
	}
}

void
lmrcImageDensityInfoWrite(FILE* fpt, lmrcImageDensityInfoInfo* linfo, long mode)
{
	int i;
	if(mode&mrcImageDensityInfoPercentageVolume) {
		for(i=0; i<linfo->percentageVolume.size; i++) {
			fprintf(fpt, "%10.2f : %15.6g\n", (float)i*linfo->deltaVolume, linfo->percentageVolume.data[i]);	
		}
	}
	if(mode&mrcImageDensityInfoCentreOfGravity) {
		fprintf(fpt, "GC(NoStart):     %10.4f %10.4f %10.4f\n", linfo->centreOfGravity.x, linfo->centreOfGravity.y, linfo->centreOfGravity.z);	
		fprintf(fpt, "GC(Start):       %10.4f %10.4f %10.4f\n", linfo->centreOfGravity.x+linfo->in.HeaderStartN.x, 
												linfo->centreOfGravity.y+linfo->in.HeaderStartN.y,
												linfo->centreOfGravity.z+linfo->in.HeaderStartN.z);	
		fprintf(fpt, "GC(NoStart[A]):  %10.4f %10.4f %10.4f\n", 
												linfo->centreOfGravity.x*linfo->in.HeaderLength.x, 
												linfo->centreOfGravity.y*linfo->in.HeaderLength.y, 
												linfo->centreOfGravity.z*linfo->in.HeaderLength.z);	
		fprintf(fpt, "GC(Start[A]):    %10.4f %10.4f %10.4f\n", 
												(linfo->centreOfGravity.x+linfo->in.HeaderStartN.x)*linfo->in.HeaderLength.x, 
												(linfo->centreOfGravity.y+linfo->in.HeaderStartN.y)*linfo->in.HeaderLength.y,
												(linfo->centreOfGravity.z+linfo->in.HeaderStartN.z)*linfo->in.HeaderLength.z);	
	}
	if(mode&mrcImageDensityInfoCentreOfGravitySD) {
		fprintf(fpt, "GCSD([pixel]):   %10.4f %10.4f %10.4f\n", linfo->centreOfGravitySD.x, linfo->centreOfGravitySD.y, linfo->centreOfGravitySD.z);	
		fprintf(fpt, "GCSD([A]):       %10.4f %10.4f %10.4f\n", 
												linfo->centreOfGravitySD.x*linfo->in.HeaderLength.x, 
												linfo->centreOfGravitySD.y*linfo->in.HeaderLength.y, 
												linfo->centreOfGravitySD.z*linfo->in.HeaderLength.z);	
	}

	if(mode&mrcImageDensityInfoMax) {
		fprintf(fpt, "Max(NoStart):    %10.4f %10.4f %10.4f\n", linfo->max.x, linfo->max.y, linfo->max.z);	
		fprintf(fpt, "Max(Start):      %10.4f %10.4f %10.4f\n", linfo->max.x+linfo->in.HeaderStartN.x, 
												linfo->max.y+linfo->in.HeaderStartN.y,
												linfo->max.z+linfo->in.HeaderStartN.z);	
		fprintf(fpt, "Max(NoStart[A]): %10.4f %10.4f %10.4f\n", 
												linfo->max.x*linfo->in.HeaderLength.x, 
												linfo->max.y*linfo->in.HeaderLength.y, 
												linfo->max.z*linfo->in.HeaderLength.z);	
		fprintf(fpt, "Max(Start[A]):   %10.4f %10.4f %10.4f\n", 
												(linfo->max.x+linfo->in.HeaderStartN.x)*linfo->in.HeaderLength.x, 
												(linfo->max.y+linfo->in.HeaderStartN.y)*linfo->in.HeaderLength.y,
												(linfo->max.z+linfo->in.HeaderStartN.z)*linfo->in.HeaderLength.z);	
	}
	if(mode&mrcImageDensityInfoMaxSD) {
		fprintf(fpt, "MaxSD([pixel]):  %10.4f %10.4f %10.4f\n", linfo->maxSD.x, linfo->maxSD.y, linfo->maxSD.z);	
		fprintf(fpt, "MaxSD([A]):      %10.4f %10.4f %10.4f\n", 
												linfo->maxSD.x*linfo->in.HeaderLength.x, 
												linfo->maxSD.y*linfo->in.HeaderLength.y, 
												linfo->maxSD.z*linfo->in.HeaderLength.z);	
	}
}

floatVector
lmrcImageDensityInfoPercentageVolume(mrcImage* in, double delta, long mode)
{
	mrcImage dst;
	long npixel;
	long nlevel; 
	long index, i, j, jstart;
	floatVector v;
	float*      fdata;
	double sum, now, flevel;
	extern int __lmrcImageDensityInfoPercentageVolumeCompareFloat(const void*, const void*);

	npixel = in->HeaderN.x*in->HeaderN.y*in->HeaderN.z;
	nlevel = (int)(100.0/delta+0.5);

	floatVectorInit(&v, nlevel);

	switch(in->HeaderMode) {
		case mrcFloatImage: {
			dst = *in;
			break;
		}
		default: {
			lmrcImageModeChange(&dst, in, mrcFloatImage);
			break;
		}
	}
	switch(dst.HeaderMode) {
		case mrcFloatImage: {
			fdata = (float*)memoryAllocate(sizeof(float)*npixel, "in lmrcImageDensityInfoPercentageVolume");
			memcpy(fdata, dst.ImageFloatImage, npixel*sizeof(float));
			qsort(fdata, npixel, sizeof(float), __lmrcImageDensityInfoPercentageVolumeCompareFloat); 
			break;
		}
		default: {
			fprintf(stderr, "Not supported mode in lmrcImageDensityInfoPercentageVolume \n");
			exit(EXIT_FAILURE);
			break;
		}
	}

	sum = 0.0;
	for(i=0; i<npixel; i++) {
		if(fdata[i]<0) {
			fdata[i] = 0;
		}
		sum += fdata[i];	
	}

	now = 0;
	jstart = 0;
	for(i=0; i<v.size; i++) {
		flevel = sum/nlevel*i;
		for(j=jstart; j<npixel; j++) {
			now += fdata[j]; 
			if(flevel <= now) {
				v.data[i] = fdata[j];
				jstart = j+1;
				break;
			}
		}
	}
	return v;
}

mrcImageParaTypeRealCoord
lmrcImageDensityInfoCentreOfGravity(mrcImage* in, long mode)
{
	mrcImageParaTypeRealCoord c;
	float x, y, z;
	double data;
	double cx, cy, cz;
	double sum;

	cx = cy = cz = 0;
	sum = 0;
	for(x=0; x<in->HeaderN.x; x++) {
	for(y=0; y<in->HeaderN.y; y++) {
	for(z=0; z<in->HeaderN.z; z++) {
		mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);	
		if(0<data) {
			cx  += x*data;
			cy  += y*data;
			cz  += z*data;
			sum += data;
		}
	}
	}
	}
	c.x = cx/sum;
	c.y = cy/sum;
	c.z = cz/sum;

	return c;
}

mrcImageParaTypeRealCoord
lmrcImageDensityInfoMax(mrcImage* in, long mode)
{
	mrcImageParaTypeRealCoord p;
	float x, y, z;
	double data;
	double max;

	p.x = p.y = p.z = 0;
	max = 0;
	for(x=0; x<in->HeaderN.x; x++) {
	for(y=0; y<in->HeaderN.y; y++) {
	for(z=0; z<in->HeaderN.z; z++) {
		mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);	
		if(max<data) {
			p.x = x;	
			p.y = y;	
			p.z = z;	
			max = data;
		}
	}
	}
	}

	return p;
}

mrcImageParaTypeRealCoord
lmrcImageDensityInfoSD(mrcImage* in, mrcImageParaTypeRealCoord p0, long mode)
{
	mrcImageParaTypeRealCoord p;
	float x, y, z;
	double weight;
	double sumWeight;

	p.x = p.y = p.z = 0;
	sumWeight = weight = 0.0;
	for(x=0; x<in->HeaderN.x; x++) {
	for(y=0; y<in->HeaderN.y; y++) {
	for(z=0; z<in->HeaderN.z; z++) {
		mrcPixelDataGet(in, x, y, z, &weight, mrcPixelRePart, mrcPixelHowNearest);	
		p.x += (SQR(x-p0.x)*weight);	
		p.y += (SQR(y-p0.y)*weight);	
		p.z += (SQR(z-p0.z)*weight);	
		sumWeight += weight;
	}
	}
	}
	p.x = sqrt(p.x/sumWeight);
	p.y = sqrt(p.y/sumWeight);
	p.z = sqrt(p.z/sumWeight);
	return p;
}

int 
__lmrcImageDensityInfoPercentageVolumeCompareFloat(void* a, void* b)
{
	if(*(float*)a<*(float*)b) {
		return 1;
	} else if(*(float*)a==*(float*)b) {
		return  0;
	} else {
		return -1;
	}
}


