/* @(#)mosslin.c	8.10 (ESO-IPG) 11/30/94 15:58:36 */
/* @(#)mosslin.c      1.2 (ESO-IPG) 2/11/94 12:11:46                       */
/*                                                                         */
/*           Search spectral lines with median: (val - med(val)) > thres   */
/*           Center lines: center of gravity or Gaussian fitting           */
/*                                                                         */
/*           Otmar Stahl, Aug. 17, 1993                                    */


#include <tbldef.h>
#include <midas_def.h>
#include <ok.h>
#include <proto_mutil.h>
#include <proto_nrutil.h>
#include <proto_fitnol.h>

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

#define MAXLINES 1000
#define MAXSLITS 100

#define GRAVITY 0
#define GAUSSIAN 1

double Start[2], Step[2];
float Thres, Rnull;
int Npix[2], Method, Nseq, Inull;
int Ybin, Ystep, Iwin, Width, Tid;
int YbinTMP;
double *Xgaus, *Ygaus, *A, Dnull;

#ifdef __STDC__
  void center_lines (float[], float[], float[], int[],
		   int[], int[], int, int []);
  void search_lines (float[], int[], int *);
  void fit_lines (float[], int[], float[], int, int []);
#else
  void center_lines ();
  void search_lines ();
  void fit_lines ();
#endif

main ()

{
  int i, iav;
  int nax, kun, knul, tid2, imno;
  int col, ncol, nslits, acol, arow, nsort, param[4], select;
  int null[3], slit[MAXSLITS], upper[MAXSLITS], lower[MAXSLITS];

  int icol[5];

  char method[21], wlc_image[60], mos_table[60];
  char output_table[60], unit[64], inst[72], line[80], linem[80];

  float moscol[3];
  float *buff1, *buff2, *image;

/*         get into Midas and get parameters             */

  SCSPRO ("mosslin");

  strcpy(inst,""), strcpy(unit,"");

  SCKGETC ("IN_A", 1, 60, &iav, wlc_image);
  SCKGETC ("IN_B", 1, 60, &iav, mos_table);
  SCKGETC ("OUT_A", 1, 60, &iav, output_table);
  SCKRDI ("INPUTI", 1, 4, &iav, param, &kun, &knul);
  SCKGETC ("INPUTC", 1, 3, &iav, method);

/*          map input file                              */

  SCIGET (wlc_image, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE,
	  2, &nax, Npix, Start, Step, inst, unit, (char **)&image, &imno);

  Method = GRAVITY;
  sprintf (linem, "centering method:    Gravity");
  if (!strncmp (method, "GAU", 3) || !strncmp (method, "gau", 3))
    {
      Method = GAUSSIAN;
      sprintf (linem, "centering method:    Gaussian");
    }
/*           read table MOS                            */

  TCTOPN (mos_table, F_I_MODE, &tid2);
  TCIGET (tid2, &col, &ncol, &nsort, &acol, &arow);
  TCCSER (tid2, ":slit",&icol[0]);
  TCCSER (tid2, ":ystart",&icol[1]);
  TCCSER (tid2, ":yend",&icol[2]);

  nslits = 0;
  for (i = 1; i <= ncol; i++)
    {
      TCSGET(tid2, i, &select);
      if (select) {
	TCRRDR (tid2, i, 3, icol, moscol, null);
	slit[nslits] = (int) moscol[0];
	lower[nslits] =  (int) ((moscol[1]-Start[1])/Step[1]);
	upper[nslits] =  (int) ((moscol[2]-Start[1])/Step[1]) + 1;
	nslits++;
      }
    }

  TCTCLO (tid2);

/*          open output table                         */

  TCTINI (output_table, F_TRANS, F_O_MODE, 5, MAXLINES, &Tid);
  SCDWRD (Tid, "Pixel", Step, 1, 1, &kun);
  TCCINI (Tid, D_R4_FORMAT, 1, "F10.2", "Pixel", "X", &icol[0]);
  TCCINI (Tid, D_R4_FORMAT, 1, "F10.2", "Pixel", "Y", &icol[1]);
  TCCINI (Tid, D_R4_FORMAT, 1, "E12.3", "Pixel", "Peak", &icol[2]);
  TCCINI (Tid, D_R4_FORMAT, 1, "F6.0", "None ", "Slit", &icol[3]);

  Thres = (float) param[0];
  Iwin  = param[1]/2;
  Width = 2*Iwin+1;
  Ystep = param[2];
  Ybin  = param[3]*2+1;

  SCTPUT ("search lines ");
  SCTPUT ("------------\n");
  sprintf (line, "Input image:         %s ", wlc_image);
  SCTPUT (line);
  sprintf (line, "Input table:         %s ", mos_table);
  SCTPUT (line);
  sprintf (line, "Output table:        %s\n ", output_table);
  SCTPUT (line);
  SCTPUT ("input parameters: ");
  sprintf (line, "search window:       %i pixels", Width);
  SCTPUT (line);
  sprintf (line, "detection threshold: %6.2f DN", Thres);
  SCTPUT (line);
  SCTPUT (linem);

  sprintf (line, "\naverage on:          %i scan lines", Ybin);
  SCTPUT (line);
  sprintf (line, "step:                %i scan lines\n", Ystep);
  SCTPUT (line);

  buff1  = (float *) osmmget (Npix[0] * sizeof (float));
  buff2  = (float *) osmmget (Npix[0] * sizeof (float));
  Xgaus = dvector (1, Width);
  Ygaus = dvector (1, Width);
  A     = dvector (1, 3);

  TCMNUL(&Inull,&Rnull,&Dnull);

  center_lines (image, buff1, buff2, slit, upper, lower, nslits, icol);

  sprintf (line, "\nTotal no of detections: %i lines\n", Nseq);
  SCTPUT (line);

/*          epilogue ...                 */

  TCSINI (Tid);
  TCTCLO (Tid);

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

  SCSEPI ();

}

