/* @(#)mo_offset.c	17.1.1.1 (ES0-DMD) 01/25/02 17:49: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
===========================================================================*/

/* @(#)mo_offset.c	17.1.1.1 (ESO-DAG) 01/25/02 17:49:07 */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENTIFICATION: routine mo_offset
.KEYWORDS:       bulk data frames, mosaicing, offset
.PURPOSE:        Construct a mosaic of a number of image frames, the result
                 will a frame with nx subframes in x and ny subframes in y.
                 The subframes must have the same size in x and y.
                 Blank subframes can be indicated on the command line.
.ALGORITHM:      Finds the relative offsets in background intensity between
                 every pair of overlapping frames in a mosaic.
                 The format of the command line is : findoff input_file 
                 min_pixels.
.VERSION:        950904 RHW Created; original from Mike Regan (Univ. MarylanD)
-----------------------------------------------------------------------------*/
/*
 * Define _POSIX_SOURCE to indicate
 * that this is a POSIX program
 */
#define _POSIX_SOURCE 1

/*
 * definition of the used functions
 */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <ftoc.h>
#include <midas_def.h>
#include <computer.h>
#include <ccd_def.h>

#define NCOLIN     8
#define NCOLOUT    4

/*
 ++++++++++++++++++++++++++++++++++++++++++++++++++
 * here starts the code of the function
 ---------------------------------------------------
*/

MO_OFFSET(inim, inmsk, tid, oc1, ol1, npair, offset, xcount, 
          iref, jref, mnpx, verb) 

int     inim;
int     inmsk;
int     tid;   
int     *oc1;
int     *ol1;
int     *npair;
float   *offset;
float   *xcount;
int     *iref;
int     *jref;
int     mnpx;
int     verb;

