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

int    w1,w2,w10,w20,mesh_x,mesh_y;                  
double *temp,*filter_x,*filter_y,**mat0,*vec0,sumvec[1000],sigma,sat,pixmin;                               
main(int argc,char *argv[])
{
  FILE   *infile,*fileout,*normf,*fp;
  char   header[144][80],outname[64],outnorm[64];
  int    i,j,k,n,buff,degre1,degre2,degre3,nc,*indx,nx,ny,cx,cy,nhead,
         deg_fond,nload_x,nload_y,iload_x,iload_y,offset,kernel_mesh;
  float  *px,*buff_imc,dummy;
  double *im,*imc,*imref,*imref2,**vectors,**mat,*vec,d,norm_kernel,sigm;

  if(argc<3){
    printf("\n\t\t******* alard2 98,10 ********");
    printf("\n\n\tUsage: alard mask.fits reference.fits");
    printf(" {list of images to subtract}\n\n");
    exit();
  }

/******************************************************************************/
/*************                List of Paramaters                ***************/
/******************************************************************************/

//  w10=w20=128;
//  w1 =w2 =128;
/**  w10=128; w20=128;  The dimensions of the image (nx,ny) ****/
/**  w1=128; w2=128;    The size of the sub-image where we  ****/
                   /** will find the constant kernel solution **/
//kernel_mesh=11;  /******** Half size of the Kernel ***********/
//sigma=1.2;       /* Size of the first gaussian core of the Kernel expansion */
//degre1=6;        /***********************************************************/
//degre2=4;   /* degrees of polynomial expension associated with each gaussian*/
//degre3=3;        /***********************************************************/
//deg_fond=0;      /****** Degree of polynomial background fitting ************/
//sat=29000.0;     /********* Maximum pixel value to fitted *******************/
//pixmin=4.0;      /********* Minimum pixel value to fitted *******************/
  
  fp=fopen("alard.par","r");
  if(fp==0){
    printf("\n\talard.par not found!\n");
    printf("\n\tcp /qso5/bin/alard.par .\n");
    printf("\n\tAnd vi it\n\n\n");
    exit();
  }
  fscanf(fp,"%d %d",&w10,&w20); while(fgetc(fp)!=10);
  fscanf(fp,"%d %d",&w1,&w2);   while(fgetc(fp)!=10);
  fscanf(fp,"%d",&kernel_mesh); while(fgetc(fp)!=10);
  fscanf(fp,"%lf",&sigma);      while(fgetc(fp)!=10);
  fscanf(fp,"%d %d %d",&degre1,&degre2,&degre3); while(fgetc(fp)!=10);
  fscanf(fp,"%d",&deg_fond);    while(fgetc(fp)!=10);
  fscanf(fp,"%lf",&sat);        while(fgetc(fp)!=10);
  fscanf(fp,"%lf",&pixmin);     
  fclose(fp);

  printf("size: x,y    %d %d\n",w10,w20);
  printf("sub_blocks   %d %d\n",w1,w2);
  printf("half_Kernel  %d\n",kernel_mesh);
  printf("sigma:       %lf\n",sigma);
  printf("polynomial   %d %d %d\n",degre1,degre2,degre3);
  printf("poly_backgr. %d\n",deg_fond);
  printf("max value    %lf\n",sat);
  printf("min value    %lf\n",pixmin);

  dummy=2*kernel_mesh;
  w1=(w10-dummy)/w1+dummy;
  w2=(w20-dummy)/w2+dummy;

/****************************************************************************/

  mesh_x=mesh_y=2*kernel_mesh+1;
  mesh_x=mesh_y=2*kernel_mesh+1;

  nc=(degre1+1)*(degre1+2)/2;
  nc+=(degre2+1)*(degre2+2)/2;
  nc+=(degre3+1)*(degre3+2)/2;
  nc+=(deg_fond+1)*(deg_fond+2)/2;

  buff=w1*w2*sizeof(double);

  im=(double *)malloc(buff);
  imc=(double *)malloc(buff);
  imref=(double *)malloc(buff);
  imref2=(double *)malloc(buff);
  temp=(double *)malloc(buff);
  filter_x=(double *)malloc(mesh_x*sizeof(double));
  filter_y=(double *)malloc(mesh_y*sizeof(double));
  vectors=(double **)malloc(nc*sizeof(double *));
  for (i=0;i<nc;i++) vectors[i]=(double *)malloc(buff);
  mat=(double **)malloc((nc+1)*sizeof(double *));
  for (i=0;i<=nc;i++) mat[i]=(double *)malloc((nc+1)*sizeof(double));
  mat0=(double **)malloc((nc+1)*sizeof(double *));
  for (i=0;i<=nc;i++) mat0[i]=(double *)malloc((nc+1)*sizeof(double));
  vec=(double *)malloc((nc+1)*sizeof(double));
  vec0=(double *)malloc((nc+1)*sizeof(double));
  indx=(int *)malloc((nc+1)*sizeof(double));
  for(n=0;n<w1*w2;n++) imc[n]=0.0;
  buff_imc=(float *)malloc((w1-mesh_x/2)*4);

  for(i=3;i<argc;i++){
    sprintf(outname,"sub_%s",argv[i]);
    printf("Initializing: %s\n", outname); 
    blank_image(outname,argv[i]);
  }

  nload_x=w10/(w1-mesh_x+1); nload_y=w20/(w2-mesh_y+1);
                                        /*nload_x = 2; nload_y = 2;*/
  for(iload_y=0;iload_y<nload_y;iload_y++)
  for(iload_x=0;iload_x<nload_x;iload_x++) {
    cx=(iload_x*(w1-mesh_x+1));
    cy=(iload_y*(w2-mesh_y+1));

    for(i=1;i<argc;i++) {
      sprintf(outname,"sub_%s",argv[i]);
      sprintf(outnorm,"nor_%s",argv[i]);
      infile=fopen(argv[i],"r");

      if(i==1) {
        printf("**** cx,cy=%d %d\n",cx,cy);
        nhead=hfread(header,infile);
        read_image(imref2,infile,cx,cy);
      }
      if(i==2) {
        nhead=hfread(header,infile);
        read_image(imref,infile,cx,cy);
        make_vectors(imref,imref2,vectors,mat,degre1,degre2,degre3,deg_fond);
        sprintf(outname,"ref0");
        fileout=fopen(outname,"wb");
        fwrite(im,buff,1,fileout);
        fclose(fileout);
      }
      if(i>2)  {
        nhead=hfread(header,infile);
        read_image(im,infile,cx,cy);
        printf("pouet1: %lf\n", im[64+w1*64]);
        for(j=0;j<=nc;j++) vec[j]=0.0;
        for(nx=mesh_x/2;nx<w1-mesh_x/2;nx++)
        for(ny=mesh_y/2;ny<w2-mesh_y/2;ny++){
          n=nx+w1*ny;
          if(imref[n]<sat&&imref[n]>=pixmin&&imref2[n]<sat && imref2[n]>=pixmin)
          for(j=0;j<nc;j++) vec[j+1] +=  vectors[j][n]*im[n]/sqrt(imref[n]);
        }

        for(n=0;n<=nc;n++){
          vec0[n]=vec[n];
          for(j=0;j<=nc;j++) mat0[n][j]=mat[n][j];
        }
        clip0(im,imref,imref2,imc,mat,vec,vectors,indx,nc);
                                        /*model(im,imref,imc,vec,vectors,nc);*/
        sigm=1.0E22;
        clip(im,imref,imref2,imc,mat,vec,vectors,indx,nc,&sigm);
        clip(im,imref,imref2,imc,mat,vec,vectors,indx,nc,&sigm);
        clip(im,imref,imref2,imc,mat,vec,vectors,indx,nc,&sigm);
        clip(im,imref,imref2,imc,mat,vec,vectors,indx,nc,&sigm);
        n=(deg_fond+1)*(deg_fond+2)/2;
        for(j=1+n,norm_kernel=0.0;j<=nc;j++) norm_kernel += vec[j]*sumvec[j-1];
        norm_kernel = 1.0/norm_kernel;

        write_kernel(vec,degre1,degre2,degre3,deg_fond,argv[i]);

        for(n=0;n<=nc;n++){
          vec[n]=vec0[n];
          for(j=0;j<=nc;j++) mat[n][j]=mat0[n][j];
        }
                 /*for(nx=mesh_x/2;nx<w1-mesh_x/2;nx++)
                   for(ny=mesh_y/2;ny<w2-mesh_y/2;ny++)
                   imc[nx+ny*w1] *= norm_kernel;*/
        printf("norm_kernel: %lf %s\n", norm_kernel,outnorm);
        normf=fopen(outnorm,"w");
        fprintf(normf,"%lf\n",  norm_kernel);
        fclose(normf);

        if(!(fileout=fopen(outname,"r+")))
        {printf("Cannot create output File: %s\n", outname); exit(0);}

        offset=(cx+(cy+mesh_y/2)*w10+mesh_x/2)*sizeof(float);
	offset+=nhead*80;

        fseek(fileout,offset,SEEK_SET);
        offset=(w10-w1+mesh_x/2)*sizeof(float);

        for(n=mesh_y/2;n<w2;n++){
          for(j=0;j<w1-mesh_x/2;j++)buff_imc[j]=(float)imc[n*w1+mesh_x/2+j];
          swap4(buff_imc,(w1-mesh_x/2)*4);
          fwrite(buff_imc,(w1-mesh_x/2)*4,1,fileout);
          fseek(fileout,offset,SEEK_CUR);
        }
        fclose(fileout);
        printf("Done \n");
      }
    }
  }
  free(im);
  free(imc);
  free(filter_x);
  free(filter_y);
}

