/* @(#)mo_links.c	17.1.1.1 (ES0-DMD) 01/25/02 17:49:06 */
/*===========================================================================
  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) 1995 European Southern Observatory
.IDENTIFIER  mo_links.c
.AUTHOR      R.H. Warmels IPG-ESO Garching
.KEYWORDS    alignment software
.LANGUAGE    C
.PURPOSE     Routine to compute the shift for each subframe
.ENVIRONment MIDAS
             #include <ccd_def.h>      Symbols used by the ccd package
.VERSION     1.0     16-May-1995   creation
------------------------------------------------------------*/
/*
 * Define _POSIX_SOURCE to indicate
 * that this is a POSIX program
 */
#define  _POSIX_SOURCE 1

/*
 * definition of the used functions in this module
 */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#include <ccd_def.h>

/*
 MO_LINK -- Routine to compute the shifts directly. Input is given
            by an input table containing the coordinates of common objects
 */

MO_LINKS(tblco, col, xrshift, yrshift, xcshift, ycshift, nrshift,
        ncshift, ncols, nrows, nxrsub, nyrsub, nxsub, nysub, nxoverlap,
        nyoverlap, order, nshifts)
int     tblco;
int     *col;                   /* columns containing the shifts */
float   (*xrshift)[MAXFRM];
float   (*yrshift)[MAXFRM];
float   (*xcshift)[MAXFRM];
float   (*ycshift)[MAXFRM];
int     (*nrshift)[MAXFRM];     /* number of row shifts */
int     (*ncshift)[MAXFRM];     /* number of column shifts */
int     ncols;                  /* number of columns per subraster */
int     nrows;                  /* number of rows per subraster */
int     nxrsub;                 /* column index of reference subraster */
int     nysub;                  /* number of subrasters in y */
int     nxoverlap;              /* number of columns of overlap */
int     nyoverlap;              /* number of rows of overlap */
char    *order;                 /* row or column order */
int     *nshifts;

