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

#undef DEBUG
#include "genUtil.h"
#include "Random.h"
#include "mrcImage.h"

typedef struct lmrcImageAbnormalValueRemoveInfo {
	float cutValueOverSigma; 
	int   flagUnsigned;
	float UnsignedMax;
	float UnsignedHalf;
} lmrcImageAbnormalValueRemoveInfo;

typedef enum lmrcImageAbnormalValueRemoveMode {
	lmrcImageAbnormalValueModeDefault = 0,
	lmrcImageAbnormalValueModeNear    = 1,  
	lmrcImageAbnormalValueModeRandom  = 2,  
	lmrcImageAbnormalValueModeMean    = 3 
} lmrcImageAbnormalValueRemoveMode;

extern void
lmrcImageAbnormalValueRemove(mrcImage* out, 
							 mrcImage* in, 
							 lmrcImageAbnormalValueRemoveInfo* linfo,
							 lmrcImageAbnormalValueRemoveMode mode);

int
main(int argc, char* argv[]) 
{
	mrcImage in;
	mrcImage out;
	mrcImageAbnormalValueRemoveInfo info;
	lmrcImageAbnormalValueRemoveInfo linfo;
	
	init0(&info);
    argCheck(&info, argc, argv);
    init1(&info);

	linfo.cutValueOverSigma = info.cutValue;
	linfo.flagUnsigned = info.flagUnsignedMax;
	linfo.UnsignedMax = info.UnsignedMax;
	linfo.UnsignedHalf= info.UnsignedHalf;
	
	DEBUGPRINT("Program Start\n");
	mrcFileRead(&in, info.In, "in main", 0);

	out.Header = in.Header;
	out.HeaderMode = mrcFloatImage; 
	mrcInit(&out, NULL);
	
	lmrcImageAbnormalValueRemove(&out, &in, &linfo, info.mode);

	mrcFileWrite(&out, info.Out, "in main", 0);

	exit(EXIT_SUCCESS);
}

void
additionalUsage()
{
	fprintf(stderr, "----- Additional Usage -----\n");
	fprintf(stderr, "----- Mode: \n");
	fprintf(stderr, "    %d: Random Value \n", lmrcImageAbnormalValueModeDefault);
	fprintf(stderr, "    %d: Near Value \n", lmrcImageAbnormalValueModeNear);
	fprintf(stderr, "    %d: Random Value using SD and mean within Normal Area\n", lmrcImageAbnormalValueModeRandom);
	fprintf(stderr, "    %d: Mean Value using SD and mean within Normal Area\n", lmrcImageAbnormalValueModeMean);
}