/*** my routines ***/

hfread(head,fp)
char head[][80];
FILE *fp;
{
  int i=0;
l10:
  if(i>=144){
    printf(" not a fits file!\n");
    exit(0);
  }
  fread(&head[i][0],36,80,fp);
  i+=36;
  if(indexpos(head,"END     ",i)==i)goto l10;
  return i;
} 

/***:::::::::::::: blank_image.c ::::::::::::::***/

blank_image(char *fileout, char *infile)
{
  FILE   *outfile,*fp;
  float  *line;
  int    i,buff;
  char   head[36][80];

  buff=w10*sizeof(float);
  line=(float *)malloc(buff);
  for(i=0;i<w10;i++) line[i]=0.0;
  outfile=fopen(fileout,"wb");
  fp=fopen(infile,"rb");
  i=0;
l10:
  i++;
  if(i>4){
    printf("%s not a fits file!\n",infile);
    exit(0);
  }
  fread(head,36,80,fp);
  fwrite(head,36,80,outfile);
  if(indexpos(head,"END     ",36)==36)goto l10;

  for(i=0;i<w20;i++) fwrite(line,buff,1,outfile);
  free(line);
  fclose(outfile);
}

/***:::::::::::::: clip.c ::::::::::::::***/

clip(double *im,double *imref,double *imref2,double *imc,double **mat,
          double *vec,double **vectors,int *indx,int nc,double *sigm)
{
  int   nx,ny,i,j,n,ncp;
  double  sigp,d,clipc,nclip,sig;

  sig=0.0;
  ncp=-1;
  for(nx=mesh_x/2;nx<w1-mesh_x/2;nx++)
  for(ny=mesh_y/2;ny<w2-mesh_y/2;ny++) {
    n=nx+ny*w1; 
    if(im[n] >=pixmin && im[n]<sat && imref2[n]<sat && imref2[n] >=pixmin){ 
      sigp = imc[n]*imc[n]/im[n];
/*printf("sigp: %lf %lf %lf %i %i\n", imc[n]/sqrt(im[n]),imc[n],im[n],nx,ny); */
      if(sigp<16.0*(*sigm));                  /* ??????????? */
      {sig += sigp; ++ncp;}                   /* ??????????? */
    }
  }

  *sigm = sig/(double)ncp;
  clipc=16.0*(*sigm);
  printf("sigm: %f\n",*sigm);

  for(i=0;i<nc+1;i++){ vec[i]=vec0[i];
                       for(j=0;j<nc+1;j++) mat[i][j]=mat0[i][j];}
  nclip=0;
  for(nx=mesh_x/2;nx<w1-mesh_x/2;nx++)
  for(ny=mesh_y/2;ny<w2-mesh_y/2;ny++){
    n=nx+ny*w1;
    if(im[n]>0.0) sigp=imc[n]*imc[n]/im[n]; else sigp=1.0E22;
    if((imref[n]<sat && imref[n]>=pixmin && imref2[n]<sat && imref2[n]>=pixmin)
    && (im[n]>sat || sigp>clipc || im[n]<pixmin)) {
      for(i=0;i<nc;i++){
        vec[i+1] -= vectors[i][n]*im[n]/sqrt(imref[n]);
        for(j=0;j<=i;j++) mat[i+1][j+1] -= vectors[i][n]*vectors[j][n];
      }
    }
    ++nclip;
  }
  for(i=1;i<=nc;i++) for(j=1;j<i;j++) mat[j][i] = mat[i][j];
  ludcmp(mat,nc,indx,&d);
  lubksb(mat,nc,indx,vec);
  model(im,imref,imc,vec,vectors,nc);
  return;
}

