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

/* @(#)ccdmosaic.c	17.1.1.1 (ESO-DAG) 01/25/02 17:49:36 */
/* 
++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENTIFICATION: program ccdmosaic 
.KEYWORDS:       bulk data frames, mosaicing
.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:      Rather tricky rearrangement of data array and pointers.
                 The actual algorithm was taken from the IRAF ir_mosaic code.
.INPUT/OUTPUT:   the following keys are used:
.VERSION:        930311 RHW Created; original version from Klaus (G.) Banse
.VERSION:        940118 RHW long int -> int, unit -> *int
.VERSION:        940221 RHW Modification for mosacing; original ccdcomb.c    
-------------------------------------------------- 
*/
/*
 * 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 some macros and constants
 */
#define NCOL       11
#define MAXIMS     80
#define MAXIMSONE  MAXIMS+1
#define MAXLEV     50                  /* maximum num of contour levels */
#define MAXPIX     512    /* max frame dimension (X,Y) accessed at once */
#define MAXSIZ     (MAXPIX * MAXPIX)

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

main()

{
int          uni;
int          imnoc, imnox, imnos, imnol[MAXFRM];
int          tid, col[NCOL];
int          nnull;
int          null_input[MAXFRM];
int          sizec;
int          naxis[MAXFRM], naxisc;
int          npx[2], npix[MAXFRM][3], npixc[3];
int          index[MAXFRM], isnull[MAXFRM];
int          snaxis, snpix[3];

char         *cbuff, input[72], subtr[11];

int          actvals;
int          allcol, allrow;
int          begin;
int          c1[MAXFRM], c2[MAXFRM];
int          cc;
int          chunkc;
int          colnum;
int          exist;
int          frmcnt;
int          iaux, inp, ind, indi, indo, i, ii;
int          j, jj, jjj;
int          imsize[2], nimcol, nimrow;
int          ll, lr;
int          rr, rr1, rr2;
int          iseq, iav, icu, icl, isub, ifirst, ilast, irow;
int          ijunk;
int          m, maxpix[3];
int          noutcols, noutrows;
int          nover[2];
int          noelem, nout, nolin, n, nd, nr, n1, n2, nndf;
int          nxsub, nysub;
int          nrcol, nrrow, nsort, ndum[2];
int          nsub[2];
int          nulo;
int          oline;
int          r1[MAXFRM], r2[MAXFRM];
int          stat;
int          sublo[2], subhi[2];
int          ul, ur;
int          verbose;
int          ismoot = 0;
int          inull;

float        image[4];
float        *p_img, *o_img, *pntrc;
float        usrnul, eps[3], dif, median[MAXFRM];
float        rnull;
 
double       step[MAXFRM][3], start[MAXFRM][3];
double       stepc[3], startc[3], endc[3];
double       aostep;
double       ostart[3], oldend[3], ostep[3];
double       sstart[3], sstep[3];
double       dnull;
 
char         *wpntr;
char         catfil[61];
char         cunit[61];
char         frame[MAXFRM][61], framec[61], framesig[61], framecnt[61];
char         tabnam[61];
char         ident[73];
char         line[81];
char         med_area[40], im_area[41], otime[20];
char         defaul[5], null[2];
char         output[64];

static float minundef  = -99999;
static float maxundef  =  99999;
static int   tblsiz[NCOL] = {20, 20, 1, 1, 1, 1, 1, 1, 1, 1, 1};

static char  error1[]    = "*** FATAL: operands do not match in stepsize... ";
static char  error2[]    = "*** FATAL: operands do not overlap... ";
static char  error3[]    = "*** FATAL: step sizes of different sign... ";
static char  error4[]    = "*** FATAL: catalogue empty... ";
static char  error5[]    = "*** FATAL: Unknown direction option";
static char  error6[]    = "*** FATAL: Wrong exposure type in catalogue";
static char  error7[]    = "*** FATAL: error detected in USRINP";

static char tbluni[][17] = { " ", " ", " ", " ", 
                             " ", " ", " ", " ",
                             " ", " "," "}, 

            tbllab[][17] = { "In_frame",  "im_area", 
                             "Xstartpix", "Xendpix", 
                             "Ystartpix", "Yendpix", 
                             "Median",    "Skycorr",
                             "X_offpix",   "Y_offpix",
                             "Offset"},

            tblfmt[][17] = { "A15",       "A15", 
                             "I6",        "I6",      
                             "I6",        "I6", 
                             "G12.6",     "G12.6",
                             "G12.6",     "G12.6",
                             "G12.6"};
 
static int  tbltyp[NCOL] = { D_C_FORMAT,  D_C_FORMAT,   
                             D_I2_FORMAT, D_I2_FORMAT, 
                             D_I2_FORMAT, D_I2_FORMAT, 
                             D_R4_FORMAT, D_R4_FORMAT,
                             D_R4_FORMAT, D_R4_FORMAT,			       
                             D_R4_FORMAT};  
/* 
 set up MIDAS environment + enable automatic error abort 
 */

SCSPRO("mosaic");

/*
 Get the table null value
 */
TCMNUL(&inull, &rnull, &dnull);
MO_NULL = rnull;

for (n=0; n<3; n++)
   {
   startc[n] = 0.0;
   stepc[n] = 1.0;
   endc[n] = 0.0;  
   npixc[n] = 1;
   ostart[n] = 0.0;
   oldend[n] = 0.;
   ostep[n] = 1.0;
   npix[0][n] = 1;
   start[0][n] = 0.0;
   step[0][n] = 1.0;
   }

/* 
 Get input frame list, table and output frame
 --------------------------------------------
 */
stat = SCKGETC("CCDIN",1,80,&iav,line);
stat = SCKGETC("CCDTBL",1,60,&iav,tabnam);           /* database table input */
tabnam[61] = '\0';
stat = SCKGETC("CCDOUT",1,60,&iav,framec);             /* final output file */
framec[61] = '\0';

stat = SCKGETC( "NULL_IN", 1, 72, &actvals, input );
if (strlen(input) != (size_t) 0)   
   {
   cbuff  = (char *) null_input;
   if ( USRINP( 'i', input, MAXLEV, cbuff, &nnull) != ERR_NORMAL )
      SCETER( 1, error7 );
   sorti( nnull, null_input );                    /* sort missing subrasters */
   }
else
   nnull = 0;

stat = SCKGETC("SUBTR",1,10,&iav,output);          /* get subtraction option */
CGN_UPCOPY(subtr,output,10);	                   /* upper case -> */
if (subtr[0] == 'Y') 
   isub = 1;
else
   isub = 0; 

stat = SCKGETC("VERBOSE",1,3,&iav,output);         /* Get the verbose option */
CGN_UPSTR(output);		                    /* convert to upper case */
output[4] = '\0';
if (strcmp(output,"YES") == 0)		                         /* set flag */
   {
   verbose = 1;
   strcpy(MO_DEFAULT,"NYFXN");
   }
else
   {
   verbose = 0;
   strcpy(MO_DEFAULT,"NYFNN");
   }

/* 
 Get the mosaicing parameters
 */
stat   = SCKRDI("OSIZE",1,2,&iav,imsize,&uni,&nulo);    /* size of output */
nimcol = imsize[0];
nimrow = imsize[1];

stat   = SCKRDI("NSUB",1,2,&iav,nsub,&uni,&nulo);    /* number of subimages */
MO_NXSUB = nsub[0];
MO_NYSUB = nsub[1];

stat   = SCKGETC("CORNER",1,3,&iav,output);
CGN_UPSTR(output);   		                    /* convert to upper case */
if  (strcmp(output,MO_LL) == 0 )
    MO_CORNER = 1;
else if (strcmp(output,MO_LR) == 0 )
    MO_CORNER = 2;
else if (strcmp(output,MO_UL) == 0 )
    MO_CORNER = 3;
else if (strcmp(output,MO_UR) == 0 )
    MO_CORNER = 4;
else
   SCETER(4,"*** FATAL: Unknown corner identification");         

stat = SCKGETC("DIRECT",1,3,&iav,output);                      /* direction */
CGN_UPSTR(output);		                    /* convert to upper case */
if (output[0] == 'R')
   strcpy(MO_ORDER,"ROW");
else if (output[0] == 'C')
   strcpy(MO_ORDER,"COLUMN");
else
   SCETER(4,error5);                             /* unknown combining option */

stat = SCKGETC("RASTER",1,3,&iav,output);                         /* raster */
CGN_UPSTR(output);		                    /* convert to upper case */
if (output[0] == 'Y')
   strcpy(MO_RASTER,"YES");
else 
   strcpy(MO_RASTER,"NO");

stat  = SCKRDI("NOVER",1,2,&iav,nover,&uni,&nulo);           /* overlap area */
MO_NXOVERLAP = nover[0];
MO_NYOVERLAP = nover[1];

/* get Null value */
stat = SCKGETC("BLANK",1,20,&iav,output);
if ((output[0] == '+') && (output[1] == '\0'))
   iaux = 1;                                    /*  use `last' value as Null */
else
  {
  iav = CGN_CNVT(output,2,1,npixc,&usrnul,&aostep);
  if (iav < 1)
     SCETER(19,"*** FATAL: Invalid Null value ... ");
  MO_BLANK = usrnul;
  }

/*
 Check the input catalogue for valid contents and get all inputs
 ---------------------------------------------------------------
 */
exist = 0;
inp   = 0;                                        /* default, list of frames */
icl = CGN_INDEXS(line,".cat");
icu = CGN_INDEXS(line,".CAT");
if ((icl > 0) || (icu > 0)) inp = 1;                  /* input was catalogue */
icl = CGN_INDEXS(line,".tbl");
icu = CGN_INDEXS(line,".TBL");
if ((icl > 0) || (icu > 0)) inp = 2;                      /* input was table */

/* handle input from a catalog - get name of first image file in catalog */
if (inp == 1)                                             /* catalogue input */
   {
   if ((int) strlen(line) > 60)
      SCETER(3,"*** INFO: Catalogue name too long...");
   else
      strcpy(catfil,line);

   iseq=0;
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)             /* max. MAXIMS frames */
      {
      stat = SCCGET(catfil,0,frame[frmcnt],output,&iseq);
      if (frame[frmcnt][0] == ' ') break;          /*  indicates the end ... */
      stat = SCFOPN(frame[frmcnt],D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
      }
   if (frmcnt <= 0) SCETER(4,error4);     /* there must be at least 1 image  */
   }

else                                              /* input was frame string */
   {
   begin = 0;
   m = (int) strlen(line);
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)
      {
      ll = CGN_EXTRSS(line,m,',',&begin,output,60);
      if (ll <= 0) break;
       
      CGN_FRAME(output,(int)F_IMA_TYPE,frame[frmcnt],0);  /* convert frames */
      stat = SCFOPN(frame[frmcnt],D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
      }
   if (frmcnt <= 0) SCETER(4,error4);     /* there must be at least 1 image  */
   }


if (MO_NXSUB * MO_NYSUB != frmcnt + nnull)
   SCETER(3,
          "*** FATAL: Number of input frames not equal to nxsub * nysub");

/*
 Get initial area from 1. frame and check the start and descriptors 
 */
stat = SCDRDI(imnol[0],"NAXIS",1,1,&iav,&naxis[0],&uni,&nulo);
stat = SCDRDI(imnol[0],"NPIX",1,3,&iav,&npix[0][0],&uni,&nulo);
stat = SCDRDD(imnol[0],"START",1,3,&iav,&start[0][0],&uni,&nulo);
stat = SCDRDD(imnol[0],"STEP",1,3,&iav,&step[0][0],&uni,&nulo);

for (m=0; m<3; m++)
   {
   ostart[m] = start[0][m];
   ostep[m] = step[0][m];
   aostep = ostep[m];
   if (aostep < 0.0) aostep = - aostep;
   eps[m] = 0.0001 * aostep;	   		   /* take 0.01% of stepsize */
   oldend[m] = ostart[m] + (npix[0][m]-1)*ostep[m];
   maxpix[m] = npix[0][m];
   }

/* 
 Now loop through the other input frames  
 */
for (n=0; n<frmcnt; n++)
   {
   for (m=0; m<3; m++)                                     /* init values... */
      {
      start[n][m] = 0.0;
      step[n][m] = 1.0;
      npix[n][m] = 1;
      }
  
   stat = SCDRDI(imnol[n],"NAXIS",1,1,&iav,&naxis[n],&uni,&nulo);
   stat = SCDRDI(imnol[n],"NPIX",1,3,&iav,&npix[n][0],&uni,&nulo);
   stat = SCDRDD(imnol[n],"START",1,3,&iav,&start[n][0],&uni,&nulo);
   stat = SCDRDD(imnol[n],"STEP",1,3,&iav,&step[n][0],&uni,&nulo);

/* 
 Stepsizes should have same sign and not differ too much...  
 */
   for (m=0; m<3; m++)
      {
      if ((ostep[m]*step[n][m]) <= 0.) SCETER(1,error3);
      aostep = step[n][m] - ostep[m];
      if (aostep < 0.) aostep = -aostep;
      if (aostep > eps[m]) SCETER(5,error1);

      dif = start[n][m] + (npix[n][m]-1)*step[n][m];
      if (ostep[m] < 0.) 
         {
         if (ostart[m] < start[n][m]) ostart[m] = start[n][m];     /* MAX */
         if (oldend[m] > dif) oldend[m] = dif;		           /* MIN */
         }
      else
         {
         if (ostart[m] > start[n][m]) ostart[m] = start[n][m];     /* MIN */
         if (oldend[m] < dif) oldend[m] = dif;		           /* MAX */
         }
      maxpix[m] = MYMAX( npix[n][m], maxpix[m] );
      }
   }

for (m=0; m<3; m++)                     /* test, if something is left...  */
   if (ostep[m]*(oldend[m]-ostart[m]) < 0.) SCETER(2,error2);

stat = SCKGETC("IM_SEC",1,40,&iav,im_area);           /* section to extract */

stat = SCKGETC("SECTION",1,40,&iav,med_area);      /* section for the median */

/* 
 Check for correct dimensionality; we only work on max 2-dim. frames...   
 */
naxisc = naxis[0];
for (n=1; n<frmcnt; n++)
   if (naxisc < naxis[n]) naxisc = naxis[n];		  	  /* MAX */

if (naxisc > 2) 
   {
   SCTPUT("*** INFO: Currently only 1-dim and 2-dim frames supported...");
   naxisc = 2;
   }

/* 
 get the image area from the first image
 */
npx[0]  = npix[0][0];			/* still needed? */
npx[1]  = npix[0][1];
stat = Convcoo(1,imnol[0],im_area,2,ndum,sublo,subhi);
if ( stat != ERR_NORMAL ) 
   SCETER( 2, "*** FATAL: Problems with input coordinates" );
image[0] = sublo[0] + 1;
image[1] = subhi[0] + 1;
image[2] = sublo[1] + 1;
image[3] = subhi[1] + 1;
nxsub    = (int) ( image[1] - image[0] ) + 1;
nysub    = (int)( image[3] - image[2] ) + 1; 

/* 
 get column and rows for first image
 */
MO_NCOLS = npix[0][0];
MO_NROWS = npix[0][1];

/*
 Compute the size of the output frame 
 */
ijunk = MO_NXSUB * maxpix[0] - (MO_NXSUB - 1) * MO_NXOVERLAP;
if (nimcol == 0)
   npixc[0] = ijunk;
else
   npixc[0] = MYMAX(nimcol,ijunk);

ijunk = MO_NYSUB * maxpix[1] - (MO_NYSUB - 1) * MO_NYOVERLAP;
if (nimrow == 0)
   npixc[1] = ijunk;
else
   npixc[1] = MYMAX(nimrow,ijunk);

/* 
 Set up standard stuff for result frame framec 
 */
sizec = 1;
for (m=0; m<naxisc; m++)
   {
   startc[m] = ostart[m];
   stepc[m]  = ostep[m];
   sizec     = sizec * npixc[m];
   }
stat = SCFCRE(framec,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnoc);
stat = SCDWRI(imnoc,"NAXIS",&naxisc,1,1,&uni);
stat = SCDWRI(imnoc,"NPIX",npixc,1,naxisc,&uni);
stat = SCDWRD(imnoc,"START",startc,1,naxisc,&uni);
stat = SCDWRD(imnoc,"STEP",stepc,1,naxisc,&uni);
stat = SCDCOP(imnol[0], imnoc, 4, "CUNIT");
sprintf(ident,"Mosaic of input catalogue %s", catfil);
stat = SCDWRC(imnoc,"IDENT",1,ident,1,72,&uni);

/*
  Set up the data base table for the frames
 */
(void) MO_SETUP(imnol, naxis, npix, start, step, nnull,
                null_input, im_area, med_area, npixc, 
                index, c1, c2, r1, r2, isnull, median);

/*
 Make the mosaic output image
 */
sortii(MO_NXSUB * MO_NYSUB, r1, index);              /* sort index array */

noutcols = npixc[0];
noutrows = npixc[1];
ifirst   = 1;
oline    = 1;

/* 
 Here we do the real work, that is reading and writing the data
 First, allocate size of one full strip with noutcols
 */
SCTPUT(" ");
sprintf(line,"Input catalogue=%s, Output frame=%s",catfil,framec);
SCTPUT(line);
sprintf(line,"Number of input subrasters= %d; number of empty subraster= %d",
        frmcnt, nnull);
SCTPUT(line);
if (isub == 1 ) 
   SCTPUT("Subrasters corrected for median value");
else
   SCTPUT("Subraster not corrected for median value");

chunkc = npixc[0] * npixc[1];                      /* size of 1 result strip */
ll = chunkc * sizeof(float);                      /* here for the output map */
wpntr = osmmget((unsigned int)ll);
if (wpntr == (char *) 0)                         /* get pointer output frame */
   SCETER(66,"*** FATAL: Could not allocate virtual memory ...");  
else
   pntrc = (float *) wpntr;

p_img = (float *) osmmget( npixc[0] * npixc[1] * sizeof( float )); 

/* 
 Now loop over all subframes
 */
for ( i = 0; i < MO_NXSUB * MO_NYSUB; i += MO_NXSUB )   /* loop over y range */
   {
   rr1 = r1[index[i]];
   rr2 = r2[index[i]];

   for (j = i; j < i + MO_NXSUB; j++)                   /* loop over x range */
      {
      ind = j * nxsub * nysub;
      if (isnull[index[j]] < 0) 
         for (ii = 0; ii < nxsub; ii++)
            *(p_img + ind + ii) = MO_BLANK;
      else
         GETDAT(imnol[isnull[index[j]]], MAXSIZ, npx, image, 
                ismoot, &p_img[ind]);
      }

/*
 open the nxsub input frames
 */
   while (oline < rr1)                              /* first the blank lines */
      {
      for (j = 0; j <noutcols; j++)
         *(pntrc + (oline-1) * noutcols + j) = MO_BLANK;
      oline++;
      }

/* 
 write the output frame with the data
 */   
   irow = 0;
   for (rr = rr1; rr <= rr2; rr++)
      {
      ind = (rr - 1) * noutcols;                   /* offset in output frame */
      for (j = 0; j < noutcols; j++)                    /* blank lines first */
         *(pntrc + ind + j) = MO_BLANK;

      for (jj = 0; jj < MO_NXSUB; jj++)         /* loop through nx subframes */
         {
         jjj  = index[jj + i];                          /* index of subframe */
         indi = jjj * nxsub * nysub;                /* offset in p_img array */
         if ( *(p_img + indi) != MO_BLANK )                   /* not a blank */
	    {
            indi = indi + irow * nxsub;                    /* offset for row */
            for (cc = 0; cc < (c2[jj] - c1[jj] + 1); cc++)
               { 
               indo =  ind + c1[jj] + cc - 1;
               if (isub == 1)
                  *(pntrc + indo) = *(p_img + indi + cc) - median[jj];
               else
                  *(pntrc + indo) = *(p_img + indi + cc);
	       }
            }
         }
      irow++;   
      }
   oline += rr2;
   }
 
while (oline < noutrows)
   {
   for (ii = 0; ii < noutcols; ii++)
      *(pntrc + (oline-1) * noutcols + ii) = MO_BLANK;
   oline++;
   }            

/* 
 store the data on disk
 */
(void) SCFPUT(imnoc, 1, npixc[0]*npixc[1], (char *) pntrc);

/*
 Write the mosaicing parematers to the dat file
 */ 
(void) TCTINI( tabnam, F_TRANS, F_O_MODE, NCOL, frmcnt+nnull, &tid );
for ( ii = 0; ii < NCOL; ii++ )
    (void) TCCINI( tid, tbltyp[ii], tblsiz[ii], tblfmt[ii] , tbluni[ii],
                         tbllab[ii], &col[ii] );
(void) MO_TBLWPAR(tid, im_area, med_area);
(void) MO_TBLWINP(tid, catfil,frmcnt,framec,im_area,index,c1,c2,r1,r2,
                  isnull,median,isub,verbose);
(void) TCTCLO(tid);

/* 
 Finally, close and return: bye, bye, zwaai, zwaai
 */
return SCSEPI();
}