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

/* typedef */
typedef struct {
  mrcImage image;
  double sy;
  double ey;
  double ysize;
} lmrcpwzCorrelation;

/* prototype */
void lmrcpwzTemplateCorrelation(lmrcpwzCorrelation* out,mrcImage* in2 ,mrcImage* in,lmrcpwzCorrelation* template ,int sy, int ey ,int mode);
void lmrcRefImagepwzCorPeakFind(mrcImage* out ,double* maxx ,double* maxy ,double* maxdata ,mrcImage* in ,int N ,int dN ,double dp ,double dz ,int length ,int div ,int mode );
void lmrcMakepwzCortemplate(lmrcpwzCorrelation* out ,int N ,double dp ,double dz ,mrcpwzCorHeader* header ,mrcImage* in );
void lmrcImageMaxFind(mrcImageParaTypeRealCoord* maxpos,double* maxdata,mrcImage* in);
void  lmrcRefImagepwzRefPeakMake(double* refmaxx ,double* refmaxy ,mrcImage* template , mrcImageParaTypeRealCoord maxpos ,int N);
double llmrcRefImagepwzCorPeakFind(double* maxx ,double* maxy ,double* maxdata ,mrcImage* in ,double* refmaxx ,double* refmaxy ,int N ,int mode);
void lmrcRefImagepwzRefPeakMake2(double* refmaxx ,double* refmaxy ,int N ,double dp ,double dz ,mrcpwzCorHeader* header ,mrcImage* in, mrcImageParaTypeRealCoord maxpos );
void lmrcRefImagepwzCorPeakDisp(mrcImage* out ,double* maxx ,double* maxy ,double* maxdata ,int N ,mrcImage* in);
void lmrcRefImagepwzCorPeakAdd(mrcImage* out ,double* maxx ,double* maxy ,double* maxdata ,int N ,mrcImage* in);
void lmrcRefImagepwzCorPeakDisp2(mrcImage* out ,double** maxx ,double** maxy ,double** maxdata ,int Nx ,int Ny ,mrcImage* in);
void lmrcRefImagepwzCorPeak2to1(double* maxx,double* maxy,double *maxdata,double** refmaxx,double** refmaxy,double** refmaxdata,int dN,int N ,int tN);
void lmrcRefImagepwzCorPeak2to1MaxFind(double* maxx,double* maxy,double* maxdata,double** refmaxx,double** refmaxy,double** refmaxdata,int dN,int N,int tN ,int n,int m);
void lmrcRefImagepwzCorSmoothing(mrcImage* out ,mrcImage* in ,int N ,double dp ,double dz ,int mode);

void
main(int argc, char* argv[]) 
{
	long status;
	mrcRefImagepwzCorPeakFindInfo info;
/* variables */
	mrcImage in;
	mrcImage out;
	double* maxx;
	double* maxy;
	double* maxdata;

/* input patameters ,file open */
	init0(&info);
	argCheck(&info, argc, argv);
	init1(&info);

/* begin */
	mrcFileRead(&in ,info.In ,"in main" ,0);
	DEBUGPRINT("Read \n");
	if (info.mode2){
	  lmrcRefImagepwzCorSmoothing(&out,&in,info.N,info.dp,info.dz,info.mode);
	}else{
	  lmrcRefImagepwzCorPeakFind(&out ,maxx ,maxy ,maxdata ,&in ,info.N ,info.dN ,info.dp ,info.dz ,info.length ,info.div ,info.mode);
	}
	mrcFileWrite(&out ,info.Out ,"in main" ,0);
	exit(EXIT_SUCCESS);
}

void
additionalUsage()
{
fprintf(stdout,"------ mode ------\n");
fprintf(stdout,"0:Use full length template\n");
fprintf(stdout,"1:Use full length template and middle length template. Requires N and dN option\n");
fprintf(stdout,"2:Mode1+Average distance between templates and peaks. \n"); 
fprintf(stdout,"\n");
fprintf(stdout,"------ Mode 2 -----\n");
fprintf(stdout,"0:peak find\n");
fprintf(stdout,"1:Smoothing pwzcorrelation file using template.\n");
}