/***:::::::::::::: clip0.c ::::::::::::::***/

clip0(double *im,double *imref,double *imref2,double *imc,
           double **mat,double *vec,double **vectors,int *indx,int nc)
{
  int   nx,ny,i,j,n,nclip;
  double d;

  for(i=0;i<nc+1;i++){
    vec[i]=vec0[i];
    for(j=0;j<nc+1;j++) mat[i][j]=mat0[i][j];
  }
  nclip=0;
  for(nx=mesh_x/2;nx<w1-mesh_x/2;nx++)
  for(ny=mesh_y/2;ny<w2-mesh_y/2;ny++) {
    n=nx+ny*w1;
    if((imref[n]<sat && imref[n]>=pixmin && imref2[n]<sat && imref2[n]>=pixmin)
    && (im[n]>sat || im[n]<pixmin)){
      for(i=0;i<nc;i++){
        vec[i+1] -= vectors[i][n]*im[n]/sqrt(imref[n]);
        for(j=0;j<=i;j++) mat[i+1][j+1] -= vectors[i][n]*vectors[j][n];
      }
      ++nclip;
    }
  }
  for(i=1;i<=nc;i++) for(j=1;j<i;j++) mat[j][i] = mat[i][j];
  ludcmp(mat,nc,indx,&d);
  lubksb(mat,nc,indx,vec);
  model(im,imref,imc,vec,vectors,nc);
  return;
}

