/* @(#)mosdefine.c	8.5 (ESO-IPG) 11/30/94 15:51:43 */
/* @(#)mosdefine.c      6.1.1.4  (LSW)  07/02/94  13:48:49           */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* .COPYRIGHT   (C) 1994 Landessternwarte Heidelberg                 */
/* .IDENT        mosdefine.c                                         */
/* .AUTHOR       Sabine Moehler - LSW Heidelberg                     */
/* .KEYWORDS                                                         */
/* .PURPOSE      Search object with median: (val - med(val)) > thres */ 
/*               Center object with Gaussian fitting                 */
/*               (adapted from mosslin.c)                            */
/* .VERSION      1.0  07-FEB-1994  Implementation                    */
/*-------------------------------------------------------------------*/

#include <math.h> 
#include <tbldef.h>
#include <midas_def.h>
#include <ok.h>

#include <proto_mutil.h>
#include <proto_nrutil.h>

#define max(A,B) ((A)  > (B) ? (A) : (B)) 
#define min(A,B) ((A)  < (B) ? (A) : (B)) 

#define MAXLINES 1000
#define MAXSLITS 100

#define GAUSSIAN 1
#define MARK -9999.9

int    inull, os_col, o1_col, o2_col, ni_col, s1_col, s2_col, ss_col;
int    Ntot,Npix[2],Method,NseqO,NseqS,upper[MAXSLITS],lower[MAXSLITS];
int    iav, kun, knul, Xbin,Iwin,Width,Tid;
int    min_dist, min_sky;
int    ot_col, obj_typ;
double Start[2],Step[2], *Xgaus, *Ygaus, *A;
float  thresh;
char   qualif[4], o_typ[MAXSLITS];

#ifdef __STDC__
void  center_obj( float [], float [], int [], int ); 
void  search_obj( float [], int [], int, int *);
void  fit_obj( float [], int [], float [], int, int, float[]);
void  def_sky( float [], float [], int , int []);
float own_median(int, float [], float *);
#else
void  center_obj(); 
void  search_obj();
void  fit_obj();
void  def_sky();
float own_median();
#endif

main()