void lmrcRefImagepwzCorSmoothing(mrcImage* out ,mrcImage* in ,int N ,double dp ,double dz ,int mode)
{
/* variables */
  lmrcpwzCorrelation template,tmp;
  mrcImage in2;
  mrcpwzCorHeader header;
  mrcImageParaTypeRealCoord shift;
/* begin */
  lmrcpwzCorHeaderGet(&header ,in);
  lmrcMakepwzCortemplate(&template ,N ,dp ,dz ,&header ,in);
  lmrcpwzTemplateCorrelation(&tmp,&in2,in,&template,header.sy ,header.ey,1);
  out->Header=in->Header;
  mrcInit(out,NULL);
  lmrcpwzCorHeaderSet(out ,&header);
  lmrc3DZ1ImgPut(out,&tmp.image,1);
  mrcStatDataSet(out,0);

}

void lmrcRefImagepwzCorPeakFind(mrcImage* out ,double* maxx ,double* maxy ,double* maxdata ,mrcImage* in ,int N ,int dN ,double dp ,double dz ,int length ,int div ,int mode)
{
/* variables */
  int tN,n,i,mindisn,sn,en,sy,ey,m,j,idN;
  double lmaxdata,mindis,mindis3,syd,eyd,flag,tmpey,tmpsy;
  mrcpwzCorHeader header;
  lmrcpwzCorrelation template,template2,tmp,tmp2;
  mrcImageParaTypeRealCoord lmaxpos,shift;
  mrcImage in2,dum;
  double* refmaxx;
  double* refmaxy;
  double* maxdis;
  double** refmaxx2;
  double** refmaxy2;
  double** refmaxdata;

/* begin */
  lmrcpwzCorHeaderGet(&header ,in);
  tN=floor( (header.ey -header.sy) *in->HeaderLength.y /dz +0.1)-1 ;
  DEBUGPRINT5("tN %d %f %f %d %d \n",tN,header.ey,header.sy,div,N);
  if ( ( refmaxx = (double *)malloc(sizeof(double)*tN) ) == NULL ){
    fprintf(stderr ,"Enough memoly is not available.");
    exit(EXIT_FAILURE);
  }
  if ( ( refmaxy = (double *)malloc(sizeof(double)*tN) ) == NULL ){
    fprintf(stderr ,"Enough memoly is not available.");
    exit(EXIT_FAILURE);
  }
  if ( ( maxx = (double *)malloc(sizeof(double)*tN) ) == NULL ){
    fprintf(stderr ,"Enough memoly is not available.");
    exit(EXIT_FAILURE);
  }
  if ( ( maxy = (double *)malloc(sizeof(double)*tN) ) == NULL ){
    fprintf(stderr ,"Enough memoly is not available.");
    exit(EXIT_FAILURE);
  }
  if ( ( maxdata = (double *)malloc(sizeof(double)*tN) ) == NULL ){
    fprintf(stderr ,"Enough memoly is not available.");
    exit(EXIT_FAILURE);
  }
  
  if (mode==0){

    lmrcMakepwzCortemplate(&template ,tN ,dp ,dz ,&header ,in);
    lmrcpwzTemplateCorrelation(&tmp,&in2,in,&template,header.sy ,header.ey,1);
/*
      shift.x=shift.y=shift.z=0;
      out->Header=tmp.image.Header;
      mrcInit(out,NULL);
      lmrcImageCopy(out,&tmp.image,shift);
*/
    lmrcImageMaxFind(&lmaxpos,&lmaxdata,&tmp.image);
    DEBUGPRINT3("%f,%f,%f\n",lmaxpos.x,lmaxpos.y,lmaxpos.z);
    DEBUGPRINT("make peak\n");
    lmrcRefImagepwzRefPeakMake(refmaxx ,refmaxy ,&template.image ,lmaxpos ,tN);
    llmrcRefImagepwzCorPeakFind(maxx ,maxy ,maxdata ,&in2 ,refmaxx ,refmaxy ,tN ,0);
    lmrcRefImagepwzCorPeakDisp(out,maxx ,maxy ,maxdata,tN ,&in2); 
    mrcImageFree(&tmp.image,"");

  }else{
    DEBUGPRINT("malloc start ");
    if ( (refmaxx2 = (double **)malloc(sizeof(double*)*tN*div ) )== NULL ){
      fprintf(stderr ,"Enough memoly is not available.");
      exit(EXIT_FAILURE);
    }
    if ( (refmaxy2 = (double **)malloc(sizeof(double*)*tN*div ) )== NULL ){
      fprintf(stderr ,"Enough memoly is not available.");
      exit(EXIT_FAILURE);
    }
    if ( (refmaxdata = (double **)malloc(sizeof(double*)*tN*div ) )== NULL ){
      fprintf(stderr ,"Enough memoly is not available.");
      exit(EXIT_FAILURE);
    }
    for(n=0;n<tN*div;n++){
      if ( (refmaxx2[n] = (double *)malloc(sizeof(double)*N ) )== NULL ){
	fprintf(stderr ,"Enough memoly is not available.");
	exit(EXIT_FAILURE);
      }
      if ( (refmaxy2[n] = (double *)malloc(sizeof(double)*N ) )== NULL ){
	fprintf(stderr ,"Enough memoly is not available.");
	exit(EXIT_FAILURE);
      }
      if ( (refmaxdata[n] = (double *)malloc(sizeof(double)*N ) )== NULL ){
	fprintf(stderr ,"Enough memoly is not available.");
	exit(EXIT_FAILURE);
      }
    }
    DEBUGPRINT("malloc end\n");
    lmrcMakepwzCortemplate(&template2 ,N/div ,dp ,dz ,&header ,in);
    lmrcpwzTemplateCorrelation(&tmp,&in2,in,&template2,header.sy ,header.ey,1);
    lmrcMakepwzCortemplate(&template ,tN-N/div ,dp ,dz ,&header ,in);
    DEBUGPRINT5("template %g %g %d %d %d\n" ,tmp.sy ,tmp.ey , tN,N,div );
    lmrcpwzTemplateCorrelation(&tmp2,&dum,&tmp.image,&template,tmp.sy ,tmp.ey,2);
    DEBUGPRINT("mode2end\n");
    /*    mrcImageFree(&dum,""); */
    DEBUGPRINT("dumfree\n");
    lmrcImageMaxFind(&lmaxpos,&lmaxdata,&tmp2.image);     
    DEBUGPRINT3("%f,%f,%f\n",lmaxpos.x,lmaxpos.y,lmaxpos.z);
    lmrcRefImagepwzRefPeakMake(refmaxx ,refmaxy ,&template.image ,lmaxpos ,tN-N/div);
    llmrcRefImagepwzCorPeakFind(maxx ,maxy ,maxdata ,&tmp.image ,refmaxx ,refmaxy ,tN-N/div ,0);      
    n=0;
    i=0;
    if (mode==1){
      while(n<tN-N){
	lmaxpos.x=maxx[n];
	lmaxpos.y=maxy[n];
	lmaxpos.z=0;
	lmrcRefImagepwzRefPeakMake(refmaxx ,refmaxy,&template2.image ,lmaxpos ,N);
	llmrcRefImagepwzCorPeakFind(refmaxx2[i] ,refmaxy2[i] ,refmaxdata[i] ,&in2 ,refmaxx ,refmaxy ,N ,0);
	n += dN;
	i++;
      }
      lmrcRefImagepwzCorPeak2to1(maxx,maxy,maxdata,refmaxx2,refmaxy2,refmaxdata,dN,N,tN);
      lmrcRefImagepwzCorPeakDisp(out,maxx ,maxy ,maxdata,tN,&in2); 
      
    }else{
      
      if ( (maxdis = (double *)malloc(sizeof(double)*tN ) )== NULL ){
	fprintf(stderr ,"Enough memoly is not available.");
	exit(EXIT_FAILURE);
      }
      mindisn=0;
      mindis3=1000;
      if (tN-N-1<2*dN) dN=(tN-N-1)/2;
      DEBUGPRINT1("dN = %d ",dN);
      idN=N/div/2;
      while(n<tN-N){
	maxdis[i]=0;
	for(m=n;m<=n+N-2*idN;m+=idN){
	  lmaxpos.x=maxx[m];
	  lmaxpos.y=maxy[m];
	  lmaxpos.z=0;
	  lmrcRefImagepwzRefPeakMake(refmaxx ,refmaxy,&template2.image ,lmaxpos ,N/div);
	  maxdis[i]+=llmrcRefImagepwzCorPeakFind(refmaxx2[i] ,refmaxy2[i] ,refmaxdata[i] ,&in2 ,refmaxx ,refmaxy ,N/div ,1);
	  DEBUGPRINT2("%d %f  ",m,maxdis[i]);
	}
	if (i>1 && mindis3>maxdis[i]+maxdis[i-1]+maxdis[i-2]){
	  mindis3=maxdis[i]+maxdis[i-1]+maxdis[i-2];
	  mindisn=i;
	}
	DEBUGPRINT3("%d %f %f\n",n,maxdis[i],mindis3);
	n += dN;
	i++;
      }
      (maxdis[mindisn]>maxdis[mindisn-1])?(mindis=maxdis[mindisn-1],i=mindisn-1):(mindis=maxdis[mindisn],i=mindisn);
      (mindis>maxdis[mindisn-2])?(mindis=maxdis[mindisn-2],i=mindisn-2):(mindis=maxdis[mindisn]);
      mindisn=i;
      if (mindisn>=(tN-N)/dN-1){
	sn=(mindisn-1)*dN;
	en=tN-N-1;
      }else if (mindisn==0){
	sn=0;
	en=dN-1;
      }else {
	sn=(mindisn-1)*dN;
	en=(mindisn+1)*dN;
      }
      i=0;
      for(n=sn;n<=en;n++){
	maxdis[i]=0;
	j=0;tmpsy=tmpey=0;
	for(m=n;m<=n+N-2*idN;m+=idN){
	  lmaxpos.x=maxx[m];
	  lmaxpos.y=maxy[m];
	  lmaxpos.z=0;
	  DEBUGPRINT3("m:%d %f %f\n",m,maxx[m],maxy[m]);
	  lmrcRefImagepwzRefPeakMake(refmaxx ,refmaxy,&template2.image ,lmaxpos ,N/div);
	  if(j==0){
	    tmpsy=refmaxy[0];
	    tmpey=tmpsy+(N-1)*dz/in2.HeaderLength.y;
	  }
	  maxdis[i]+=llmrcRefImagepwzCorPeakFind(refmaxx2[i+j] ,refmaxy2[i+j] ,refmaxdata[i+j] ,&in2 ,refmaxx ,refmaxy ,N/div ,1);
	  j += en-sn+1;
	}	    
	if (mindis>=maxdis[i]){
	  mindis=maxdis[i];
	  mindisn=i;
	  eyd=tmpey;
	  syd=tmpsy;
	}
	DEBUGPRINT3("%d %f %f \n",n,maxdis[i],mindis);
	i++;
      }
      lmrcRefImagepwzCorPeakDisp(out,refmaxx2[mindisn] ,refmaxy2[mindisn] ,refmaxdata[mindisn],N/div,&in2); 	
      j=0;
      for(m=0;m<=N-2*idN;m+=idN){
	lmrcRefImagepwzCorPeakAdd(out,refmaxx2[mindisn+j],refmaxy2[mindisn+j] ,refmaxdata[mindisn+j],N/div,&in2);
	j += en-sn+1;
      }
      (header.inverse==1)?(flag-=0.5):(flag+=0.5);
      sy=floor( (eyd+syd-length)/2 +0.5+flag);
      if (sy<0) sy=0;
      ey=sy+length-1;
      if (ey>=in2.HeaderN.y){
	ey=in2.HeaderN.y-1;
	sy=ey-length+1;
	if (sy<0) sy=0;
      }
      fprintf(stdout,"%d %d %f\n",sy,ey,mindis);
    }
  }
}