{	
int     i, j, nxsize, nysize, ilimit, olimit;
int     ns;
float   *xcolavg, *ycolavg, *xrowavg, *yrowavg;
int     *nrowavg, *ncolavg;
float   isign, jsign, xrmed, yrmed, xcmed, ycmed;

if (strcmp(order,MO_COLUMN) == 0)
   {
   ilimit = nysub;
   olimit = nxsub;
   }
else 
   {
   ilimit = nxsub;
   olimit = nysub;
   }
        
/*
 Accumulate the shifts.
 */
nxsize  = ncols - nxoverlap;
nysize  = nrows - nyoverlap;
MO_GET_SHIFTS(tblco, col, xrshift, yrshift, nrshift, xcshift,
                        ycshift, ncshift, nxsub, nysub, nxrsub, nyrsub, 
                        nxoverlap, nyoverlap, nxsize, nysize, &ns);
if (ns == 0)
   return (0);

xcolavg = (float *) osmmget( olimit * sizeof( float ));   
ycolavg = (float *) osmmget( olimit * sizeof( float ));   
xrowavg = (float *) osmmget( olimit * sizeof( float ));   
yrowavg = (float *) osmmget( olimit * sizeof( float ));   
nrowavg = (int *) osmmget( olimit * sizeof( int ));   
ncolavg = (int *) osmmget( olimit * sizeof( int ));   

for (i = 0; i < olimit; i++) *(xcolavg+i) = 0.0;
for (i = 0; i < olimit; i++) *(ycolavg+i) = 0.0;
for (i = 0; i < olimit; i++) *(xrowavg+i) = 0.0;
for (i = 0; i < olimit; i++) *(yrowavg+i) = 0.0;
for (i = 0; i < olimit; i++) *(nrowavg+i) = 0.0;
for (i = 0; i < olimit; i++) *(ncolavg+i) = 0.0;

/*
 Compute the row or column sums.
 */
if (strcmp(order,MO_COLUMN) == 0)
   {
   for (i = 0; i < nxsub; i++) 
      {
      for (j = 0; j < nysub; j++)
          {
          if (nrshift[i][j] > 0) 
             {
             *(xrowavg+i) = *(xrowavg+i) + fabs(xrshift[i][j]);
             *(yrowavg+i) = *(yrowavg+i) + fabs(yrshift[i][j]);
             *(nrowavg+i) = *(nrowavg+i) + 1;
             } 
                    
          if (ncshift[i][j] > 0) 
             {
             *(xcolavg+i) = *(xcolavg+i) + fabs(xcshift[i][j]);
             *(ycolavg+i) = *(ycolavg+i) + fabs(ycshift[i][j]);
             *(ncolavg+i) = *(ncolavg+i) + 1;
             }
          }
      }
   } 

else 
   {
   for (i = 0; i < nysub; i++) 
      {
      for (j = 0; j < nxsub; j++)
          {
          if (nrshift[j][i] > 0) 
             {
             *(xrowavg+i) = *(xrowavg+i) + fabs(xrshift[j][i]);
             *(yrowavg+i) = *(yrowavg+i) + fabs(yrshift[j][i]);
             *(nrowavg+i) = *(nrowavg+i) + 1;
             } 
                   
          if (ncshift[j][i] > 0) 
             {
             *(xcolavg+i) = *(xcolavg+i) + fabs(xcshift[j][i]);
             *(ycolavg+i) = *(ycolavg+i) + fabs(ycshift[j][i]);
             *(ncolavg+i) = *(ncolavg+i) + 1;
             }
          }
      }
   }
        
/*
 Compute the averages.
 */
for (i = 0; i < olimit; i++)
   {
   if (*(nrowavg+i) > 0) 
      {
      *(xrowavg+i) = *(xrowavg+i) / *(nrowavg+i);
      *(yrowavg+i) = *(yrowavg+i) / *(nrowavg+i);
      }
   if ( (*ncolavg+i) > 0) 
      { 
      *(xcolavg+i) = *(xcolavg+i) / *(ncolavg+i);
      *(ycolavg+i) = *(ycolavg+i) / *(ncolavg+i);
      }
   }

/*
 Compute the medians of the row and column averages.
 */
MO_MEDR(xrowavg, nrowavg, olimit, &xrmed);
MO_MEDR(yrowavg, nrowavg, olimit, &yrmed);
MO_MEDR(xcolavg, ncolavg, olimit, &xcmed);
MO_MEDR(ycolavg, ncolavg, olimit, &ycmed);

/*
 Use the average shifts for subrasters with no information.
 */
for (j = 1; j <= nysub; j++)
   {
   if (j == nyrsub)
      jsign = 0.0;
   else if (j < nyrsub)
      jsign = 1.0;
   else
      jsign = -1.0;

   for (i = 1; i <= nxsub; i++)
      {
      if (i == nxrsub)
         isign = 0.0;
      else if (i < nxrsub)
         isign = 1.0;
      else
         isign = -1.0;

      if (nrshift[i-1][j-1] <= 0) 
         {
         if ( *(nrowavg+i-1) <= 0) 
            {
            xrshift[i-1][j-1] = isign * xrmed;
            yrshift[i-1][j-1] = jsign * yrmed;
            } 
         else if (strcmp(order,MO_COLUMN) == 0)
            {
            xrshift[i-1][j-1] = isign * *(xrowavg+i-1);
            yrshift[i-1][j-1] = jsign * *(yrowavg+i-1);
            } 
         else 
            {
            xrshift[i-1][j-1] = isign * *(xrowavg+j-1);
            yrshift[i-1][j-1] = jsign * *(yrowavg+j-1);
            }
         }

      if (ncshift[i-1][j-1] <= 0) 
         {
         if ( *(ncolavg+i-1) <= 0) 
            {
            xcshift[i-1][j-1] = isign * xcmed;
            ycshift[i-1][j-1] = jsign * ycmed;
            } 
         else if (strcmp(order,MO_COLUMN) == 0)
            {
            xcshift[i-1][j-1] = isign * *(xcolavg+i-1);
            ycshift[i-1][j-1] = jsign * *(ycolavg+i-1);
            } 
         else 
            {
            xcshift[i-1][j-1] = isign * *(xcolavg+j-1);
            ycshift[i-1][j-1] = jsign * *(ycolavg+j-1);
            }
         }
      }
   }

*nshifts = ns;
osmmfree((char *) xcolavg);
osmmfree((char *) xrowavg);
osmmfree((char *) ycolavg);
osmmfree((char *) xrowavg);
osmmfree((char *) nrowavg);
osmmfree((char *) ncolavg);
}


/*
 MO_CLINKS -- Routine to compute the shift for each subframe. 
              Input is given by a fixed shift in x and y
*/

MO_CLINKS(xrshift, yrshift, xcshift, ycshift, nxrsub, nyrsub,
          nxsub, nysub, xshift, yshift, nshifts)

float    (*xrshift)[MAXFRM];
float    (*yrshift)[MAXFRM];
float    (*xcshift)[MAXFRM];
float    (*ycshift)[MAXFRM];
int      nxrsub;
int      nyrsub;
int      nxsub;
int      nysub;
float    xshift;
float    yshift;
int      *nshifts;