{
   int   i, nax, tidslit, imno, select;
   int   col,nslits,acol,arow,nsort,param[2], null, slit[MAXSLITS];
   int   sl_col, up_col, lo_col;

   char  method[21],obj_image[60],slit_table[60];
   char  window_table[60],unit[64],inst[72],line[80],linem[80];

   float *buff,*image, rpar[1];


/* get into Midas and get parameters */

   SCSPRO("mosdefine");
   
   SCKGETC("IN_A",1,60,&iav,obj_image);
   SCKGETC("IN_B",1,60,&iav,slit_table);
   SCKGETC("OUT_A",1,60,&iav,window_table);
   SCKGETC("INPUTC",1,3,&iav,method);
   SCKGETC ("QUALIF",  1, 4, &iav, qualif);

   SCKRDI("MIN_DIST",1,1,&iav,&min_dist,&kun,&knul);
   SCKRDI("MIN_SKY",1,1,&iav,&min_sky,&kun,&knul);
   SCKRDI("INPUTI",1,2,&iav,param,&kun,&knul);
   SCKRDR("INPUTR",1,1,&iav,rpar,&kun,&knul);


/* map input file */
   
   strcpy (inst," ");
   strcpy (unit," ");
   SCIGET(obj_image,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
		2,&nax,Npix,Start,Step,inst,unit,(char **)&image,&imno);

   Method = GAUSSIAN;
   sprintf(linem,"centering method:    Gaussian");

/* read table MOS to get limits of slitlets  */

   TCTOPN(slit_table,F_I_MODE,&tidslit);
   TCIGET(tidslit,&col,&nslits,&nsort,&acol,&arow);
   TCLSER(tidslit, "slit", &sl_col);
   TCLSER(tidslit, "ystart", &lo_col);
   TCLSER(tidslit, "yend", &up_col);
   if (toupper (qualif[0]) == 'S')
      TCLSER(tidslit, "ray_typ", &obj_typ);

   for (i=1;i<=nslits;i++) 
     {
      TCSGET(tidslit, i, &select);
      if (select)
        {
         TCERDI(tidslit, i, sl_col, &slit[i-1], &null);
         if (null) slit[i-1] = inull;
         TCERDI(tidslit, i, lo_col, &lower[i-1], &null);
         if (null) lower[i-1] = inull;
         TCERDI(tidslit, i, up_col, &upper[i-1], &null);
         if (null) upper[i-1] = inull;
         if (toupper (qualif[0]) == 'S')
           TCERDC(tidslit, i, obj_typ, &o_typ[i-1], &null);
	}
     }
   TCTCLO(tidslit);

/* open output table */

   if (toupper (qualif[0]) == 'S')
     TCTINI(window_table,F_TRANS,F_O_MODE,8,MAXLINES,&Tid);
   if (toupper (qualif[0]) == 'M')
     TCTINI(window_table,F_TRANS,F_O_MODE,7,MAXLINES,&Tid);
   SCDWRD(Tid,"Pixel",Step,1,1,&kun);
   TCCINI(Tid,D_R4_FORMAT,1,"F6.0" , "None ","Obj_Slit",&os_col);
   TCCINI(Tid,D_R4_FORMAT,1,"F10.2", "Pixel","Obj_Strt",&o1_col);
   TCCINI(Tid,D_R4_FORMAT,1,"F10.2", "Pixel","Obj_End",&o2_col);
   TCCINI(Tid,D_R4_FORMAT,1,"E12.3", "Pixel","Net_Intens",&ni_col);
   if (toupper (qualif[0]) == 'S')
     TCCINI(Tid,D_C_FORMAT,1,"A1", "","Obj_Typ",&ot_col);
   TCCINI(Tid,D_R4_FORMAT,1,"F10.2", "Pixel","Sky_Strt",&s1_col);
   TCCINI(Tid,D_R4_FORMAT,1,"F10.2", "Pixel","Sky_End",&s2_col);
   TCCINI(Tid,D_R4_FORMAT,1,"F6.0" , "None ","Sky_Slit",&ss_col);

/* set parameters */

   thresh  =  rpar[0];
   Width  =   param[0];
   Iwin  =   (param[0]-1)/2;
   Xbin  =   param[1];

   SCTPUT("search object ");
   SCTPUT("------------\n");
   sprintf(line,"Input image:         %s ",obj_image);
   SCTPUT(line);
   sprintf(line,"Input table:         %s ",slit_table);
   SCTPUT(line);
   sprintf(line,"Output table:        %s\n ",window_table);
   SCTPUT(line);
   SCTPUT("input parameters: ");
   sprintf(line,"search window:                                  %i pixels",Width); 
   SCTPUT(line);
   sprintf(line,"detection threshold:                            %6.2f DN",thresh);
   SCTPUT(line);
   sprintf(line,"minimum distance between object limits and sky: %d pixels",min_dist);
   SCTPUT(line);
   sprintf(line,"minimum number of CCD rows for valid sky:       %d rows",min_sky);
   SCTPUT(line);
   SCTPUT(linem);
   sprintf(line,"\nmedian on:                                   %i scan lines",Xbin);
   SCTPUT(line);

   Ntot = Start[1]+Step[1]*Npix[1];

   buff = (float *) osmmget(Ntot * sizeof (float));

   Xgaus = dvector(1,Width);
   Ygaus = dvector(1,Width);	
   A     = dvector(1,4);

/* center objects and derive sky limits */


   center_obj(image, buff, slit, nslits);
   NseqO--;
   NseqS--;

   SCDWRI(Tid, "NOBJ", &NseqO, 1, 1, &kun);        
   SCDWRI(Tid, "NSKY", &NseqS, 1, 1, &kun);        

/* epilogue ... */

   TCSINI(Tid);    
   TCTCLO(Tid);

   osmmfree((char *)buff);
   free_dvector(Xgaus, 1, Width);
   free_dvector(Ygaus, 1, Width);
   free_dvector(A, 1, 3);

   SCSEPI();

}
/*--------------------------------------------------------------------------*/
#ifdef __STDC__
void center_obj(float image[], float rval[], int slit[], int nslits)
#else
void center_obj( image, rval, slit, nslits)
  float image[],rval[];
  int   slit[],nslits;
#endif