void lmrcMakepwzCortemplate(lmrcpwzCorrelation* out ,int N ,double dp ,double dz ,mrcpwzCorHeader* header ,mrcImage* in)
{
/* variables */
  double z,p;
  int oz,op,n;
  
  union uniheader{
    mrcRefHeader ref;
    mrcpwzCorHeader cor;
  };   
  union uniheader tmpheader;

/* initialization */
  out->image.Header = in->Header;
  out->image.HeaderN.z = 1;
  out->image.HeaderN.x = in->HeaderN.x -1 ;
  mrcInit(&out->image ,NULL);
  tmpheader.cor = *header;
  p=0;
  z=dz/2;
  n=0;
  out->sy = 0;
  out->ey = floor(N*dz / in->HeaderLength.y+0.5);
  DEBUGPRINT1("%f ",out->ey);
  
/* begin */
  while (n < N){
    p=NormalDegree(p);
    op=lmrcRefGetzCordinate(p,header->iw ,header->ia,&(tmpheader.ref) ) - 1; 
    oz=floor(z / in->HeaderLength.y +0.5);
    mrcPixelDataSet(&out->image ,op ,oz ,0 ,1 ,mrcPixelRePart);
    n ++ ;
    z += dz ;
    p += dp ;
  }
  mrcStatDataSet(&out->image,0);
}