{
int      i, j;
int      isign, jsign;

for (j = 0; j < nysub; j++)
    {
    if (j+1 == nyrsub)
       jsign = 0;
    else if (j+1 < nyrsub)
       jsign = 1;
    else
       jsign = -1;

    for (i = 0; i < nxsub; i++)
        {
        if (i+1 == nxrsub)
           isign = 0;
        else if (i+1 < nxrsub)
           isign = 1;
        else
           isign = -1;

        xrshift[i][j] = isign * fabs (xshift);
        yrshift[i][j] = 0.0;
        xcshift[i][j] = 0.0; 
        ycshift[i][j] = jsign * fabs (yshift);
        }
    }
*nshifts = 1;
}

/*
 MO_FLINKS -- Routine to fetch the shifts directly from the table
 */

MO_FLINKS(tblco,col,deltax, deltay, deltai, max_nshifts, nshifts)
int     tblco;                 /* table containing the shifts */
int     *col;                  /* columns containing the shifts */
float   *deltax;               /* x shifts */
float   *deltay;               /* y shifts */
float   *deltai;               /* intensity shifts */
int     max_nshifts;           /* maximum number of shifts */
int     *nshifts;

{
int     i;
int     ns;
int     ncolco;
int     nrowco;
int     nsortco;
int     allcolco;
int     allrowco;
float   data[3];
int     null;

(void) TCIGET(tblco, &ncolco, &nrowco, &nsortco, &allcolco, &allrowco);
ns = 0; 

for (i = 1; i <= nrowco; i++)
   {
   if (ns >= max_nshifts) 
      break;
   if (col[0] > 0 && col[1] > 0)
      {
      TCERDR(tblco,i,col[0],&data[0],&null);
      TCERDR(tblco,i,col[1],&data[1],&null);
      deltax[ns] = data[0];
      deltay[ns] = data[1];
      }
   else
      continue;

   if (col[2] == -1)
      {
      deltai[ns] = 0.0;
      }
   else
      TCERDR(tblco,i,col[2],&data[2],&null);
      deltai[ns] = data[2];
   ns++;
   }
*nshifts = ns;
}


/*
 MO_GET_SHIFTS -- Procedure to accumulate shifts for each subraster.
 */

MO_GET_SHIFTS(tblco, col, xrshift, yrshift, nrshift, xcshift,
        ycshift, ncshift, nxsub, nysub, nxrsub, nyrsub, nxoverlap,
        nyoverlap, nxsize, nysize, nshift)

int     tblco;
int     col[3];
float   (*xrshift)[MAXFRM];
float   (*yrshift)[MAXFRM];
int     (*nrshift)[MAXFRM];
float   (*xcshift)[MAXFRM];
float   (*ycshift)[MAXFRM];
int     (*ncshift)[MAXFRM];     /* number of column shifts */
int     nxsub;                  /* number of subrasters in x */
int     nysub;                  /* number of subrasters in y */
int     nxrsub;                 /* column index of reference subraster */
int     nyrsub;                 /* column index of reference subraster */
int     nxoverlap;              /* number of columns of overlap */
int     nyoverlap;              /* number of rows of overlap */
int     nxsize;                 /* number of columns per subraster */
int     nysize;                 /* number of rows per subraster */
int     *nshift;

