/*
# %M% %Y% %I%
# The latest update : %G% at %U%
#
#%Z% mrcImage3DPad ver %I%
#%Z% Created by 
#%Z%
#%Z% Usage : mrcImage3DPad
#%Z% Attention
#%Z%
*/
static char __sccs_id[] = "%Z%mrcImage3DPad 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 "mrcImage.h"

/* prototypes */
void lmrcImage3DPadfilament(mrcImage* out ,mrcImage* in ,float rmax ,int Nx ,int Ny ,double w ,int mode ,double v ,long flagv);
double lmrcImage3DPadfilamentCalcavr(mrcImage* in ,double rmax);


void
main(int argc, char* argv[]) 
{
	mrcImage3DPadInfo info;

/* variables */
	mrcImage in;
	mrcImage out;

/* input patameters ,file open */
	init0(&info);
	argCheck(&info, argc, argv);
	init1(&info);
	mrcFileRead(&in ,info.In ,"in main" ,0);
	
/* begin */
	DEBUGPRINT("Program Start\n");
	if (!info.flagr){
	  info.r=in.HeaderN.x*0.5;
	}
	if (!info.flagNx){
	  info.Nx=in.HeaderN.x;
	}
	if (!info.flagNy){
	  info.Ny=in.HeaderN.y;
	}
	if (!info.flagNz){
	  info.Nz=in.HeaderN.z;
	}

	lmrcImage3DPadfilament(&out ,&in ,info.r ,info.Nx ,info.Ny ,info.w ,info.mode ,info.v ,info.flagv);
	mrcFileWrite(&out ,info.Out ,"in main" ,0);
	exit(EXIT_SUCCESS);
}

void
additionalUsage()
{
  fprintf(stderr, "----- mode -----\n");
  fprintf(stderr, "0:Out of threshold is avrage of threshold plain.\n");
  fprintf(stderr, "1:Mulplying cos function.\n");
  fprintf(stderr, "2:Mulplying cos function and out of threshold is avrage of threshold plain.\n");

  /*
	fprintf(stderr, "----- mode2 -----\n");
	fprintf(stderr, "0:For filament.R is calculated from its axis.\n");
	fprintf(stderr, "1:For particle.R is calculated from its center.\n");
	*/	
}

void lmrcImage3DPadfilament(mrcImage* out ,mrcImage* in ,float rmax ,int Nx ,int Ny ,double w ,int mode ,double v ,long flagv)
{
/* variables */
  int x,y,z;
  double data,cx,cy,sx,sy,r2,av,c,dr2;

/* begin */
  DEBUGPRINT1("%f ",rmax);
  out->Header = in->Header;
  out->HeaderN.x = Nx;
  out->HeaderN.y = Ny;
  mrcInit(out,NULL);
  cx=(in->HeaderN.x-1)*0.5;
  cy=(in->HeaderN.y-1)*0.5;
  sx=(out->HeaderN.x-in->HeaderN.x)*0.5;
  sy=(out->HeaderN.y-in->HeaderN.y)*0.5;
  r2=SQR(rmax);
  switch(mode){
  case 0:
    if (!flagv){
      av=lmrcImage3DPadfilamentCalcavr(in ,rmax);
    } else {
      av=v;
    }
    DEBUGPRINT1("%f \n",av);
    for (z=0 ;z<out->HeaderN.z;z++){
      for (x=0;x<out->HeaderN.x;x++){
	for (y=0;y<out->HeaderN.y;y++){
	  mrcPixelDataSet(out ,x ,y ,z ,av, mrcPixelRePart);
	}
      }
    }
    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);
	  dr2=SQR(x-cx)+SQR(y-cy);
	  if (r2>=dr2 ){
	    mrcPixelDataSet(out ,x+sx ,y+sy ,z ,data, mrcPixelRePart);
	  }
	}
      }
    }
    break;
  case 1:
    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);
	  dr2=SQR(x-cx)+SQR(y-cy);
	  if (SQR(rmax-w)>dr2 ){
	    mrcPixelDataSet(out ,x+sx ,y+sy ,z ,data, mrcPixelRePart);
	  }else if (SQR(rmax+w)<dr2 ){
	    mrcPixelDataSet(out ,x+sx ,y+sy ,z ,0, mrcPixelRePart);
	  }else {
	    c= (cos( ( rmax-sqrt(dr2)) / w *PI/2 )+1 )/2;
	    mrcPixelDataSet(out ,x+sx ,y+sy ,z ,data*c, mrcPixelRePart);
	  }
	    
	}
      }
    }
    break;
  case 2:
    av=lmrcImage3DPadfilamentCalcavr(in ,rmax+w);
    fprintf(stdout,"av :%g\n",av); 
    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);
	  dr2=SQR(x-cx)+SQR(y-cy);
	  if (SQR(rmax-w)>dr2 ){
	    mrcPixelDataSet(out ,x+sx ,y+sy ,z ,data, mrcPixelRePart);
	  }else if (SQR(rmax+w)<dr2 ){
	    mrcPixelDataSet(out ,x+sx ,y+sy ,z ,av, mrcPixelRePart);
	  }else {
	    c= (cos( ( rmax-sqrt(dr2)) / w *PI/2 )+1 )/2;
	    mrcPixelDataSet(out ,x+sx ,y+sy ,z ,(data-av)*c+av , mrcPixelRePart);
	  }
	}
      }
    }  
  }
  mrcStatDataSet(out ,0);
}

double lmrcImage3DPadfilamentCalcavr(mrcImage* in ,double rmax)
{
/* variables */
  int x,y,z,n;
  double cx,cy,r2,sum,data,av;

/* begin */
  n=0;
  sum=0;
  cx=(in->HeaderN.x-1)*0.5;
  cy=(in->HeaderN.y-1)*0.5; 
  for (x=0;x<in->HeaderN.x;x++){
    for (y=0;y<in->HeaderN.y;y++){
      r2=SQR(x-cx)+SQR(y-cy);
      if (r2>=SQR(rmax) && r2<SQR(rmax+1) ){
	for (z=0 ;z<in->HeaderN.z;z++){
	  mrcPixelDataGet(in ,x ,y ,z ,&data, mrcPixelRePart, mrcPixelHowNearest);
	  sum+=data;
	  n++;
	}
      }
    }
  }
  av=sum/n;
  return(av);
}

