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

void
main(int argc, char* argv[]) 
{
	long status;
	mrcImageNoiseAddInfo info;
	mrcImage in;
	mrcImage noise;
	mrcImage out;
	mrcImageParaTypeReal x, y, z;
	mrcImageInformation  inInfo;
	long num;
	double data;

	init0(&info);
    argCheck(&info, argc, argv);
    init1(&info);
	
	mrcFileRead(&in, info.In, "in main", 0);
	lmrcImageInformation(&inInfo, &in);
	noise.Header = in.Header;
	mrcInit(&noise, NULL);
		
	for(x=0; x<in.HeaderN.x; x++) {
		for(y=0; y<in.HeaderN.y; y++) {
		  for(z=0; z<in.HeaderN.z; y++) {
			data = info.SD*randomNormalGet(2);
			mrcPixelDataSet(&noise, x, y, z, data, mrcPixelRePart); 
		  }
		}
	}
	num = 0;
	lmrcImageAdd(&out, &in, &num);
	lmrcImageAdd(&out, &noise, &num);
	mrcStatDataSet(&out, 0);
	mrcFileWrite(&out, info.Out, "in main", 0);
	exit(EXIT_SUCCESS);
}

void
additionalUsage()
{
}

double lmrcImageStandardDeviationwithLH(mrcImage* in, double H, double L)
{
/* variables */
	long x,y,z,n;
	double sum2,sum,data;

/* begin */
	sum2 = 0.0;
	sum = 0.0;
	n=0;
	for(z=0 ;z < in->HeaderN.z ;z++) {
	  for(x=0 ;x < in->HeaderN.x ;x++) {
	    for(y=0 ;y < in->HeaderN.y ;y++){
	      mrcPixelDataGet(in ,x ,y ,z ,&data ,mrcPixelRePart ,mrcPixelHowNearest);
	      if (data > L && data < H){
		sum2 += data*data;
		sum += data;
		n++;
	      }
	    }
	  }
	}
	sum /= n;
	sum2 /= n;
	sum2 -= sum*sum;
	return(sum2);
}