{
int     allcolco, allrowco;
int     i, irow;
int     j;
int     ncolco, nrowco, nsortco;
int     ns;
int     nx1, ny1, nx11, ny11;
int     nx2, ny2, nx21, ny21;
int     r21, r22;
int     stat, stat1, stat2;
float   x1, y1, x2, y2;
float   xdif, xdifm, ydif, ydifm;
int     null;

ns = 0;
(void) TCIGET(tblco, &ncolco, &nrowco, &nsortco, &allcolco, &allrowco);
irow = 1;

do {
   if (col[0] == -1 || col[1] == -1)
      continue;
   stat1 = TCERDR(tblco,irow,col[0],&x1,&null);
   stat2 = TCERDR(tblco,irow,col[1],&y1,&null);
   if (stat1 !=0 || stat2 != 0) 
      continue;
   
/*
 Compute which subraster 1 belongs to.
 */
   if (fmod( (int) x1, nxsize) == 0)
      nx1 = (int) x1 / nxsize;
   else
      nx1 = (int) x1 / nxsize + 1;

   if (fmod( (int) y1, nysize) == 0)
      ny1 = (int) y1 / nysize;
   else
      ny1 = (int) y1 / nysize + 1;

/*
 Get the second coordinate pair.
 */
   irow++;
   do {
      if (irow > nrowco) break;
      stat1 = TCERDR(tblco,irow,col[0],&x2,&null);
      stat2 = TCERDR(tblco,irow,col[1],&y2,&null);

/*
 Compute which subraster 2 belongs to.
 */
      if (stat1 == 0 & stat2 == 0)
         {
         if (fmod( (int) x2, nxsize) == 0)
            nx2 = (int) x2 / nxsize;
         else
            nx2 = (int) x2 / nxsize + 1;

         if (fmod( (int) y2, nysize) == 0)
            ny2 = (int) y2 / nysize;
         else
            ny2 = (int) y2 / nysize + 1;
         }
      }
   while (stat1 != 0 || stat2 != 0);


   if ((irow > nrowco) || (stat1 != 0 || stat2 != 0))
      break;

   r21 = (nx1 - nxrsub) * (nx1 - nxrsub) + (ny1 - nyrsub) * (ny1 - nyrsub);
   r22 = (nx2 - nxrsub) * (nx2 - nxrsub) + (ny2 - nyrsub) * (ny2 - nyrsub);

/*
 Compute the shift for the first subraster.
 */
   if (r21 == r22) 
      continue;

   else if (r21 > r22) 
      {
      xdif = x2 - x1;
      if (nxoverlap < 0) 
         {
         if (xdif < 0.0)
            xdifm = xdif - nxoverlap;
         else if (xdif > 0.0)
            xdifm = xdif + nxoverlap;
         } 
      else
         xdifm = xdif;

      ydif = y2 - y1;
      if (nyoverlap < 0) 
         {
         if (ydif < 0.0)
            ydifm = ydif - nyoverlap;
         else if (ydif > 0.0)
            ydifm = ydif + nyoverlap;
         } 
      else
         ydifm = ydif;

      nx11 = nx1 - 1;
      ny11 = ny1 - 1;
      if (nx1 == nx2) 
         {
         xcshift[nx11][ny11] = xcshift[nx11][ny11] + xdif;
         ycshift[nx11][ny11] = ycshift[nx11][ny11] + ydifm;
         ncshift[nx11][ny11] = ncshift[nx11][ny11] + 1;
         } 
      else if (ny1 == ny2) 
         {
         xrshift[nx11][ny11] = xrshift[nx1][ny11] + xdifm;
         yrshift[nx11][ny11] = yrshift[nx1][ny11] + ydif;
         nrshift[nx11][ny11] = nrshift[nx1][ny11] + 1;
         } 
      else
         continue;

      ns++;
      }

/*
 Compute the shift for the second subraster.
 */
   else 
      {
      xdif = x1 - x2;
      if (nxoverlap < 0) 
         {
         if (xdif < 0.0)
            xdifm = xdif - nxoverlap;
         else if (xdif > 0.0)
            xdifm = xdif + nxoverlap;
         } 
      else
         xdifm = xdif;

      ydif = y1 - y2;
      if (nyoverlap < 0) 
         {
         if (ydif < 0.0)
            ydifm = ydif - nyoverlap;
         else if (ydif > 0.0)
            ydifm = ydif + nyoverlap;
         } 
      else
         ydifm = ydif;


      nx21 = nx2 - 1;
      ny21 = ny2 - 1;
      if (nx1 == nx2) 
         {
         xcshift[nx21][ny21] = xcshift[nx21][ny21] + xdif;
         ycshift[nx21][ny21] = ycshift[nx21][ny21] + ydifm;
         ncshift[nx21][ny21] = ncshift[nx21][ny21] + 1;
         } 
    
      else if (ny1 == ny2) 
         {
         xrshift[nx21][ny21] = xrshift[nx21][ny21] + xdifm;
         yrshift[nx21][ny21] = yrshift[nx21][ny21] + ydif;
         nrshift[nx21][ny21] = nrshift[nx21][ny21] + 1;
         } 
      else
         continue;
      }
      irow++;
      ns++;
   }
while (irow <=  nrowco);

/*
 Compute the final shifts.
 */
for (j = 0; j < nysub; j++)
   {
   for (i = 0; i < nxsub; i++)
      {
      if (nrshift[i][j] > 0) 
         {
         xrshift[i][j] = xrshift[i][j] / nrshift[i][j];
         yrshift[i][j] = yrshift[i][j] / nrshift[i][j];
         }
 
      if (ncshift[i][j] > 0) 
         {
         xcshift[i][j] = xcshift[i][j] / ncshift[i][j];
         ycshift[i][j] = ycshift[i][j] / ncshift[i][j];
         }
      }
   }
*nshift = ns;
}