/*                                         */ 
/* Search slitlets with gradient technique */ 
/* Center lines: center of gravity         */ 
/*                                         */ 
/* Otmar Stahl,              Aug. 17, 1994 */
/*                                         */ 
/* Adapt for SPOL (additional column for   */
/* object type, offsets always valid       */
/* for two slitlets)                       */
/*                                         */ 
/* Sabine Moehler Jan. 3, 1995             */

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

#define MAXSLITS 1000

double Start[2], Step[2];
float Thres;
int Npix[2], Nseq;
int Iwin, Width;

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

main ()

{
  int i, jj, actval;
  int nax, kun, knul, tid, imno;
  int ystart[MAXSLITS], yend[MAXSLITS];

  int icol[5], bytelem, noelem, noman;

  char qualif[4], ff_image[60], mos_table[60];
  char original[60], unit[64], instruma[72], line[80], type;

  float tab[5], offset[MAXSLITS], param[4], low;
  float *buff, *image;

/*         get into Midas and get parameters             */

  SCSPRO ("mosfind");

  strcpy(instruma,""); 
  strcpy(unit,"");

  SCKGETC ("QUALIF",  1, 4, &actval, qualif);
  SCKGETC ("ORI",   1, 60, &actval, original);
  SCKGETC ("IN_A",   1, 60, &actval, ff_image);
  SCKGETC ("OUT_A",  1, 60, &actval, mos_table);
  SCKRDR  ("INPUTR", 1, 4, &actval, param, &kun, &knul);
  SCKRDI  ("INPUTI", 1, 1, &actval, &noman, &kun, &knul);
/*          map 1-dimensional input file                             */

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

  SCDFND (imno, "XPOS", &type, &noelem, &bytelem);
  if (noelem == 0)
    {
      sprintf(line,"Descriptor XPOS not found, set to zero \n");
      SCTPUT (line);
    }
  else
    {
      sprintf(line,"Descriptor XPOS found, read it\n");
      SCTPUT (line);
      SCDRDR (imno, "XPOS", 1, 100, &actval, offset, &kun, &knul);
    }

  Thres = param[0];
  Iwin  = (int) param[1]/2;
  Width = 2*Iwin+1;
  low = param[2];

  SCTPUT ("search slitlets ");
  SCTPUT ("------------\n");
  sprintf (line, "Input image:         %s ", original);
  SCTPUT (line);
  sprintf (line, "Output table:        %s\n ", mos_table);
  SCTPUT (line);
  SCTPUT ("input parameters: ");
  sprintf (line, "search window:       %i pixels", Width);
  SCTPUT (line);
  sprintf (line, "detection threshold: %6.2f ", Thres);
  SCTPUT (line);
  sprintf (line, "no mans land at edge of slitlet: %i ", noman);
  SCTPUT (line);

  buff  = (float *) osmmget (Npix[0] * sizeof (float));

/* do the real work */

  center_lines (image, buff, ystart, yend, low); 

  SCTPUT (" ----------------------- ");
  sprintf (line, "\nTotal no of detections: %i slitlets\n", Nseq);
  SCTPUT (line);
 
/* write output table MOS */

  if (Nseq == 0)
    {
      SCTPUT ("Slit detection failed");
      if (Iwin < 1) 
	SCTPUT ("width must be >= 2."); 
      else 
	SCTPUT ("Use a lower or higher threshold");
      SCKWRI ("NSLIT", &Nseq, 1, 1, &kun);
      SCSEPI ();
      exit (1);
    }

  SCKWRI ("NSLIT", &Nseq, 1, 1, &kun);
  sprintf (line, "  y_start         y_end");
  SCTPUT (line);
  sprintf (line, "  -------         -----");
  SCTPUT (line);

  switch (toupper(qualif[0]))
    {
     case ('M') :
       {
        if (TCTINI (mos_table, F_TRANS, F_O_MODE, 4, Nseq, &tid))
          SCTPUT ("Error creating table");

        TCCINI (tid, D_R4_FORMAT, 1, "F6.0", "unitless", "slit", &icol[0]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "ystart", &icol[1]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "yend", &icol[2]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "xoffset", &icol[3]);
        for (i = 1; i < Nseq + 1; i++)
          {
           jj = i - 1;
           tab[0] = i;
           tab[1] = Start[0] + (ystart[jj] + noman - 1) * Step[0];
           tab[2] = Start[0] + (yend[jj] - noman - 1) * Step[0];
	   if ( (tab[2]-tab[1]) > 0.5) 
	     sprintf (line, " %9.2f     %9.2f", tab[1], tab[2]);
	   else
	     sprintf (line, " %9.2f     %9.2f  WARNING: y_start = y_end", 
		      tab[1], tab[2]);
           SCTPUT (line);
           tab[3] = offset[jj];
           TCRWRR (tid, i, 4, icol, tab);
          }
        break;
       }
     case ('S') :
       {
        if (TCTINI (mos_table, F_TRANS, F_O_MODE, 5, Nseq, &tid))
          SCTPUT ("Error creating table");

        TCCINI (tid, D_R4_FORMAT, 1, "F6.0", "unitless", "slit", &icol[0]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "ystart", &icol[1]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "yend", &icol[2]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "xoffset", &icol[3]);
        /* Additional column ray_typ (O(rdinary) or E(xtra-ordinary)) */
        TCCINI (tid, D_C_FORMAT, 1, "A1", "", "RAY_TYP", &icol[4]);
        /* The offsets that are read in are always valid for two slitlets due to
           the splitting of the incoming light */

        for (i = 1; i < Nseq + 1; i=i+2)
          {
           jj = i - 1;
           tab[0] = i;
           tab[1] = Start[0] + (ystart[jj] - 1) * Step[0];
           tab[2] = Start[0] + (yend[jj] - 1) * Step[0];
           sprintf (line, " %9.2f     %9.2f", tab[1], tab[2]);
           SCTPUT (line);
           tab[3] = offset[jj/2];
           TCRWRR (tid, i, 4, icol, tab);
          /* Additional column ray_typ (O(rdinary)) */
           TCEWRC (tid, i, icol[4], "O");
           jj = i;
           tab[0] = i+1;
           tab[1] = Start[0] + (ystart[jj] - 1) * Step[0];
           tab[2] = Start[0] + (yend[jj] - 1) * Step[0];
           sprintf (line, " %9.2f     %9.2f", tab[1], tab[2]);
           SCTPUT (line);
           tab[3] = offset[(jj-1)/2];
           TCRWRR (tid, i+1, 4, icol, tab);
           /* Additional column ray_typ (E(xtra-ordinary)) */
           TCEWRC (tid, i+1, icol[4], "E");
          }
        break;
       }
     case ('I') :
       {
        if (TCTINI (mos_table, F_TRANS, F_O_MODE, 4, Nseq, &tid))
          SCTPUT ("Error creating table");

        TCCINI (tid, D_R4_FORMAT, 1, "F6.0", "unitless", "slit", &icol[0]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "ystart", &icol[1]);
        TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "Pixel", "yend", &icol[2]);
        /* Additional column ray_typ (O(rdinary) or E(xtra-ordinary)) */
        TCCINI (tid, D_C_FORMAT, 1, "A1", "", "RAY_TYP", &icol[4]);

        for (i = 1; i < Nseq + 1; i=i+2)
          {
           jj = i - 1;
           tab[0] = i;
           tab[1] = Start[0] + (ystart[jj] - 1) * Step[0];
           tab[2] = Start[0] + (yend[jj] - 1) * Step[0];
           sprintf (line, " %9.2f     %9.2f", tab[1], tab[2]);
           SCTPUT (line);
           TCRWRR (tid, i, 3, icol, tab);
          /* Additional column ray_typ (O(rdinary)) */
           TCEWRC (tid, i, icol[3], "O");
           jj = i;
           tab[0] = i+1;
           tab[1] = Start[0] + (ystart[jj] - 1) * Step[0];
           tab[2] = Start[0] + (yend[jj] - 1) * Step[0];
           sprintf (line, " %9.2f     %9.2f", tab[1], tab[2]);
           SCTPUT (line);
           TCRWRR (tid, i+1, 3, icol, tab);
           /* Additional column ray_typ (E(xtra-ordinary)) */
           TCEWRC (tid, i+1, icol[3], "E");
          }
        break;
       }
     default:
       {
        sprintf(line,"Error: Unknown qualifier %s\n",qualif);
        SCTPUT(line);
        SCETER(9, "Exiting...");
        break;
       }
    }

  strcpy (instruma, "Slit positions");
  SCDWRC (tid, "IDENT", 1, instruma, 1, 72, &kun);