void
lmrcImageAbnormalValueRemove(mrcImage* out, 
							 mrcImage* in, 
							 lmrcImageAbnormalValueRemoveInfo* linfo,
							 lmrcImageAbnormalValueRemoveMode mode)
{
	mrcImageInformation inInfo;
	mrcImage tmp;
	mrcImageParaTypeReal x, y, z;
	double data;

	if(linfo->flagUnsigned){
		tmp.Header = in->Header;	
		tmp.HeaderMode = mrcFloatImage;
		mrcInit(&tmp, NULL);
		for(z=0; z<in->HeaderN.z; z++) {
		for(y=0; y<in->HeaderN.y; y++) {
		for(x=0; x<in->HeaderN.x; x++) {
			mrcPixelDataGet(in,   x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
			if(linfo->UnsignedHalf <= data) {
				data -= linfo->UnsignedMax;	
			}
			mrcPixelDataSet(&tmp, x, y, z,  data, mrcPixelRePart);
		}
		}
		}
		in = &tmp;
	}
	inInfo.mode = meanOfAll;
	lmrcImageInformation(&inInfo, in);
	DEBUGPRINT2("%f +- %f\n", inInfo.mean, inInfo.sd);

	switch(mode) {
		case lmrcImageAbnormalValueModeDefault: {
			for(z=0; z<in->HeaderN.z; z++) {
			for(y=0; y<in->HeaderN.y; y++) {
			for(x=0; x<in->HeaderN.x; x++) {
				mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
				if(data<0) {
					fprintf(stderr, "data is negative: %lf\n", data);
				}
				if(data < inInfo.mean - linfo->cutValueOverSigma*inInfo.sd 
				  ||      inInfo.mean + linfo->cutValueOverSigma*inInfo.sd < data) {
					data = inInfo.mean + inInfo.sd*randomNormalGet(2);
					switch(out->HeaderMode) {
						case mrcCharImage:
						case mrcShortImage:
						case mrcLongImage: {
							if(data<0) {
								data = 0;
							}
							break;
						}
					}
				}
				mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
			}
			}
			}
			break;
		}
		case lmrcImageAbnormalValueModeRandom: {
			double sum, mean, SD;
			int n;
			sum = n = 0;
			for(z=0; z<in->HeaderN.z; z++) {
			for(y=0; y<in->HeaderN.y; y++) {
			for(x=0; x<in->HeaderN.x; x++) {
				mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
				if(inInfo.mean - linfo->cutValueOverSigma*inInfo.sd < data 
				 && data < inInfo.mean + linfo->cutValueOverSigma*inInfo.sd ) {
					sum += data;
					n++;
				}
			}
			}
			}
			if(0<n) {
				mean = sum/n;
			} else {
				fprintf(stderr, "No data within proper area\n");
			}
			for(z=0; z<in->HeaderN.z; z++) {
			for(y=0; y<in->HeaderN.y; y++) {
			for(x=0; x<in->HeaderN.x; x++) {
				mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
				if(inInfo.mean - linfo->cutValueOverSigma*inInfo.sd < data 
				 && data < inInfo.mean + linfo->cutValueOverSigma*inInfo.sd ) {
					sum += SQR(data - mean);
				}
			}
			}
			}
			SD = sqrt(sum/n);

			for(z=0; z<in->HeaderN.z; z++) {
			for(y=0; y<in->HeaderN.y; y++) {
			for(x=0; x<in->HeaderN.x; x++) {
				mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
				if(data<0) {
					fprintf(stderr, "data is negative: %lf\n", data);
				}
				if(data < mean - linfo->cutValueOverSigma*inInfo.sd
				  ||      mean + linfo->cutValueOverSigma*inInfo.sd < data) {
					data = mean + SD*randomNormalGet(2);
					switch(out->HeaderMode) {
						case mrcCharImage:
						case mrcShortImage:
						case mrcLongImage: {
							if(data<0) {
								data = 0;
							}
							break;
						}
					}
				}
				mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
			}
			}
			}
			break;
		}
		case lmrcImageAbnormalValueModeMean: {
			double sum, mean, SD;
			int n;
			sum = n = 0;
			for(z=0; z<in->HeaderN.z; z++) {
			for(y=0; y<in->HeaderN.y; y++) {
			for(x=0; x<in->HeaderN.x; x++) {
				mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
				if(inInfo.mean - linfo->cutValueOverSigma*inInfo.sd < data 
				 && data < inInfo.mean + linfo->cutValueOverSigma*inInfo.sd ) {
					sum += data;
					n++;
				}
			}
			}
			}
			if(0<n) {
				mean = sum/n;
			} else {
				fprintf(stderr, "No data within proper area\n");
			}

			for(z=0; z<in->HeaderN.z; z++) {
			for(y=0; y<in->HeaderN.y; y++) {
			for(x=0; x<in->HeaderN.x; x++) {
				mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
				if(data<0) {
					fprintf(stderr, "data is negative: %lf\n", data);
				}
				if(data < mean - linfo->cutValueOverSigma*inInfo.sd
				  ||      mean + linfo->cutValueOverSigma*inInfo.sd < data) {
					data = mean;
					switch(out->HeaderMode) {
						case mrcCharImage:
						case mrcShortImage:
						case mrcLongImage: {
							if(data<0) {
								data = 0;
							}
							break;
						}
					}
				}
				mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
			}
			}
			}
			break;
		}
		case lmrcImageAbnormalValueModeNear: {
			mrcImageParaTypeReal ox, oy, oz;
			mrcImageParaTypeReal dx, dy, dz;
			mrcImageParaTypeReal min, max;
			mrcImageParaTypeReal minHalf, maxHalf;
			double sum; 
			mrcImageParaTypeInteger count;
			mrcImageParaTypeInteger flag;

			min = inInfo.mean - linfo->cutValueOverSigma*inInfo.sd;
			minHalf = inInfo.mean - linfo->cutValueOverSigma*inInfo.sd/2.0;
			max = inInfo.mean + linfo->cutValueOverSigma*inInfo.sd;
			maxHalf = inInfo.mean + linfo->cutValueOverSigma*inInfo.sd/2.0;

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

				if(data < min  || max < data) {
					DEBUGPRINT6("data %f out of range %f <-> %f at (%f, %f, %f) \n", data, min, max, x, y, z);
					flag  = 0;
					dx = dy = dz = 1;
					while(flag<4) {	
					 	count = 0;	
						sum = 0;
						for(oz=z-dz; oz<=z+dz; oz++) {
						for(oy=y-dy; oy<=y+dy; oy++) {
						for(ox=x-dx; ox<=x+dx; ox++) {
							if(SQR((ox-x)/dx)+SQR((oy-y)/dy)+SQR((oz-z)/dz)<=1) {
						    	mrcPixelDataGet(in, ox, oy, oz, &data, mrcPixelRePart, mrcPixelHowNearest);
						    	if( min < data && data < max) {
									sum += data;
							    	count++;
						  		}
							}
						}
						}
						}
						if(4<=count) {
							flag ++;
							data = sum/count;  
							DEBUGPRINT5("data %f count %d at (%f, %f, %f) \n", data, count, x, y, z);
							if(data<minHalf || maxHalf < data ) {
								DEBUGPRINT("wider range +1 because of wide abnormal value\n");
								dx++;
								dy++;
								dz++;
							} else {
								flag = 4;	
							}

						} else {
							DEBUGPRINT("wider range +1\n");
							flag=0;
							dx++;
							dy++;
							dz++;
						}
					}

				}
				mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
			}
			}
			}
			break;
		}
		default: {
			fprintf(stderr, "Not supported mode: %d\n", mode);
			exit(EXIT_FAILURE);
			break;
		}
	}
	return;	
}