{
   int    nact, nslit, j, k, jj, objlist[MAXLINES];
   int    center, isum;
   float  *sval, sum, obj[4],sky[3],dum;
   double scan_pos,*tmp;
   char   text[80];

   NseqO = 1;
   NseqS = 1;
   SCKRDD("SCAN", 1, 1, &iav, &scan_pos, &kun, &knul); 
   center = (int) ((scan_pos-Start[0])/Step[0]); /* transform world coordinates
                                                 to pixels */

   sval = (float *) osmmget(Ntot * sizeof (float));
   tmp = (double *) osmmget(Ntot * sizeof (double));
   
   /* loop over all slitlets */
   for (nslit=0; nslit <= nslits-1; nslit++) 
     {
       if (slit[nslit] > inull)
        {
	  sprintf(text,"Now searching slit nr. %d", slit[nslit]);
	  SCTPUT(text);
	  obj[0] = slit[nslit];
	 
	 /* median on XBIN columns at chosen position */
	  for (j=0;j<Ntot;j++) rval[j] = 0.0;

	  for (k=lower[nslit]-1; k<upper[nslit]-1; k++) 
	    {
	      isum = 0;
	      for (j=center-Xbin/2-1; j<=center+Xbin/2-0.5; j++) 
		{           /* MIDAS starts arrays with 1, C with 0 */
		  jj = j+Npix[0]*(k-(int)(Start[1])+1);
		  /*                jj = Npix[0]*k+j;   */
       		  tmp[isum+1] = (double) (image[jj]);
		  isum++;
		}
	      dum = (float) select_pos ( (int)(isum+1)/2 , isum, tmp); 
	      rval[k] = dum;
	      sval[k] = rval[k];
	    }


	 /* now search for object ...  */
	 search_obj(rval,objlist,nslit,&nact);
	 /* ... fit object ...  */
	 fit_obj(rval,objlist,obj,nact,nslit, sval);	
	 /* ... and define sky */
	 def_sky(sval, sky, nslit, slit);
	}
     }
   osmmfree((char *) sval);
   osmmfree((char *) tmp);
}
/*---------------------------------------------------------------------------*/
#ifdef __STDC__
void search_obj( float rval[], int objlist[], int nslit, int *nact)
#else
void search_obj( rval, objlist, nslit, nact )
  float rval[];
  int   objlist[],nslit,*nact;
#endif

{
   int number = 0,j,i,imax;
   float med,diff,max_intens,Thresh;

   for (j=lower[nslit]-1+Iwin;j<upper[nslit]-1-Iwin;j++) 
     {
      
      diff = rval[j]-own_median(Width,&rval[j-Iwin],&med);                  
      if (med<0.) med = -med;
      Thresh = (thresh>0.) ? thresh : -thresh*med;
      if (diff > Thresh) 
        {
         imax = j;
         max_intens = rval[j];
         for (i=j-Iwin;i<=j+Iwin;i++) 
           if (rval[i] > max_intens )  
             {
              max_intens = rval[i];
              imax = i;
             }
           objlist[number] =  imax;
           number++;
        }
     }
   *nact = number;
}
/*---------------------------------------------------------------------------*/
#ifdef __STDC__
void fit_obj( float rval[], int objlist[], float obj[], int nact, int nslit,
              float sval[])
#else
void fit_obj( rval, objlist, obj, nact, nslit, sval )
  float rval[],obj[],sval[];
  int   objlist[],nact,nslit;