void lmrcpwzTemplateCorrelation(lmrcpwzCorrelation* out,mrcImage* in2 ,mrcImage* in,lmrcpwzCorrelation* template ,int sy, int ey ,int mode)
{
/* variables */
  mrcImage tmp;
  mrcpwzCorHeader header;
  int x,y;
/* begin */
  if (mode/2){
    lmrcImageCorrelation(&out->image,in,&template->image,0);
    for (x=0;x<in->HeaderN.x ;x++){
      for (y=0 ;y<sy-(int)template->sy ;y++){
	mrcPixelDataSet(&out->image,x,y,0,in->HeaderAMin,mrcPixelRePart);
      }
      for (y=ey-(int)template->ey+1 ;y<in->HeaderN.y ;y++){
	mrcPixelDataSet(&out->image,x,y,0,in->HeaderAMin,mrcPixelRePart);
      }
    }	
    DEBUGPRINT4("mode2:%d %d %g %g \n",sy-(int)template->sy,ey-(int)template->ey+1,sy-(int)template->sy ,ey-(template->ey-template->ey) );
    DEBUGPRINT4("%d %g %d %g \n",sy,template->sy,ey,template->ey );

  }else{
    lmrc3DZ1ImgGet(&tmp ,in ,1);
    lmrcImageRectangleGet(in2 ,&tmp ,0 ,0 ,in->HeaderN.x-2 ,in->HeaderN.y-1 );
    if (mode){
      lmrcImageCorrelation(&out->image ,in2 ,&template->image ,0);
      /*      DEBUGPRINT5("%d %d %f %f %d ",sy,ey,template->sy,template->ey,ey-(int)(template->ey-template->sy) ); */
      for (x=0;x<in->HeaderN.x ;x++){
	for (y=0 ;y<sy-(int)template->sy ;y++){
	  mrcPixelDataSet(&out->image,x,y,0,in2->HeaderAMin,mrcPixelRePart);
	}
	for (y=ey-(int)template->ey+1 ;y <in->HeaderN.y ;y++){
	  mrcPixelDataSet(&out->image,x,y,0,in2->HeaderAMin,mrcPixelRePart);
	}
      }	
    } else {
      lmrcpwzCorHeaderGet(&header ,in);
      lmrcImageCorrelation(&tmp ,in2 ,&template->image ,0);
      for (x=0;x<in->HeaderN.x ;x++){
	for (y=0 ;y<sy-(int)template->sy ;y++){
	  mrcPixelDataSet(&tmp,x,y,0,in2->HeaderAMin,mrcPixelRePart);
	}
	for (y=ey-(int)template->ey+1 ;y <in->HeaderN.y ;y++){
	  mrcPixelDataSet(&tmp,x,y,0,in2->HeaderAMin,mrcPixelRePart);
	}
      }	
      out->image.Header = in->Header;
      out->image.HeaderN.x = in->HeaderN.x-1;
      mrcInit(&out->image,NULL);
      lmrc3DZ1ImgPut(&out->image,&tmp ,1);
      lmrcpwzCorHeaderSet(&out->image,&header);
    }
  }
  out->ey=ey-(template->ey-template->sy);
  out->sy=sy-(int)template->sy;
  mrcStatDataSet(&out->image,0);
}