/*---------------------------------------------------------------------------*/
#ifdef __STDC__
  void 
  center_lines (float image[], float rval[], float mval[], int slit[],
       	        int upper[], int lower[], int nslits, int icol[])
#else
  void 
  center_lines ( image, rval, mval, slit, upper, lower, nslits, icol)
  float image[],rval[], mval[];
  int   slit[],upper[],lower[],nslits,icol[];
#endif
{

  float obj[6];
  /*float buff[5];  replaced by TS-0896*/
  float *buff;
  int nact, nslit, ii, middle;
  int j, jj, k, linelist[MAXLINES];
  int iord;
  char line[80];

  Nseq = 1;

  if (Npix[1] > 1)
    {
      SCTPUT (" slit no.  detected lines ");
      SCTPUT (" --------  -------------- ");
    }

/*          loop over all slitlets                         */

  for (nslit = 0; nslit < nslits ; nslit++)
    {
      obj[3] = slit[nslit];

      /*        median Ybin rows in middle                */
      /*        changed from 5 rows to Ybin rows  TS-0896 */
      for (j = 0; j < Npix[0]; j++)
	rval[j] = 0.0;

      if (Ybin <= upper[nslit]-lower[nslit]+1) 
	{
	  middle = (int) (lower[nslit]+upper[nslit])/2 - (Ybin-1)/2; 
	  YbinTMP = Ybin;
	}
      else
	{
	  middle = 0;
	  YbinTMP = upper[nslit]-lower[nslit]+1;
	  sprintf(line," warning:  Ybin > slitlet");
	  SCTPUT(line);
	}
      buff = (float *) osmmget((YbinTMP+1)*sizeof(float));

      for (k = 0; k < Npix[0]; k++)
	{
	  for (j = middle ; j < middle + YbinTMP ; j++)
	    {
	      jj = j * Npix[0] + k;
	      ii = j-middle+1;
	      buff[ii] = image[jj];
	    }
	  mval[k] = pik_median(YbinTMP,buff);
	}
      
      /*       now search for lines   ...  */

      search_lines (mval, linelist, &nact);

      sprintf (line, "    %4i            %4i", slit[nslit], nact);
      SCTPUT (line);

      /*      now loop through the rows of one slitlet           */
      /*                  note that Ybin = 2 * radius + 1        */
      for (iord = lower[nslit]; iord <= upper[nslit]-YbinTMP+1; iord += Ystep)
	{
	  obj[1] =  (float) (iord + (YbinTMP-1)/2.) * Step[1] + Start[1];
	  /*      average scan  lines                                */

	  for (j = 0; j < Npix[0]; j++)
	    rval[j] = 0.0;
	  for (j = iord; j < iord + YbinTMP; j++)
	    {
	      for (k = 0; k < Npix[0]; k++)
		{
		  jj = j * Npix[0] + k;
		  rval[k] = rval[k] +  image[jj];
		}
	    }
	  for (k = 0; k < Npix[0]; k++) 
	    {
	      rval[k] = rval[k] / YbinTMP;
	    }

	  fit_lines (rval, linelist, obj, nact, icol);
	}
    }

  SCTPUT (" ----------------------- ");
}

