#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "../inc/config.h"
#undef DEBUG
#include "genUtil.h"
#include "mrcImage.h"

typedef struct infoDataList {
	mrcImage 			datXSection;	
	mrcImage 			datYSection;	
	mrcImage 			datXProj;	
	mrcImage 			datYProj;	
	double*             datHistgram;
	mrcImageInformation info;
} infoDataList;

extern void lmrcImageInfo(infoDataList* data, mrcImage* img, mrcImageInfoInfo* info);
extern void outputMrcImageInfo(FILE* fpt, infoDataList* data, mrcImage* img, mrcImageInfoInfo* info);

int
main(int argc, char* argv[])
{
    mrcImage img;
	mrcImageInfoInfo  info;
	infoDataList data;

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

#ifdef DEBUG
	mrcFileRead(&img, info.In, "from Main Routine", 1);
#else
    mrcFileRead(&img, info.In, "from Main Routine", 0);
#endif

    lmrcImageInfo(&data, &img, &info);
    outputMrcImageInfo(info.fptOut, &data, &img, &info);
	exit(EXIT_SUCCESS);
}

void
additionalUsage()
{
	fprintf(stderr, "----- Additional Usage -----\n");
	fprintf(stderr, "The -I option uses -rmax or -m options\n");
	fprintf(stderr, "-m 0: mean of all density\n");
	fprintf(stderr, "   1: mean of sphere surface density at rmax\n");
	fprintf(stderr, "   2: mean of cylinder surface density at rmax\n");
	fprintf(stderr, "   3: mean of edge pixels \n");
	fprintf(stderr, "   4: mean of Centre\n");
	fprintf(stderr, "   5: mean of Sparse \n");
	fprintf(stderr, "   6: mean of 2DEdge\n");
	fprintf(stderr, "   7: RMS of all pixels \n");
}

void
lmrcImageInfo(infoDataList* data, mrcImage* img, mrcImageInfoInfo* info)
{
    if(info->flagX) {
		lmrcImageXSection(&(data->datXSection), img, info->X);
	}
    if(info->flagY) {
		lmrcImageYSection(&(data->datYSection), img, info->Y);
    }
    if(info->xProjection) {
		lmrcImageXProjection(&(data->datXProj), img);
    }
    if(info->yProjection) {
		lmrcImageYProjection(&(data->datYProj), img);
    }
    if(info->flagH) {
		lmrcImageHistgram(&(data->datHistgram), info->H, img);
    }
    if(info->flagH2) {
		lmrcImageHistgram2(&(data->datHistgram), info->H2, img);
    }
    if(info->flagH3) {
		lmrcImageHistgram3(&(data->datHistgram), info->H3, img);
    }
    if(info->flagH4) {
		lmrcImageHistgram4(&(data->datHistgram), info->H4, img);
    }
    if(info->flagH5) {
		lmrcImageHistgram5(&(data->datHistgram), info->H5, info->Low, info->High, img);
    }
    if(info->flagH6) {
		lmrcImageHistgram5(&(data->datHistgram), info->H6, info->Low, info->High, img);
    }
    if(info->Info) {
		data->info.rmax     = info->rmax;
		data->info.flagrmax = info->flagrmax;
		data->info.mode     = info->mode;
		lmrcImageInformation(&(data->info), img);
	}
}