void lmrcImageMaxFind(mrcImageParaTypeRealCoord* maxpos,double* maxdata,mrcImage* in)
{
/* variables */
  int x,y,z;
  double data;

/* begin */
  mrcPixelDataGet(in ,0 ,0 ,0 ,maxdata ,mrcPixelRePart , mrcPixelHowNearest);
  maxpos->x=maxpos->y=maxpos->z=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 (data > *maxdata ){
	  *maxdata = data;
	  maxpos->x=x;
	  maxpos->y=y;
	  maxpos->z=z;
	}
      }
    }
  }
}
	

void  lmrcRefImagepwzRefPeakMake(double* refmaxx ,double* refmaxy ,mrcImage* template , mrcImageParaTypeRealCoord maxpos ,int N)
{
  int x,y,z,n;
  double data;
  mrcImage tmp;

  n=0;
  lmrcImageShift(&tmp ,template ,maxpos ,mrcPixelHowNearest);
  for (z=0 ;z<tmp.HeaderN.z ;z++){
    for (y=0 ;y<tmp.HeaderN.y ;y++){
      for (x=0 ;x<tmp.HeaderN.x ;x++){
	mrcPixelDataGet(&tmp ,x ,y ,z ,&data ,mrcPixelRePart , mrcPixelHowNearest);
	if (data != 0.0 && n<N){
	  refmaxx[n]=x;
	  refmaxy[n]=y;
	  n++;
	}
      }
    }
  }
  mrcImageFree(&tmp,"");
}
 

