/* @(#)search.c	17.1.1.1 (ES0-DMD) 01/25/02 17:54:07 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT     (c)  1994  European Southern Observatory
.LANGUAGE      C
.IDENT         search
.AUTHOR        Preben Grosbol  (ESO-IPG)
.KEYWORDS      Object search
.ENVIRONMENT   MIDAS
.COMMENTS      Routines used for Measuring Machine facility
.PURPOSE       Fast search of objects in large frames
.VERSION       1.0   1988-Nov-18 :   Creation,     PJG
.VERSION       1.1   1994-Feb-11 :   Update to new SC-routines,     PJG
------------------------------------------------------------------------*/
#include   <math.h>                  /* mathematical macros             */
#include   <midas_def.h>             /* MIDAS definitions               */
#include   <mmfsrch.h>

main()
{
  char    *pc, *buf, meth, mflag[3], fname[41], tname[41], text[81];
  int     err, na, i, n, nx, ny, iva, inull;
  int     nnew, nold, iobj, ntp, ntl;
  int     idf, idx, idi, idt, lsd, lev, size, iy, mm, ioff, ip;
  int     npix[2], kunit[4], iparm[4];
  int     *pi, *ppi;
  float   *ppf;
  float   dif, adif, sden, bg, dmax, grad;
  float   x, y, dx, dy, xs, ys, xx, yy, xy, xp, yp, xpos, ypos, xoff, yoff;
  double  start[2], step[2];
  OBJ     *nobj, *oobj, newobj[MAX_OBJ], oldobj[MAX_OBJ];

  SCSPRO("SEARCH");                     /* initiate MIDAS environment   */

  /* read FRAME name and open it  */

  err = SCKGETC("IN_A", 1, 40, &iva, fname);
  err = SCFOPN(fname, D_R4_FORMAT, 0, 1, &idf);
  err = SCDRDI(idf, "NAXIS", 1, 1, &iva, &na, kunit, &inull);
  if (na!=2)                            /* search frame must be 2-dim.  */
     SCETER(1, "Error: Frame not 2-dimentional");
  err = SCDRDI(idf, "NPIX", 1, 2, &iva, npix, kunit, &inull);
  nx = npix[0]; ny = npix[1];
  err = SCDRDD(idf, "START", 1, 2, &iva, start, kunit, &inull);
  xs = start[0]; ys = start[1];
  err = SCDRDD(idf, "STEP", 1, 2, &iva, step, kunit, &inull);
  dx = step[0]; dy = step[1];

  /* read search method and parameters  */

  err = SCKGETC("INPUTC", 1, 1, &iva, mflag);
  meth = mflag[0];

  err = SCKRDI("INPUTI", 1, 3, &iva, iparm, kunit, &inull);
  switch (meth) {                   /* define parameters for methods    */
     case 'L' : mm = iparm[0]; lev = iparm[1];
		xoff = 0.5*(mm+1) - 2.0; yoff = -1.0;
		break;
     case 'G' : mm = iparm[0]; lev = iparm[1]; lsd = iparm[2];
		xoff = 0.5*(mm+1) - 2.0; yoff = -1.0;
		break;
     case 'P' : mm = iparm[0]; lev = iparm[1];
		xoff = 0.5*(mm+1) - 2.0; yoff = -1.0;
		break;
     default  : mm = iparm[0]; lev = iparm[1];
		xoff = 0.5*(mm+1) - 2.0; yoff = -1.0;
  }

  /* read table name and create it      */

  err = SCKGETC("OUT_A", 1, 40, &iva, tname);

  /* get memory and initiate variables  */

  err = SCFCRE("DUMR4", D_R4_FORMAT, F_X_MODE, 1, nx, &idx);
  err = SCFMAP(idx, F_X_MODE, 1, nx, &size, &pc);
  buf = pc;
  err = SCFCRE("DUMI4", D_I4_FORMAT, F_X_MODE, 1, nx, &idi);
  err = SCFMAP(idi, F_X_MODE, 1, nx, &size, &pc);
  pi  = (int *) pc;

  ioff = 1; iobj = 0; nold = 0;
  for (iy=0; iy<ny; iy++) {          /* search of objects line by line  */
      err = SCFGET(idf, ioff, nx, &size, buf);
      ppf = (float *) buf; ppi = pi; i = size;
      while (i--) *ppi++ = *ppf++;           /* convert from R4 to I4   */
      switch (meth) {                        /* chose detection method  */
	 case 'L' : nnew = slslv(pi,size,newobj,mm,lev);
		    break;
	 case 'G' : nnew = slsgd(pi,size,newobj,mm,lev,lsd);
		    break;
	 case 'P' : nnew = slslv(pi,size,newobj,mm,lev);
		    break;
	 default  : nnew = slslv(pi,size,newobj,mm,lev);
      }
      if (MAX_OBJ<nnew) {                   /* too many objects found   */
	 SCTPUT("Error: Too many objects in single line detected");
	 break;
      }
      i  = nold; oobj = oldobj;
      ip = nnew; nobj = newobj;
      while (i--) {                             /* compare with old obj. */
       	  while (ip) {
             dif = nobj->x - oobj->x/oobj->den + nobj->ix - oobj->ix;
             adif = 2.0 * ( (dif<0.0) ? -dif : dif );
             if (adif<nobj->npix) {                /* update object     */
                xp = nobj->ix - oobj->ix;
                xx = nobj->den * nobj->x;
                xy = nobj->den * xp;
                yp = oobj->nline + 1;
                nobj->x = oobj->x + xx + xy;
                nobj->y = oobj->y + nobj->y*yp;
                nobj->xs = oobj->xs + nobj->xs + xp*(2.0*xx+xy);
                nobj->ys = oobj->ys + yp*yp*nobj->ys;
                nobj->xy = oobj->xy + yp*(xy+xx);
                nobj->den += oobj->den;
                if (nobj->max<oobj->max) nobj->max = oobj->max;
                nobj->bg += oobj->bg;
                if (nobj->grad<oobj->grad) nobj->grad = oobj->grad;
                nobj->ix = oobj->ix;
                nobj->iy = oobj->iy;
                nobj->npix += oobj->npix;
                nobj->nline += oobj->nline;
                ip--; nobj++;
                goto nexti;
             }
             else if (0.0<dif) break;
             ip--; nobj++;
          }

          if (1<oobj->nline) {                     /* terminate object  */
             iobj++;
             xp = oobj->x / oobj->den;
             yp = oobj->y / oobj->den;
             xpos = dx*(xp+xoff+oobj->ix) + xs;
             ypos = dx*(yp+yoff+oobj->iy) + ys;
             xx = oobj->xs/oobj->den - xp*xp;
             xx = (xx<0.0) ? 0.0 : sqrt(xx);
             yy = oobj->ys/oobj->den - yp*yp;
             yy = (yy<0.0) ? 0.0 : sqrt(yy);
             xy = (0.0<xx && 0.0<yy) ?
                    (oobj->xy/oobj->den-xp*xp)/(xx*yy) : 0.0;
             xx *= 2.35482*dx; yy *= 2.35482*dy;
             sden = oobj->den;
             dmax = oobj->max;
             bg   = oobj->bg/oobj->nline;
             grad = oobj->grad;
             ntp  = oobj->npix;
             ntl  = oobj->nline;
             sprintf(text,
               "%5d:%9.1f %9.1f,%5.1f %5.1f,%9.1f %7.1f %6.1f, %6d %5d",
                iobj,xpos,ypos,xx,yy,sden,dmax,bg,ntp,ntl);
             SCTPUT(text);
          }
	  nexti: oobj++;
      }

      nold = nnew;                            /* swap newobj to oldobj  */
      nobj = newobj; oobj = oldobj;
      while (nnew--) {
        oobj->x = nobj->x;
        oobj->y = nobj->y;
        oobj->xs = nobj->xs;
        oobj->ys = nobj->ys;
        oobj->xy = nobj->xy;
        oobj->den = nobj->den;
        oobj->max = nobj->max;
        oobj->bg = nobj->bg;
        oobj->grad = nobj->grad;
        oobj->ix = nobj->ix;
        oobj->iy = nobj->iy;
        oobj->npix = nobj->npix;
        oobj->nline = nobj->nline;
	if (oobj->nline<2) {
	   oobj->x *= oobj->den;
	   oobj->iy = iy;
	}
        oobj++; nobj++;
     }
     ioff += size;
  }
  sprintf(text,"Total no. of objects found : %d",iobj);
  SCTPUT(text);

  SCFCLO(idf);                    /* finished - close frames and exit   */
  SCFCLO(idx);
  SCFCLO(idi);
  SCSEPI();
}