{
char   filename[MAXFRM][61];
char   line[80];

float  *cube;
float  *mask;
float  deltax, deltay;
float  image[4];  
float  reloff[MAXFRM][MAXFRM];
float  vec1[MAXSIZ];
float  vec2[MAXSIZ];
float  xoff[MAXFRM];
float  yoff[MAXFRM];
float  med1, med2, moffset;

int    colin[NCOLIN];
int    count1, count11;
int    count2, count22;
int    ex1, ex2;
int    ey1, ey2;
int    idata[3];
int    imnom;
int    i, ii, ic, iav, irow;
int    indexi, indexj, indexjj, indexii;
int    ismoot = 0;
int    j, jj;
int    k;
int    mask1, mask2;
int    numfiles;
int    npixx, npixy;
int    null;
int    nulo;
int    ncoldb, nrowdb, nsortdb, allcoldb, allrowdb;
int    naxis, npix[3];
int    naxism, npixm[3];
int    size[MAXFRM],sizem;
int    stat;
int    sx1, sx2;
int    sy1, sy2;
int    uni;
int    xsta[MAXFRM], xend[MAXFRM];
int    ysta[MAXFRM], yend[MAXFRM];

double start[3], step[3], startm[3], stepm[3];

static char  im_col[]       = "In_frame";
static char  xsta_col[]     = "XStartpix";
static char  xend_col[]     = "XEndpix";
static char  ysta_col[]     = "Ystartpix";
static char  yend_col[]     = "YEndpix";
static char  xoff_col[]     = "X_offpix";
static char  yoff_col[]     = "Y_offpix";
static char  ioff_col[]     = "Offset";

/* 
 Get some data
 */
(void) TCIGET(tid, &ncoldb, &nrowdb, &nsortdb, &allcoldb, &allrowdb);
TCCSER(tid, im_col, &colin[0]);
TCCSER(tid, xsta_col, &colin[1]);
TCCSER(tid, xend_col, &colin[2]);
TCCSER(tid, ysta_col, &colin[3]);
TCCSER(tid, yend_col, &colin[4]);
TCCSER(tid, xoff_col, &colin[5]);
TCCSER(tid, yoff_col, &colin[6]);
TCCSER(tid, ioff_col, &colin[7]);

stat  = SCDRDI(inim,"NPIX",1,3,&iav,npix,&uni,&nulo);

/*
 Allocate memory for the arrays
 */
cube  = (float *) osmmget(MAXSIZ * nrowdb * sizeof( float ));
mask  = (float *) osmmget(MAXSIZ * sizeof( float ));

/*
 Run through the subrasters and store data into cube
 */
numfiles = 0;
for (i = 0; i < nrowdb; i++)
   {
   indexi = i * MAXSIZ;
   irow = i + 1;
   if (strncmp(filename[i],"null",4) == 0) 
      *(cube+indexi) = MO_NULL;

   else
      {
      stat =  TCERDC(tid, irow, colin[0], filename[i], &null);
      stat =  TCERDI(tid, irow, colin[1], &xsta[i], &null);
      stat =  TCERDI(tid, irow, colin[2], &xend[i], &null);
      stat =  TCERDI(tid, irow, colin[3], &ysta[i], &null);
      stat =  TCERDI(tid, irow, colin[4], &yend[i], &null);
      stat =  TCERDR(tid, irow, colin[5], &xoff[i], &null);
      stat =  TCERDR(tid, irow, colin[6], &yoff[i], &null);

      image[0] = (float) xsta[i];
      image[1] = (float) xend[i];
      image[2] = (float) ysta[i];
      image[3] = (float) yend[i];
      npixx    = xend[i] - xsta[i] + 1;
      npixy    = xend[i] - xsta[i] + 1;
      size[i]  = npixx * npixy;
      if (i > 0 && size[i] != size[0]) 
         SCETER(2,"*** FATAL: Subrasters have unequal sizes");
      GETDAT(inim, MAXSIZ, npix, image, ismoot, &cube[indexi]);
      numfiles++;
      }
   }

/*
 Get the mask 
 */
if (inmsk == -1)
   for (j = 0; j < MAXSIZ; j++)
      *(mask+j) = 1.0;

else
   {
   stat  = SCDRDI(inmsk,"NAXIS",1,1,&iav,&naxism,&uni,&nulo);
   stat  = SCDRDI(inmsk,"NPIX",1,3,&iav,npixm,&uni,&nulo);
   stat  = SCDRDD(inmsk,"START",1,3,&iav,startm,&uni,&nulo);
   stat  = SCDRDD(inmsk,"STEP",1,3,&iav,stepm,&uni,&nulo);
   image[0] = 1.0;
   image[1] = 1.0;
   image[2] = (float) npixm[0];
   image[3] = (float) npixm[1];
   sizem    = npixm[0] * npixm[1];
   if (sizem != size[0]) 
         SCETER(2,"*** FATAL: Mask does have different size");
   GETDAT(inmsk, MAXSIZ, npixm, image, ismoot, mask);
   }

/* 
 List the input subrasters
 */
if (verb == 1)
   {
   sprintf(line,
       "Filename                    in_area             xoffset      yoffset");
   SCTPUT(line);
   for (i = 0; i < nrowdb; i++)
       {
      if (strncmp(filename[i],"null",4) != 0) 
         {
          sprintf(line, "%-20s   [%4d,%4d:%4d,%4d]   %6.1f       %6.1f", 
                 filename[i], xsta[i], ysta[i], xend[i], yend[i], 
                 xoff[i], yoff[i]);
          SCTPUT(line);
	 }
       }
   sprintf(line,"Total number of input subrasters: %d", numfiles);
   SCTPUT(line);
   SCTPUT(" ");
   }

else
   {
   sprintf(line,"Total number of input subrasters: %d", numfiles);
   SCTPUT(line);
   SCTPUT(" ");
   }

/* 
 Loop through the input frame subrasters
 */
sprintf(line,"frame_1   frame_2    offset   #x_pix   #y_pix");
   SCTPUT(line);
for (i = 0; i < nrowdb; i++)
   {
   if (strncmp(filename[i],"null",4) == 0) 
      continue;

   indexi = i * MAXSIZ;
   for (j = i+1; j < nrowdb; j++)
      {
      if (strncmp(filename[j],"null",4) == 0) 
      continue;

      indexj = j * MAXSIZ;
      deltax = (float) (oc1[i]-oc1[j]);               /* x offset in pixels */
      deltay = (float) (ol1[i]-ol1[j]);               /* y offset in pixels */
      if (fabs (deltax) < (float)(npixx-4) && 
          fabs(deltay)  < (float)(npixy-4))
         {
         if (deltax < 0)
            {
            sx1 = 1 + fabs(deltax);
            ex1 = npixx;
            sx2 = 1;
            ex2 = npixx + deltax;
	    }

         else
            {
            sx1 = 1;
            ex1 = npixx - deltax;
            sx2 = 1 + deltax;
            ex2 = npixx;
            }

         if (deltay < 0)
            {
            sy1 = 1+fabs(deltay);
            ey1 = npixy;
            sy2 = 1;
            ey2 = npixy+deltay;
            }

         else
            {
            sy1 = 1;
            ey1 = npixy - deltay;
            sy2 = 1 + deltay;
            ey2 = npixy; 
            }

         count1   = 0;
         count11  = 0;

         for (jj = sy1; jj <= ey1; jj++)
            {
            indexjj = (jj-1)*npixx;
            for (ii = sx1; ii <= ex1; ii++)
               {
               indexii = ii-1;
               mask1   = *(mask + indexjj + indexii);
               mask2   = *(mask + indexjj + (int)deltay*npixx 
                                + indexii + (int)deltax);
               if ((mask1 == 1.0) && (mask2 == 1.0))
                  {
                  vec1[count1] = *(cube + indexi + indexjj + indexii);
                  count1++;
		  }
               else
                  count11++;
               }
	   }

         count2  = 0;
         count22 = 0;
             
         for (jj = sy2; jj <= ey2; jj++)
            {
            indexjj  = (jj-1)*npixx;
            for (ii = sx2; ii <= ex2; ii++)
               {
               indexii = ii-1;
               mask1   = *(mask + indexjj + indexii);
               mask2   = *(mask + indexjj - (int)deltay*npixx 
                                + indexii - (int)deltax);
               if ((mask1 == 1.0) && (mask2 == 1.0))
                  {
                  vec2[count2] = *(cube + indexj + indexjj + indexii);
                  count2++;
		  }
               else
                  count22++;
               }
	    }

         if (count1 > mnpx) 
            {
            MO_MEDM(count1,vec1,&med1);
            MO_MEDM(count2,vec2,&med2);
            moffset = med1-med2;
            sprintf(line, "%3d       %3d   %9.3g   %6d   %6d", 
                             i+1, j+1, moffset, count1, count2);
            SCTPUT(line);
	    }

         else
            moffset = 0.0;
         }

      else
         moffset = 0.0;

      reloff[i][j] = moffset;
      }
   }


ic = 1;
for (i = 0; i < nrowdb; i++)
   {
   for (j = i+1; j < nrowdb; j++)
      {
      *(offset+ic) = reloff[i][j];
      if (*(offset+ic) != 0) 
         {
         *(xcount+ic) = (float) ic;
         *(iref+ic)   = i+1;
         *(jref+ic)   = j+1;
         ic++;
         }
      }
   }
sprintf(line, "Number of pairs included in offset calculations: %3d",
               ic-1);
SCTPUT(line);
SCTPUT(" ");

*npair = ic-1;
}


/*
 MO_MEDM -- Routine to compute the median value of an array
 */
MO_MEDM(n, a, amed)

int        n;
float      *a;
float      *amed;

{
int        i, imed, m;
float      denom;
double     *w;
char       *osmmget();

w = (double *) osmmget(n * sizeof( double ));
m = 0;
for (i = 0; i < n; i++)
   {
   w[i] = a[i];
   m++;
   }

if (m == 0)
   *amed = 0.0;
else if (m == 1) 
   *amed = *a;
else
   {
   (void) sortd(n,w);
   imed  = m/2 + 1;
   denom = 2.0;
   if (fmod((float) n, denom) == 0) 
      *amed = .5 * (w[imed] + w[imed-1]);
   else
      *amed = w[imed];
   }
osmmfree((char *) w);
}