double llmrcRefImagepwzCorPeakFind(double* maxx ,double* maxy ,double* maxdata ,mrcImage* in ,double* refmaxx ,double* refmaxy ,int N ,int mode)
{
/* variables */
  int n;
  double avr;

/* begin */
  switch(mode){
  case 0:
    for (n=0 ;n<N ;n++){
      mrcImagePeakFind(&maxx[n] ,&maxy[n] ,&maxdata[n] ,in ,refmaxx[n] ,refmaxy[n] ,0);
      /*
      DEBUGPRINT4("%d %f %f %g\n",n,maxx[n],maxy[n],maxdata[n]);
      */
    }
    return(0);
    break;
  case 1:
    avr=0;
    for (n=0 ;n<N ;n++){
      mrcImagePeakFind(&maxx[n] ,&maxy[n] ,&maxdata[n] ,in ,refmaxx[n] ,refmaxy[n] ,0);
      /*
      DEBUGPRINT6("%d %f %f %f %f %g ",n,maxx[n],maxy[n],refmaxx[n],refmaxy[n],maxdata[n]);
      */
      if (fabs(refmaxx[n]-maxx[n]) >fabs(refmaxx[n]-maxx[n]-in->HeaderN.x)){
	avr += SQR(refmaxx[n]-maxx[n]-in->HeaderN.x);
      } else if (fabs(refmaxx[n]-maxx[n]) >fabs(refmaxx[n]-maxx[n]+in->HeaderN.x)){
	avr += SQR(refmaxx[n]-maxx[n]+in->HeaderN.x);
      }else{
	avr += SQR(refmaxx[n]-maxx[n]);
      }
      avr += SQR(refmaxy[n]-maxy[n]);
      /*
      DEBUGPRINT1("%g \n",avr);
      */
    }
    avr /= N;
    avr = sqrt(avr);
    return(avr);
  }
  return(-1);
}

void lmrcRefImagepwzCorPeakDisp(mrcImage* out ,double* maxx ,double* maxy ,double* maxdata ,int N ,mrcImage* in)
{
  int n;
  out->Header=in->Header;
  mrcInit(out,NULL);
  for (n=0 ;n<N ;n++){
    mrcPixelDataSet(out ,maxx[n] ,maxy[n] ,0 ,1 ,mrcPixelRePart);
    DEBUGPRINT2("%f %f   ",maxx[n] ,maxy[n]);
  }
  DEBUGPRINT("\n");
  mrcStatDataSet(out,0);
}

void lmrcRefImagepwzCorPeakAdd(mrcImage* out ,double* maxx ,double* maxy ,double* maxdata ,int N ,mrcImage* in)
{
  int n;
  for (n=0 ;n<N ;n++){
    mrcPixelDataSet(out ,maxx[n] ,maxy[n] ,0 ,1 ,mrcPixelRePart);
    DEBUGPRINT2("%f %f   ",maxx[n] ,maxy[n]);
  } 
  DEBUGPRINT("\n");
  mrcStatDataSet(out,0);
}