void
outputMrcImageInfo(FILE* fpt, infoDataList* dat, mrcImage* img, mrcImageInfoInfo* info)
{
    int i, ix, iy, j;
	double data;
	double dens;
		
    if(info->flagX) {
        for(iy=0; iy<img->HeaderN.y; iy++) {
            fprintf(fpt, "%5d %15.6g\n", iy, mrcPixelDataGet(&(dat->datXSection), 
															0.0,
															(mrcImageParaTypeReal)iy,
															0.0,
															&data,
															mrcPixelRePart,
															mrcPixelHowNearest));
        }
    }

    if(info->flagY) {
        for(ix=0; ix<img->HeaderN.x; ix++) {
            fprintf(fpt, "%5d %15.6g\n", ix, mrcPixelDataGet(&(dat->datYSection), 
															(mrcImageParaTypeReal)ix,
															0.0,
															0.0,
															&data,
															mrcPixelRePart,
															mrcPixelHowNearest));
        }
    }

    if(info->xProjection) {
        for(ix=0; ix<img->HeaderN.x; ix++) {
            fprintf(fpt, "%5d %15.6g\n", ix, mrcPixelDataGet(&(dat->datXProj), 
															(mrcImageParaTypeReal)ix,
															0.0,
															0.0,
															&data,
															mrcPixelRePart,
															mrcPixelHowNearest));
        }
    }

    if(info->yProjection) {
        for(iy=0; iy<img->HeaderN.y; iy++) {
            fprintf(fpt, "%5d %15.6g\n", iy, mrcPixelDataGet(&(dat->datYProj),
															0.0,
															(mrcImageParaTypeReal)iy,
															0.0,
															&data,
															mrcPixelRePart,
															mrcPixelHowNearest));
        }
    }
    if(info->flagH) {
        for(i=0; i<info->H; i++) {
			DEBUGPRINT2("MAX: %g, MIN: %g\n", img->HeaderAMax, img->HeaderAMin);
            dens = (img->HeaderAMax-img->HeaderAMin)/(info->H-1)*i+img->HeaderAMin, 
            fprintf(fpt, "%15.6g %15.6g |", 
                          dens, 
                          dat->datHistgram[i]);
			for(j=0; j<dat->datHistgram[i]*16; j++) {
				fprintf(fpt, " ");
			}
			fprintf(fpt, "*\n");
		}
	}
    if(info->flagH5) {
        for(i=0; i<info->H5; i++) {
			DEBUGPRINT2("MAX: %g, MIN: %g\n", img->HeaderAMax, img->HeaderAMin);
            dens = (info->High-info->Low)/(info->H5-1)*i+info->Low, 
            fprintf(fpt, "%15.6g %15.6g |", 
                          dens, 
                          dat->datHistgram[i]);
			for(j=0; j<dat->datHistgram[i]*16; j++) {
				fprintf(fpt, " ");
			}
			fprintf(fpt, "*\n");
		}
	}
    if(info->flagH2) {
		unsigned long n;
	
		n = (unsigned long)((img->HeaderAMax - img->HeaderAMin)/info->H2 + 1); 
        for(i=0; i<n; i++) {
            dens = img->HeaderAMin + i*info->H2;
            fprintf(fpt, "%15.6g %15.6g |", 
                          dens, 
                          dat->datHistgram[i]);
			
			for(j=0; j<dat->datHistgram[i]*16; j++) {
				fprintf(fpt, " ");
			}
			fprintf(fpt, "*\n");
        }
   	} 
    if(info->flagH6) {
		unsigned long n;
	
		n = (unsigned long)((info->High - info->Low)/info->H6 + 1); 
        for(i=0; i<n; i++) {
            dens = info->Low + i*info->H6;
            fprintf(fpt, "%15.6g %15.6g \n", 
                          dens, 
                          dat->datHistgram[i]);
			
        }
   	} 
    if(info->flagH3) {
        for(i=0; i<info->H; i++) {
			DEBUGPRINT2("MAX: %g, MIN: %g\n", img->HeaderAMax, img->HeaderAMin);
            dens = (img->HeaderAMax-img->HeaderAMin)/(info->H3-1)*i+img->HeaderAMin, 
            fprintf(fpt, "%15.6g %15.6g \n", 
                          dens, 
                          dat->datHistgram[i]);
		}
	}
    if(info->flagH4) {
		unsigned long n;
	
		n = (unsigned long)((img->HeaderAMax - img->HeaderAMin)/info->H4 + 1); 
        for(i=0; i<n; i++) {
            dens = img->HeaderAMin + i*info->H4;
            fprintf(fpt, "%15.6g %15.6g \n", 
                          dens, 
                          dat->datHistgram[i]);
			
        }
   	} 
    if(info->Info) {
		switch(info->mode){
			case meanOfAll:{
				DEBUGPRINT3("SizeCheck: int %d: mrcImageParaTypeIntegerCoord: %d = 3x%d \n", 
								sizeof(int), sizeof(mrcImageParaTypeIntegerCoord), sizeof(mrcImageParaTypeInteger));
				DEBUGPRINT3("SizeCheck: info:%d Coord%d Coord.x:%d\n", 
								sizeof(dat->info), sizeof(dat->info.minCoord), sizeof(dat->info.minCoord.x));
   	    		fprintf(fpt, "Min: %15.6g (%d, %d, %d)\n", dat->info.min, 
								(int)dat->info.minCoord.x, (int)dat->info.minCoord.y, (int)dat->info.minCoord.z);
   	     		fprintf(fpt, "Max: %15.6g (%d, %d, %d)\n", dat->info.max, 
								(int)dat->info.maxCoord.x, (int)dat->info.maxCoord.y, (int)dat->info.maxCoord.z);
   	     		fprintf(fpt, "Mean: %15.6g\n", dat->info.mean);
   	    		fprintf(fpt, "SD:  %15.6g\n", dat->info.sd);
   	     		fprintf(fpt, "SE:  %15.6g\n", dat->info.se);
   	     		fprintf(fpt, "Sum: %15.6g\n", dat->info.mean*img->HeaderN.x*img->HeaderN.y*img->HeaderN.z);
				break;
			}
			case meanOfCylinderSurface:{
   	     		fprintf(fpt, "Mean:%15.6g\n", dat->info.meanOfCylinderSurface);
   	     		fprintf(fpt, "SD:%15.6g\n", dat->info.sdOfCylinderSurface);
   	     		fprintf(fpt, "SE:%15.6g\n", dat->info.seOfCylinderSurface);
				break;
			}
			case meanOfSphereSurface:{
   	     		fprintf(fpt, "Mean:%15.6g\n", dat->info.meanOfSphereSurface);
   	     		fprintf(fpt, "SD:%15.6g\n", dat->info.sdOfSphereSurface);
   	     		fprintf(fpt, "SE:%15.6g\n", dat->info.seOfSphereSurface);
				break;
			}
			case meanOfEdge:{
   	     		fprintf(fpt, "Mean: %15.6g\n", dat->info.meanOfEdge);
   	     		fprintf(fpt, "SD:   %15.6g\n", dat->info.sdOfEdge);
   	     		fprintf(fpt, "SE:   %15.6g\n", dat->info.seOfEdge);
				break;
			}
			case RMSofAllPixels:{
   	     		fprintf(fpt, "RMS: %15.6g\n", dat->info.RMS);
				break;
			}
			default:{
				fprintf(stderr, "Not supported mode :%ld", info->mode);
				exit(EXIT_FAILURE);
			}

		}
   	} 
}
