/* @(#)mo_align.c	17.1.1.1 (ES0-DMD) 01/25/02 17:49:05 */
/*===========================================================================
  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_align.c
.AUTHOR      R.H. Warmels IPG-ESO Garching
.KEYWORDS    alignment software
.LANGUAGE    C
.PURPOSE     Routine to align all subrasters
.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_SUBALIGN -- Align all subrasters
 */

MO_SUBALIGN(inim, inpntr, outim, outpntr, area, 
            ic1, ic2, il1, il2, oc1, oc2, ol1, ol2,
            deltax, deltay, deltai, match, interp, verb)

int     inim;
char    *inpntr;
int     outim;
char    *outpntr;    
int     *area;
int     *ic1;
int     *ic2;
int     *il1;
int     *il2;
int     *oc1;
int     *oc2;
int     *ol1;
int     *ol2;
float   *deltax;
float   *deltay;
float   *deltai;
int     match;
int     interp;
int     verb;
{

char    line[81];
int     cin1, cin2;
int	i;
int     iav;
int	ideltax, ideltay;
int     index;
int     ism;
int     ll;
int     k, kk;
int     lin1, lin2;
int     lout1, lout2;
int     lxoffset, hxoffset, lyoffset, hyoffset;
int	ixoffset, iyoffset;
int     nicols, nilines;
int     nimcols, nimlines, nimages;
int     nocols, norows;
int     npixi[2], npixo[2];
int     nyout, fstline, lstline;
int     nulo;
int     size;
int     stat;
int	tlin1;
int     uni;

float   image[4];
float   dx, dy, ytemp;
float   *x, *y, *inbuf, *outbuf, *ptr;

/*
 Get smooting option of GETDAT
 */
ism = 0;
 
/* 
 Allocate memory for arrays
 */
x = (float *) osmmget( MO_NCOLS * sizeof( float ));
y = (float *) osmmget( MO_NCOLS * sizeof( float ));

/* 
 allocate memory for input data
 */
stat     = SCDRDI(inim,"NPIX",1,3,&iav,npixi,&uni,&nulo);
nimcols  = npixi[0];
nimlines = npixi[1];
size     = npixi[0] * npixi[1];
inbuf    = (float *) osmmget( size * sizeof( float ));
stat     = SCDRDI(outim,"NPIX",1,3,&iav,npixo,&uni,&nulo);

/*
 Initialize the interpolant.
 */
MO_MSIINIT(interp, nimcols, nimlines);

/*
 Extract the subrasters one by one.
 */
nimages = MO_NXSUB * MO_NYSUB;

for (i = 0; i < nimages; i++)
   {
   *inbuf = MO_NULL;

/*
 Reject and subraster which is off the image.
 */
   if (oc1[i] > nimcols || oc2[i] < 1 || ol1[i] > nimlines || ol2[i] < 1)
      continue;

/*
 Compute the integer and fractional part of the shift.
 */
   ideltax  = NINT (deltax[i]);
   ideltay  = NINT (deltay[i]);
   dx       = deltax[i] - ideltax;
   dy       = deltay[i] - ideltay;

/*
 Compute the output image limits.
 */
   lxoffset = MYMAX (1 - oc1[i], area[0]);
   hxoffset = MYMAX (oc2[i] - nimcols, area[1]);
   oc1[i]   = MYMAX (1, MYMIN (nimcols, oc1[i] + lxoffset));
   oc2[i]   = MYMIN (nimcols, MYMAX(1, oc2[i] - hxoffset));
   nocols   = oc2[i] - oc1[i] + 1;
   lyoffset = MYMAX (1 - ol1[i], area[2]);
   hyoffset = MYMAX (ol2[i] - nimlines, area[3]);
   ol1[i]   = MYMAX (1, MYMIN (nimlines, ol1[i] + lyoffset));
   ol2[i]   = MYMIN (nimlines, MYMAX (1, ol2[i] - hyoffset));
   norows   = ol2[i] - ol1[i] + 1;

/*
 Compute some input image parameters.
 */
   cin1     = MYMAX (ic1[i], MYMIN (ic1[i] + lxoffset - MO_NMARGIN, ic2[i]));
   cin2     = MYMIN (ic2[i], MYMAX (ic2[i] - hxoffset + MO_NMARGIN, ic1[i]));
   nicols   = cin2 - cin1 + 1;
    
/* 
 Compute the x offset and x interpolation coordinates.
 */
   ixoffset = MYMIN (lxoffset, MO_NMARGIN);
   for (k = 0; k < nicols; k++)
      *(x+k) = MYMAX(1.0, MYMIN ((float) nicols, (float) (k + ixoffset - dx)));

/*
 Subdivide the image and do the shifting.
 */
   for (lout1 = ol1[i]; lout1 <= ol2[i]; lout1 = lout1 + MO_NYOUT) 
      {
/*
Compute the output image limits.
 */
      lout2 = MYMIN (ol2[i], lout1 + MO_NYOUT - 1);
      nyout = lout2 - lout1 + 1;

/*
 Compute the input image limits.
 */
      tlin1   = il1[i] + lyoffset + lout1 - ol1[i];
      lin2    = MYMIN (il2[i], MYMAX (tlin1 + nyout + MO_NMARGIN - 1, il1[i]));
      lin1    = MYMAX (il1[i], MYMIN (tlin1 - MO_NMARGIN, il2[i]));
      nilines = lin2 - lin1 + 1;

/*
 Get the appropriate input image section and fit the interpolant.
 */
      if ((*inbuf == MO_NULL) || (lin1 < fstline) || (lin2 > lstline)) 
         {
	 fstline  = lin1;
	 lstline  = lin2;
         image[0] = (float) cin1;
         image[1] = (float) cin2;
         image[2] = (float) lin1;
         image[3] = (float) lin2;
	 GETDAT (inim, size, npixi, image, ism, inbuf);
	 MO_MSIFIT(inbuf, nicols, nilines, nicols);
         }
   
/*
 Get the y offset and y interpolation coordinates. 
 */
      if (lout1 == ol1[i])
	 iyoffset = MYMIN (lyoffset, MO_NMARGIN);
      else
	 iyoffset = tlin1 - lin1;

/*
 Shift the input images.
 */
      index  = (lout1 - 1) * npixo[0] + oc1[i] - 1;
      outbuf = (float *) outpntr + index;
      ptr    = outbuf;
      for (k = 1; k <= nyout; k++)
	 {
         ytemp = MYMAX (1.0, MYMIN ((float) nilines, 
                                   (float) (k + iyoffset - dy))); 
         for (kk = 0; kk < nocols; kk++) *(y+kk) = ytemp;
  	 MO_MSIVECTOR(x, y, ptr, nocols);
	 ptr += npixo[0];
         }

/*
 Shift the intensities
 */
      if (match == 1 && deltai[i] != MO_INDEFR)
         {
         outbuf = (float *) outpntr + index;
         ptr    = outbuf;
         for (k = 0; k < nyout; k++) 
            {
            for (kk = 0; kk < nocols; kk++)
               *(ptr+kk) += deltai[i];
            ptr += npixo[0];
  	    }
         }
      }
   *inbuf = MO_NULL;

/*
 Print a message.
 */
   if (verb == 1)
      {
      sprintf(line, "Input:  frame #%3d org_area =  [%4d,%4d:%4d,%4d]",
         i, lxoffset + 1, lyoffset + 1, lxoffset + nocols, lyoffset + norows);
      SCTPUT(line); 
      
      sprintf(line, "                   mos_area = [%4d,%4d:%4d,%4d]",
           ic1[i], il1[i], ic2[i], il2[i]);
      SCTPUT(line); 

      sprintf(line, "                   Xoffset = %g Yoffset=%g",
          deltax[i], deltay[i]);
      SCTPUT(line); 

      sprintf(line, 
        "Output:            out_area = [%4d,%4d:%4d,%4d] Offset = %g", 
                    oc1[i], ol1[i], oc2[i], ol2[i], deltai[i]);
      SCTPUT(line);
      }
   }
}