#endif
{
   int    j,s1,s2,k,jj,ii,index;
   double Aold[3], obj_lim, fraction;
   float  back, med; 

   for (jj=0; jj<=2; jj++) Aold[jj] = 0.0;
   for (jj=0;jj<nact;jj++) 
     {
      index = objlist[jj];

/* Determine local background (but ensure to take it within the slitlet) */

      if (index-Iwin-10 > lower[nslit] && index+Iwin+10 < upper[nslit])
        { 
         back = own_median(Width,&rval[index-Iwin-10],&med);
         back = (back + own_median(Width,&rval[index+Iwin+10],&med))/2;
        }
      else if (index-Iwin-10 <= lower[nslit]) 
          back = own_median(Width,&rval[index+Iwin+10],&med);   
      else if (index+Iwin+10 >= upper[nslit]-1) 
          back = own_median(Width,&rval[index-Iwin-10],&med);
      med = own_median(Width,&rval[index],&med);
      if (back > med) back = med;

/* GAUSSIAN fit of object profile */

      A[1] = rval[index];
/*      A[2] = Start[1] + Step[1]*objlist[jj];  */
      A[2] = objlist[jj]+1;
      A[3] = Step[1];
      for (k = 1 , j = index-Iwin ; j<= index+Iwin; j++, k++) 
        {
/*         Xgaus[k] = Start[1] + Step[1]*objlist[jj]+k-Iwin-1;    */
         Xgaus[k] = objlist[jj]+k-Iwin;
         Ygaus[k] = rval[j]-back;
        }	
      fit_gauss(Xgaus,Ygaus,Width,A);  

/* Avoid finding the same object several times  */
/*      if (Aold[0] != A[0] || Aold[1] != A[1] || Aold[2] != A[2]) */
      /*      if ( (int) Aold[0] !=  (int) A[1]    
	  || (int) (Aold[1] - A[2]) > 0 || Aold[2] !=  A[3] ) */
      /* new object if difference in x > 0.5 -- otherwise many double 
	 identifications  TS*/
      if ( (int) (Aold[1] - A[2] + 0.5 ) != 0)
        { 
/*         for (ii=0; ii<=2; ii++) Aold[ii] = A[ii]; */
	  for (ii=0; ii<=2; ii++) Aold[ii] = A[ii+1];

/* object limits are defined at position where I(x) = fraction * I(x0) */

	  SCKRDD("INT_LIM",1,1,&iav,&fraction,&kun,&knul);
	  obj_lim = sqrt(-2*log(fraction));
	  obj[1] = ((int) (A[2]-obj_lim*A[3]));               /* lower limit */
	  if (obj[1] < lower[nslit]+1) obj[1] = lower[nslit]+1;
	  obj[2] = ((int) (A[2]+obj_lim*A[3]+0.5));           /* upper limit */
	  if (obj[2] > upper[nslit]-1) obj[2] = upper[nslit]-1;
	  obj[3] = A[1];                                 /* net peak intensity */

/* mark object's position for later sky determination (ensure that to stay 
   within slitlet) */

/*         s1 = obj[1]-3;
         s2 = obj[2]+2; */
	  s1 = obj[1]-min_dist;
	  s2 = obj[2]+min_dist-1; 
	  for (j = s1; j<= s2; j++) 
	    {
	      if (j < lower[nslit]) s1++;
	      if (j > upper[nslit]) s2--;
	    } 
	  for (j = s1; j<= s2; j++) sval[j] = MARK;
	  TCEWRR(Tid, NseqO, os_col, &obj[0]);
	  TCEWRR(Tid, NseqO, o1_col, &obj[1]);
	  TCEWRR(Tid, NseqO, o2_col, &obj[2]);
	  TCEWRR(Tid, NseqO, ni_col, &obj[3]);
	  if (toupper (qualif[0]) == 'S')
	    TCEWRC(Tid, NseqO, ot_col, &o_typ[nslit]);
	  NseqO++;  
        }
     }                 
}
/*---------------------------------------------------------------------------*/
#ifdef __STDC__
void def_sky(float sval[], float sky[], int nslit, int slit[])
#else
void def_sky( sval, sky, nslit, slit )
  float sval[],sky[];
  int   nslit,slit[];
#endif

{
   int   skystrt,skyend,nsky,j,slitend,found;
   char  text[80];

/*   slitend = upper[nslit]-2-Iwin; */
   slitend = upper[nslit]-(min_dist-1)-Iwin; 
   skystrt = 0.0;
   skyend = 0.0;
   nsky = 0;
   found = 0;

/*   sprintf(text, "Now searching sky in slitlet nr %d", slit[nslit]);
   SCTPUT(text);  */
   for (j=lower[nslit]-1+Iwin;j<=upper[nslit]-(min_dist-1)-Iwin;j++) 
     {
      if (sval[j] > MARK && j < slitend)
        {
         if (nsky == 0) skystrt = j;
         nsky++;    
        }
      else if ((sval[j] <= MARK || j == slitend) && nsky!=0)
        {
         skyend = j;
         nsky++;
/*         if (nsky > 5) */
         if (nsky > min_sky) 
           {
            sky[0] = skystrt;
            sky[1] = skyend;
            sky[2] = slit[nslit];
            TCEWRR(Tid, NseqS, s1_col, &sky[0]);
            TCEWRR(Tid, NseqS, s2_col, &sky[1]);
            TCEWRR(Tid, NseqS, ss_col, &sky[2]);
	    NseqS++;
            nsky = 0;
            found = 1;
           }
/*         else if (nsky <= 5) nsky = 0; */
         else if (nsky <= min_sky) nsky = 0;
        }
     }
/*   if (nsky <= 5 && found == 0) */
   if (nsky <= min_sky && found == 0)
     {
      sprintf(text, "No sky found in slitlet %d", slit[nslit]+1);
      SCTPUT(text);
     }
}

/*----------------------------------------------------------------------------*/
#ifdef __STDC__
float own_median(int n, float arr[], float *med)
#else
float own_median(n, arr, med)
int n;
float arr[],*med;
#endif
/*              find median of distribution            */
{
   int   i,j,index;
   float a,b[21];

   for (j=0;j<n;j++) b[j] = arr[j];

   for (j=1;j<n;j++) 
     {
      a=b[j];
      i=j-1;
      while (i >= 0 && b[i] > a) 
        {
         b[i+1]=b[i];
         i--;
        }
      b[i+1]=a;
     }
   index = (n-1)/2;
   *med = b[index];
   return *med;
}