/*          epilogue ...                 */

  TCSINI (tid);
  TCTCLO (tid);
  
  osmmfree((char *) buff);

  SCSEPI ();

}

/*---------------------------------------------------------------------------*/
#ifdef __STDC__
  void 
  center_lines (float image[], float rval[], int ystart[], int yend[], 
		float low)
#else
  void 
  center_lines ( image, rval, ystart, yend, low)
  float image[], rval[], low;
  int   ystart[], yend[];
#endif
{

  float xpos[MAXSLITS];
  int i, nact, linelist[MAXSLITS];

  Nseq = 0;
  
  for (i = 2; i < Npix[0]-2; i++)    {
    rval[i] = (image[i]<low) ? 0.0 :  (image[i+1]-image[i-1]) / image[i];
  }

  search_lines (rval, linelist, &nact); 
  fit_lines (rval, linelist, xpos, nact); 
  
  for (i = 0 ; i< Nseq ; i+=2) {
    ystart[i/2] = xpos[i];
    yend[i/2] = xpos[i+1];
  }

  Nseq = (Nseq+1)/2;

}

/*---------------------------------------------------------------------------*/
#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, mode;

  mode = 1;

  for (j = Iwin; j < (Npix[0]-1) - Iwin; j++)
    {
      diff = (float) mode * (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;
	  mode = -mode;
	  number++;
	}
    }

  *nact = number; 
}

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

{

  int  jj, index;
  float factor, shift, a, b, aleft, aright;

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

      /*  center of 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[jj] = Start[0] + Step[0] * linelist[jj] + factor * shift;
      
      Nseq++;
    }
}