/***:::::::::::::: indexx.c ::::::::::::::***/

indexx(int n,double *arrin,int *indx)
{
  int l,j,ir,indxt,i;
  double q;

  for (j=1;j<=n;j++) indx[j]=j;
  l=(n >> 1) + 1;
  ir=n;
  for (;;) {
    if (l > 1) q=arrin[(indxt=indx[--l])];
    else {
      q=arrin[(indxt=indx[ir])];
      indx[ir]=indx[1];
      if (--ir == 1) {
        indx[1]=indxt;
        return;
      }
    }
    i=l;
    j=l << 1;
    while (j <= ir) {
      if (j < ir && arrin[indx[j]] < arrin[indx[j+1]]) j++;
      if (q < arrin[indx[j]]) {
        indx[i]=indx[j];
        j += (i=j);
      }
      else j=ir+1;
    }
    indx[i]=indxt;
  }
}

/***:::::::::::::: lubksb.c ::::::::::::::***/

lubksb(a,n,indx,b)
double  **a,b[];
int n,*indx;
{
  int i,ii=0,ip,j;
  double  sum;

  for (i=1;i<=n;i++) {
    ip=indx[i];
    sum=b[ip];
    b[ip]=b[i];
    if (ii) for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
    else if (sum) ii=i;
    b[i]=sum;
  }
  for (i=n;i>=1;i--) {
    sum=b[i];
    for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
    b[i]=sum/a[i][i];
  }
}

/***:::::::::::::: ludcmp.c ::::::::::::::***/

#define TINY 1.0e-20;

ludcmp(a,n,indx,d)
int n,*indx;
double  **a,*d;
{
  int    i,imax,j,k;
  double big,dum,sum,temp2;
  double *vv,*lvector();

  vv=lvector(1,n);
  *d=1.0;
  for (i=1;i<=n;i++) {
    big=0.0;
    for (j=1;j<=n;j++) if ((temp2=fabs(a[i][j])) > big) big=temp2;
    if (big == 0.0) lnrerror("Singular matrix in routine LUDCMP");
    vv[i]=1.0/big;
  }
  for (j=1;j<=n;j++) {
    for (i=1;i<j;i++) {
      sum=a[i][j];
      for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
    }
    big=0.0;
    for (i=j;i<=n;i++) {
      sum=a[i][j];
      for (k=1;k<j;k++) sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
      if ( (dum=vv[i]*fabs(sum)) >= big) {
        big=dum;
        imax=i;
      }
    }
    if (j != imax) {
      for (k=1;k<=n;k++) {
        dum=a[imax][k];
        a[imax][k]=a[j][k];
        a[j][k]=dum;
      }
      *d = -(*d);
      vv[imax]=vv[j];
    }
    indx[j]=imax;
    if (a[j][j] == 0.0) a[j][j]=TINY;
    if (j != n) {
      dum=1.0/(a[j][j]);
      for (i=j+1;i<=n;i++) a[i][j] *= dum;
    }
  }
  lfree_vector(vv,1,n);
}

lnrerror(error_text)
char error_text;
{
  fprintf(stderr," Run error....");
  fprintf(stderr,"%s\n",error_text);
  fprintf(stderr,"Goodbye ! \n");
  exit(1);
}

double  *lvector(nl,nh)
int nl,nh;
{
  double  *v;
  v=(double  *)malloc((size_t) (nh-nl+1)*sizeof(double ));
  if(!v) lnrerror("allocation failure in vector()");
  return v-nl;
}

lfree_vector(v,nl,nh)
double  *v; int nl,nh;
{
  free((char*) (v+nl));
  return;
}

/***:::::::::::::: make_kernel.c ::::::::::::::***/