void lmrcRefImagepwzCorPeakDisp2(mrcImage* out ,double** maxx ,double** maxy ,double** maxdata ,int Nx ,int Ny ,mrcImage* in)
{
  int x,y;
  out->Header=in->Header;
  mrcInit(out,NULL);
  for (x=0 ;x<Nx ;x++){
    for (y=0 ;y<Ny ;y++){
      mrcPixelDataSet(out ,maxx[x][y] ,maxy[x][y] ,0 ,1 ,mrcPixelRePart);
    }
  }
  mrcStatDataSet(out,0);
}

void lmrcRefImagepwzRefPeakMake2(double* refmaxx ,double* refmaxy ,int N ,double dp ,double dz ,mrcpwzCorHeader* header ,mrcImage* in, mrcImageParaTypeRealCoord maxpos )
{
/* variables */
  union uniheader{
    mrcRefHeader ref;
    mrcpwzCorHeader cor;
  };   
  union uniheader tmpheader;
  double z,p,dummy;
  int oz,op,n;
  
/* initialization */
  tmpheader.cor = *header;
  lmrcRefGetdegreeCoordinate(&p ,&dummy ,&dummy ,floor(maxpos.x+1.5)+lmrcRefGetzCordinate(0,0 ,0,&(tmpheader.ref) ) ,&tmpheader.ref);
  z=dz/2 + maxpos.y*in->HeaderLength.y;
  n=0;
 
/* begin */
  while (n < N){
    p=NormalDegree(p);
    op=lmrcRefGetzCordinate(p,header->iw ,header->ia,&(tmpheader.ref) ) - 1; 
    oz=floor(z / in->HeaderLength.y +0.5);
    refmaxx[n]=op;
    refmaxy[n]=oz;
    n ++ ;
    z += dz ;
    p += dp ;
  }
}

void lmrcRefImagepwzCorPeak2to1(double* maxx,double* maxy,double *maxdata,double** refmaxx,double** refmaxy,double** refmaxdata,int dN,int N ,int tN)
{
  int n,m,on;
  
  m=0;
  on=0;
  for(m=0;m<N;m++){
    lmrcRefImagepwzCorPeak2to1MaxFind(&maxx[on],&maxy[on],&maxdata[on],refmaxx,refmaxy,refmaxdata,dN,N,tN,0,m);
    on++;
  }
  n=1;
  while(n<(tN-N)/dN){
    for(m=N-dN;m<N;m++){
      lmrcRefImagepwzCorPeak2to1MaxFind(&maxx[on],&maxy[on],&maxdata[on],refmaxx,refmaxy,refmaxdata,dN,N,(tN-N)/dN,n,m);    
      on++;
    }
    n++;
  }
}

void lmrcRefImagepwzCorPeak2to1MaxFind(double* maxx,double* maxy,double* maxdata,double** refmaxx,double** refmaxy,double** refmaxdata,int dN,int N,int tN ,int n,int m)
{
  int tmpn,tmpm,tmpmaxn,tmpmaxm;
  double tmpmax;

  tmpmaxn=n;
  tmpmaxm=m;
  tmpmax=refmaxdata[n][m];
  tmpn=n+1;
  tmpm=m-dN;
  DEBUGPRINT3("%d %d %g\n",n,m,refmaxdata[n][m]);
  while(tmpn<tN && tmpm>=0 ){
    DEBUGPRINT5("%d %d %g %g %g ",tmpn,tmpm,refmaxdata[tmpn][tmpm],refmaxx[tmpmaxn][tmpmaxm] ,refmaxy[tmpmaxn][tmpmaxm] );
    if (refmaxdata[tmpn][tmpm] >tmpmax){
      tmpmax=refmaxdata[tmpn][tmpm];
      tmpmaxn=tmpn;
      tmpmaxm=tmpm;
    }
    tmpn++;
    tmpm -= dN;
    DEBUGPRINT1("%g\n",tmpmax);
  }
  *maxx=refmaxx[tmpmaxn][tmpmaxm];
  *maxy=refmaxy[tmpmaxn][tmpmaxm];
  *maxdata=tmpmax;
}