/*---------------------------------------------------------------------------*/
#ifdef __STDC__
  void 
  search_lines (float rval[], int linelist[], int *nact)
#else
  void 
  search_lines ( rval, linelist, nact )
  float rval[];
  int linelist[],*nact;
#endif

{

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

  for (j = Iwin; j < Npix[0] - Iwin; j++)
    {
      diff = rval[j] - pik_median (Width, &rval[j - Iwin]);
      if (diff > Thres)
	{
	  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;
	      }
	  linelist[number] = imax;
	  number++;
	}
    }

  /*       remove doublets                                  */

  i = 0;
  while (i < number - 1)
    {
      if ((linelist[i + 1] - linelist[i]) < Width)
	{
	  number--;
	  for (j = i + 1; j < number; j++)
	    {
	      linelist[j] = linelist[j + 1];
	    }
	  i--;			/* check again (more than two detections!) */
	}
      i++;
    }
  *nact = number;
}

/*---------------------------------------------------------------------------*/
#ifdef __STDC__
  void 
  fit_lines (float rval[], int linelist[], float obj[], int nact, int icol[])
#else
  void 
  fit_lines ( rval, linelist, obj, nact, icol )
  float rval[],obj[];
  int   linelist[],nact,icol[];
#endif

{

  int j, k, jj, index;
  float factor, shift, a, b, aleft, aright, minimum;
  double x_guess;

  for (jj = 0; jj < nact; jj++)
    {
      index = linelist[jj];

      /*  center of gravity                                */

      switch (Method)
	{

	case GRAVITY:

	  aleft = rval[index - 1];
	  aright = rval[index + 1];
	  factor = 1;
	  if (aleft >= aright)
	    {
	      aleft = rval[index + 1];
	      aright = rval[index - 1];
	      factor = -1;
	    }
	  a = rval[index] - aleft;
	  b = aright - aleft;
	  shift = (a + b == 0.0) ? 0.0 : Step[0] * b / (a + b);
	  obj[0] = Start[0] + Step[0] * linelist[jj] + factor * shift;
	  obj[2] = rval[index];
	  break;

	case GAUSSIAN:

	  for (k = 1, j = index - Iwin; j <= index + Iwin; j++, k++)
	    {
	      Xgaus[k] = Start[0] + Step[0] * linelist[jj] + k - Iwin - 1;
	      Ygaus[k] = rval[j];
	    }

	    /* subtract minimum */

	  minimum = 1.e34;
	  for (k = 1 ; k <= 2 * Iwin + 1; k++)
	    {
	      if(Ygaus[k] < minimum)
		{
		  minimum = Ygaus[k];
		}
	    }
	  for (k = 1 ; k <= 2 * Iwin + 1; k++)
	    {
	      Ygaus[k] = Ygaus[k]-minimum;
	    }

	  A[1] = Ygaus[Iwin+1];
	  A[2] = Start[0] + Step[0] * linelist[jj];
	  x_guess = A[2];
	  A[3] = 3*Step[0];


	  /*	  
		  for (k = 1 ; k <= 2 * Iwin + 1; k++)
		  {
		  if (Xgaus[k] > 0 && Xgaus[k] < 500) 
		  {
		  printf("%f %f %f %f %f\n",Xgaus[k],Ygaus[k],A[1],A[2],A[3]);
		  } 
		  }
	  */
	    
	  
	  fit_gauss (Xgaus, Ygaus, Width, A);
	  
	  obj[0] = A[2];
	  obj[2] = A[1];

	  /* check if result out of limits */

	  if (myabs (A[2] - x_guess)  > Iwin) {
	    obj[0] = Rnull;
	    obj[2] = Rnull;
	  }
	  break;
	}

      TCRWRR (Tid, Nseq, 4, icol, obj);
      Nseq++;
    }
}