write_kernel(double *vec,int degre1,int degre2,int degre3,int deg_fond,char *s1)
{
  FILE   *outfile;
  int    i,j,k,ix,iy,n,nx,ny,buff,jx;
  double sigma2,px,sgn_i,sgn_j,sum_x,sum_y,*imk;
  char   fileout[64];

  sprintf(fileout,"ker_%s",s1);
  outfile=fopen(fileout,"wb");

  sigma2=sigma*sigma;
  buff=mesh_x*mesh_y*sizeof(double);
  imk=(double *)malloc(buff);
  for(ix=0;ix<mesh_x;ix++) for(jx=0;jx<mesh_x;jx++)  imk[ix+mesh_x*jx]=0.0;
  k=(deg_fond+1)*(deg_fond+2)/2;

  sgn_i=-1.0;
  for(i=0;i<=degre1;i++) {
    sgn_i *= -1.0;
    for(j=0,sgn_j=-1.0;j<=degre1-i;j++) {
      sgn_j *= -1.0;
      for(ix=1;ix<=mesh_x/2;ix++){
        px=(double)ix;
        filter_x[ix+mesh_x/2]=exp((double)i*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_x/2;ix<0;ix++){
        px=(double)ix;  
        filter_x[ix+mesh_x/2]=sgn_i*exp((double)i*log(fabs(px))-px*px/sigma2);
      }
      filter_x[mesh_x/2]=0.0;
      if(i==0) filter_x[mesh_x/2]=1.0;
    
      for(ix=1;ix<=mesh_y/2;ix++){
        px=(double)ix;
        filter_y[ix+mesh_y/2]=exp((double)j*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_y/2;ix<0;ix++){
        px=(double)ix;  
        filter_y[ix+mesh_y/2]=sgn_j*exp((double)j*log(fabs(px))-px*px/sigma2);
      }
      filter_y[mesh_y/2]=0.0;
      if(j==0) filter_y[mesh_y/2]=1.0;
      for(ix=-mesh_x/2,sum_x=0.0;ix<=mesh_x/2;ix++)sum_x+=filter_x[ix+mesh_x/2];
      for(ix=-mesh_y/2,sum_y=0.0;ix<=mesh_y/2;ix++)sum_y+=filter_y[ix+mesh_y/2];
      for(ix=0;ix<mesh_x;ix++)for(jx=0;jx<mesh_y;jx++)
      imk[ix+mesh_x*jx] += vec[k+1]*filter_y[jx]*filter_x[ix];
      ++k;
    }
  }
 
  sigma2 *= 3.0;
  sgn_i=-1.0;
  for(i=0;i<=degre2;i++) {
    sgn_i *= -1.0;
    for(j=0,sgn_j=-1.0;j<=degre2-i;j++) {
      sgn_j *= -1.0;
      for(ix=1;ix<=mesh_x/2;ix++){
        px=(double)ix;
        filter_x[ix+mesh_x/2]=exp((double)i*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_x/2;ix<0;ix++){
        px=(double)ix;  
        filter_x[ix+mesh_x/2]=sgn_i*exp((double)i*log(fabs(px))-px*px/sigma2);
      }
      filter_x[mesh_x/2]=0.0;
      if(i==0) filter_x[mesh_x/2]=1.0;
    
      for(ix=1;ix<=mesh_y/2;ix++){
        px=(double)ix;
        filter_y[ix+mesh_y/2]=exp((double)j*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_y/2;ix<0;ix++){
        px=(double)ix;  
        filter_y[ix+mesh_y/2]=sgn_j*exp((double)j*log(fabs(px))-px*px/sigma2);
      }
      filter_y[mesh_y/2]=0.0;
      if(j==0) filter_y[mesh_y/2]=1.0;
      for(ix=-mesh_x/2,sum_x=0.0;ix<=mesh_x/2;ix++)sum_x+=filter_x[ix+mesh_x/2];
      for(ix=-mesh_y/2,sum_y=0.0;ix<=mesh_y/2;ix++)sum_y+=filter_y[ix+mesh_y/2];
      for(ix=0;ix<mesh_x;ix++) for(jx=0;jx<mesh_y;jx++)
      imk[ix+mesh_x*jx]+=vec[k+1]*filter_y[jx]*filter_x[ix];
      ++k;
    }
  }

  sigma2 *= 3.0;
  sgn_i=-1.0;
  for(i=0;i<=degre3;i++) {
    sgn_i *= -1.0;
    for(j=0,sgn_j=-1.0;j<=degre3-i;j++) {
      sgn_j *= -1.0;
      for(ix=1;ix<=mesh_x/2;ix++){
        px=(double)ix;
        filter_x[ix+mesh_x/2]=exp((double)i*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_x/2;ix<0;ix++){
        px=(double)ix;  
        filter_x[ix+mesh_x/2]=sgn_i*exp((double)i*log(fabs(px))-px*px/sigma2);
      }
      filter_x[mesh_x/2]=0.0;
      if(i==0) filter_x[mesh_x/2]=1.0;
      for(ix=1;ix<=mesh_y/2;ix++){
        px=(double)ix;
        filter_y[ix+mesh_y/2]=exp((double)j*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_y/2;ix<0;ix++){
        px=(double)ix;  
        filter_y[ix+mesh_y/2]=sgn_j*exp((double)j*log(fabs(px))-px*px/sigma2);
      }
      filter_y[mesh_y/2]=0.0;
      if(j==0) filter_y[mesh_y/2]=1.0;
      for(ix=-mesh_x/2,sum_x=0.0;ix<=mesh_x/2;ix++)sum_x+=filter_x[ix+mesh_x/2];
      for(ix=-mesh_y/2,sum_y=0.0;ix<=mesh_y/2;ix++)sum_y+=filter_y[ix+mesh_y/2];
      for(ix=0;ix<mesh_x;ix++) for(jx=0;jx<mesh_y;jx++)
      imk[ix+mesh_x*jx]+=vec[k+1]*filter_y[jx]*filter_x[ix];
      ++k;
    }
  }

  fwrite(imk,buff,1,outfile);
  free(imk);
  fclose(outfile);
  return;
}

/***:::::::::::::: make_vectors.c ::::::::::::::***/

make_vectors(double *im, double *imref2, double **vectors,
             double **mat,int degre1, int degre2, int degre3, int deg_fond)
{
  FILE    *outfile;
  int     i,j,k,ix,iy,n,nx,ny,nclip,split_x,split_y;
  double  sigma2,px,sgn_i,sgn_j,sum_x,sum_y;

  sigma2=sigma*sigma;
  k=0;
  for(i=0;i<=deg_fond;i++) for(j=0;j<=deg_fond-i;j++){
    for(nx=0;nx<w1;nx++) for(ny=0;ny<w2;ny++){
      if(i>0 && j>0) vectors[k][nx+ny*w1] = pow((double)nx,(double)i)*
      pow((double)ny,(double)j)/sqrt(fabs(im[nx+ny*w1])+1.0);
      if(i==0 && j>0) vectors[k][nx+ny*w1] = pow((double)ny,(double)j)/
                                             sqrt(fabs(im[nx+ny*w1])+1.0);
      if(j==0 && i>0) vectors[k][nx+ny*w1] = pow((double)nx,(double)i)/
                                             sqrt(fabs(im[nx+ny*w1])+1.0);
      if(i==0 && j==0)vectors[k][nx+ny*w1] = 1.0/sqrt(fabs(im[nx+ny*w1])+1.0);
    }
    ++k;
  }
  sumvec[0]=0.0;
  sgn_i=-1.0;
  for(i=0;i<=degre1;i++) {
    sgn_i *= -1.0;
    for(j=0,sgn_j=-1.0;j<=degre1-i;j++) {
      sgn_j *= -1.0;
      for(ix=1;ix<=mesh_x/2;ix++){
        px=(double)ix;
        filter_x[ix+mesh_x/2]=exp((double)i*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_x/2;ix<0;ix++){
        px=(double)ix;  
        filter_x[ix+mesh_x/2]=sgn_i*exp((double)i*log(fabs(px))-px*px/sigma2);
      }
      filter_x[mesh_x/2]=0.0;
      if(i==0) filter_x[mesh_x/2]=1.0;
      for(ix=1;ix<=mesh_y/2;ix++){
        px=(double)ix;
        filter_y[ix+mesh_y/2]=exp((double)j*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_y/2;ix<0;ix++){
        px=(double)ix;  
        filter_y[ix+mesh_y/2]=sgn_j*exp((double)j*log(fabs(px))-px*px/sigma2);
      }
      filter_y[mesh_y/2]=0.0;
      if(j==0) filter_y[mesh_y/2]=1.0;
      for(ix=-mesh_x/2,sum_x=0.0;ix<=mesh_x/2;ix++)sum_x+=filter_x[ix+mesh_x/2];
      for(ix=-mesh_y/2,sum_y=0.0;ix<=mesh_y/2;ix++)sum_y+=filter_y[ix+mesh_y/2];
      sumvec[k] = sum_x*sum_y;
      split_x=split_y=1;
      if((i/2)*2 == i) split_x=0;
      if((j/2)*2 == j) split_y=0;  
      xy_convolve(im,vectors[k],split_x,split_y);
      for(n=0;n<w1*w2;n++){if(im[n]>pixmin) vectors[k][n] /= sqrt(im[n]);} ++k;
    }
  }
  sigma2 *= 3.0;
  sgn_i=-1.0;
  for(i=0;i<=degre2;i++) {
    sgn_i *= -1.0;
    for(j=0,sgn_j=-1.0;j<=degre2-i;j++) {
      sgn_j *= -1.0;
      for(ix=1;ix<=mesh_x/2;ix++){
        px=(double)ix;
        filter_x[ix+mesh_x/2]=exp((double)i*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_x/2;ix<0;ix++){
        px=(double)ix;  
        filter_x[ix+mesh_x/2]=sgn_i*exp((double)i*log(fabs(px))-px*px/sigma2);
      }
      filter_x[mesh_x/2]=0.0;
      if(i==0) filter_x[mesh_x/2]=1.0;
      for(ix=1;ix<=mesh_y/2;ix++){
        px=(double)ix;
        filter_y[ix+mesh_y/2]=exp((double)j*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_y/2;ix<0;ix++){
        px=(double)ix;  
        filter_y[ix+mesh_y/2]=sgn_j*exp((double)j*log(fabs(px))-px*px/sigma2);
      }
      filter_y[mesh_y/2]=0.0;
      if(j==0) filter_y[mesh_y/2]=1.0;
      for(ix=-mesh_x/2,sum_x=0.0;ix<=mesh_x/2;ix++)sum_x+=filter_x[ix+mesh_x/2];
      for(ix=-mesh_y/2,sum_y=0.0;ix<=mesh_y/2;ix++)sum_y+=filter_y[ix+mesh_y/2];
      if((i/2)*2 == i && (j/2)*2==j) sumvec[k] = sum_x*sum_y;
   
      sumvec[k] = sum_x*sum_y;

      split_x=split_y=1;
      if((i/2)*2 == i) split_x=0;
      if((j/2)*2 == j) split_y=0;  
      xy_convolve(im,vectors[k],split_x,split_y);
      for(n=0;n<w1*w2;n++){if(im[n]>pixmin)  vectors[k][n] /= sqrt(im[n]);} ++k;
    }
  }
  sigma2 *= 3.0;
  sgn_i=-1.0;
  for(i=0;i<=degre3;i++) {
    sgn_i *= -1.0;
    for(j=0,sgn_j=-1.0;j<=degre3-i;j++) {
      sgn_j *= -1.0;
      for(ix=1;ix<=mesh_x/2;ix++){
        px=(double)ix;
        filter_x[ix+mesh_x/2]=exp((double)i*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_x/2;ix<0;ix++){
        px=(double)ix;  
        filter_x[ix+mesh_x/2]=sgn_i*exp((double)i*log(fabs(px))-px*px/sigma2);
      }
      filter_x[mesh_x/2]=0.0;
      if(i==0) filter_x[mesh_x/2]=1.0;
      for(ix=1;ix<=mesh_y/2;ix++){
        px=(double)ix;
        filter_y[ix+mesh_y/2]=exp((double)j*log(px)-px*px/sigma2);
      }
      for(ix=-mesh_y/2;ix<0;ix++){
        px=(double)ix;  
        filter_y[ix+mesh_y/2]=sgn_j*exp((double)j*log(fabs(px))-px*px/sigma2);
      }
      filter_y[mesh_y/2]=0.0;
      if(j==0) filter_y[mesh_y/2]=1.0;
      for(ix=-mesh_x/2,sum_x=0.0;ix<=mesh_x/2;ix++)sum_x+=filter_x[ix+mesh_x/2];
      for(ix=-mesh_y/2,sum_y=0.0;ix<=mesh_y/2;ix++)sum_y+=filter_y[ix+mesh_y/2];

      sumvec[k] = sum_x*sum_y;
      split_x=split_y=1;
      if((i/2)*2 == i) split_x=0;
      if((j/2)*2 == j) split_y=0;  
      xy_convolve(im,vectors[k],split_x,split_y);
      for(n=0;n<w1*w2;n++){if(im[n]>pixmin) vectors[k][n] /= sqrt(im[n]);} ++k;
    }
  }

  for(i=0;i<k+1;i++) for(j=0;j<k+1;j++) mat[i][j]=0.0;
  nclip=0;
  for(nx=mesh_x/2;nx<w1-mesh_x/2;nx++) for(ny=mesh_y/2;ny<w2-mesh_y/2;ny++) {
    n=nx+w1*ny;
    if(im[n]<sat && im[n]>=pixmin && imref2[n]<sat && imref2[n]>=pixmin) { 
      ++nclip;
      for(i=0;i<k;i++)for(j=0;j<=i;j++)
      mat[i+1][j+1]+=vectors[i][n]*vectors[j][n];
    }
  }
  for(i=1;i<=k;i++) for(j=1;j<i;j++) mat[j][i] = mat[i][j];
  return;
}

/***:::::::::::::: model.c ::::::::::::::***/

model(double *im,double *imref,double *imc,double *vec,double **vectors,int nc)
{
  int nx,ny,n,j;
     
  for(nx=mesh_x/2;nx<w1-mesh_x/2;nx++)
  for(ny=mesh_y/2;ny<w2-mesh_y/2;ny++){
    n=nx+w1*ny; imc[n]=0.0;
    for(j=2;j<=nc;j++) imc[n] += vec[j]*vectors[j-1][n];
    imc[n] = imc[n]*sqrt(imref[n])+vec[1]-im[n];
  }
  return;
}

/***:::::::::::::: read_image.c ::::::::::::::***/

read_image(double *im, FILE *infile, int cx, int cy)
{
  int  i,j,k,offset;
  char cp[4],cp2[4];
  float *px,px2;

  offset=(cx+cy*w10)*4;
  fseek(infile,offset,SEEK_CUR);
  offset=(w10-w1)*4;
  for(k=0;k<w2;k++) {
    for(j=0;j<w1;j++){
      fread(cp2,4,1,infile);
      cp[0]=cp2[3]; cp[1]=cp2[2]; cp[2]=cp2[1]; cp[3]=cp2[0];
      px=(float *)&cp[0]; im[j+k*w1]=(double)*px;
    }
    if(im[j+k*w1]<1.0) im[j+k*w1]=1.0;
    fseek(infile,offset,SEEK_CUR);
  }
  fclose(infile);
}

/***:::::::::::::: xy_convolve.c ::::::::::::::***/

xy_convolve(double *im, double *imc, int split_x, int split_y)
{
  int    i,j,ij,buff;
  double norm,tmp,norm01x,norm02x,norm01y,norm02y,norm1x,norm1y;

  norm01x=norm02x=norm01y=norm02y=0.0;
  for(j=-mesh_x/2;j<0;j++) norm01x += filter_x[mesh_x/2-j];
  for(j=0;j<=mesh_x/2;j++) norm02x += filter_x[mesh_x/2-j];
  for(j=-mesh_y/2;j<0;j++) norm01y += filter_y[mesh_x/2-j];
  for(j=0;j<=mesh_y/2;j++) norm02y += filter_y[mesh_x/2-j];

  norm1x=norm01x+norm02x; norm1y=norm01y+norm02y;

  for(i=0;i<w1*w2;i++) {temp[i]=0.0; imc[i]=0.0;}

  if(split_x) {
    for(i=0;i<w2;i++) {
      buff=w1*i;
      for(ij=mesh_x/2;ij<w1-mesh_x/2;ij++) {
        temp[ij+buff]=0.0; norm = 0.0;
        if(im[ij+buff] >= pixmin && im[ij+buff] <sat){
          for(j=-mesh_x/2;j<0;j++){
            if(im[ij+j+buff] >= pixmin && im[ij+j+buff] <sat){
              temp[ij+buff] += im[ij+j+buff]*filter_x[mesh_x/2-j];
              norm += filter_x[mesh_x/2-j];
            }
          } 
          if(norm != 0.0) temp[ij+buff] /= norm/norm01x;
          tmp=0.0; norm = 0.0;
          for(j=0;j<=mesh_x/2;j++){
            if(im[ij+j+buff] >= pixmin && im[ij+j+buff] <sat){
              tmp += im[ij+j+buff]*filter_x[mesh_x/2-j];
              norm += filter_x[mesh_x/2-j];
            }
          } 
          if(norm != 0.0) temp[ij+buff] += tmp*norm02x/norm;
        }
      }
    }
  }
  else {
    for(i=0;i<w2;i++) {
      buff=w1*i;
      for(ij=mesh_x/2;ij<w1-mesh_x/2;ij++) { 
        temp[ij+buff]=0.0; norm = 0.0;
        if(im[ij+buff] >= pixmin && im[ij+buff] <sat){
          for(j=-mesh_x/2;j<=mesh_x/2;j++){
            if(im[ij+j+buff] >= pixmin && im[ij+j+buff] <sat){
              temp[ij+buff] += im[ij+j+buff]*filter_x[mesh_x/2-j];
              norm += filter_x[mesh_x/2-j];
            }
          } 
          if(norm != 0.0)  temp[ij+buff] /= norm/norm1x;
        }
      }
    }
  }
  
  if(split_y) {
    for(i=0;i<w1;i++) {   
      for(ij=mesh_y/2;ij<w2-mesh_y/2;ij++) {
        buff=ij*w1;
        imc[buff+i]=0.0; norm=0.0;
        if(im[buff+i]>= pixmin && im[buff+i] <sat){
          for(j=-mesh_y/2;j<0;j++){
            if(im[buff+j*w1+i] >= pixmin && im[buff+j*w1+i] <sat){
              imc[buff+i] += temp[buff+j*w1+i]*filter_y[mesh_y/2-j];
              norm += filter_y[mesh_y/2-j];
            }
          } 
          if(norm != 0.0) imc[buff+i] /= norm/norm01y;
          norm=0.0; tmp=0.0;
          for(j=0;j<=mesh_y/2;j++){
            if(im[buff+j*w1+i] >= pixmin && im[buff+j*w1+i] <sat){
              tmp += temp[buff+j*w1+i]*filter_y[mesh_y/2-j];
              norm += filter_y[mesh_y/2-j];
            }
          } 
          if(norm != 0.0)  imc[buff+i] += tmp*norm02y/norm;
        }
      }
    }
  }
  else {
    for(i=0;i<w1;i++) {   
      for(ij=mesh_y/2;ij<w2-mesh_y/2;ij++) {
        buff=ij*w1;
        imc[buff+i]=0.0; norm=0.0;
        if(im[buff+i]>= pixmin && im[buff+i] <sat){
          for(j=-mesh_y/2;j<=mesh_y/2;j++){
            if(im[buff+j*w1+i] >= pixmin && im[buff+j*w1+i] <sat){
              imc[buff+i] += temp[buff+j*w1+i]*filter_y[mesh_y/2-j];
              norm += filter_y[mesh_y/2-j];
            }
          }
          if(norm != 0.0) imc[buff+i] /= norm/norm1y;
        }
      }
    }
  }
  return;
}