/* @(#)ccdalign.c 17.1.1.1 (ES0-DMD) 01/25/02 17:49:35 */ /*=========================================================================== 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 ===========================================================================*/ /* @(#)ccdalign.c 17.1.1.1 (ESO-DAG) 01/25/02 17:49:35 */ /* ++++++++++++++++++++++++++++++++++++++++++++++++++ .IDENTIFICATION: program ccdalign .KEYWORDS: bulk data frames, mosaicing, aligment .PURPOSE: Align the individual subraster elements in the input image. In order to run this program the user should have created the output image and the database file with the IRMOSAIC task. In addition the user should supply a coordinate list consisting of pairs of coordinates of identical objects or features in two adjacent subrasters. The actual algorithm was taken from the IRAF ir package. -------------------------------------------------- */ /* * Define _POSIX_SOURCE to indicate * that this is a POSIX program */ #define _POSIX_SOURCE 1 /* * definition of the used functions */ #include #include #include #include #include #include #include #include /* * 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 imnoi, imnoc; int tdb, tco; int sizec; int naxis, naxisc; int npix[3], npixc[3]; int allcoldb, allrowdb; int colc[3], colo[3]; int iaux, i, ialign; int imsize[2], nimcol, nimrow; int iav; int ll; int nimages; int nxsub, nysub, nshift; int nrsub[2], nrxsub, nrysub; int ncoldb, nrowdb, nsortdb, ndum[2]; int nulo; int oline; int stat; int sublo[2]; int subhi[2]; int tr_area[4]; int xyref[2]; int verbose; int match; int interp; int inull; char *pntri, *pntrc; float data[3]; float usrnul; float xyshift[2]; float xshift, yshift; float rnull; double step[3], start[3]; double stepc[3], startc[3], endc[3]; double aostep; double dnull; char cunit[61]; char input[72]; char framei[61], framec[61]; char tabdb[61], tabco[61]; char ident[72]; char line[81]; char med_area[41], im_area[41]; char defaul[5], null[2]; char output[64], align[11]; static float minundef = -99999; static float maxundef = 99999; static int tblsiz[NCOL] = {20, 15, 1, 1, 1, 1, 1, 1, 1, 1, 1}; static char x_colo[] = "X_offpix"; static char y_colo[] = "Y_offpix"; static char i_colo[] = "Offset"; static char x_colc[] = "X_coordpix"; static char y_colc[] = "Y_coordpix"; static char i_colc[] = "Value"; 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 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("align"); /* Get the table null value */ TCMNUL(&inull, &rnull, &dnull); MO_NULL = 0.0; /* Get input frame list, tables and output frame -------------------------------------------- */ stat = SCKGETC("CCDIN",1,60,&iav,framei); /* input mosaic frame */ stat = SCIGET(framei,D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, MO_DIM2, &naxis, npix, start, step, ident, cunit, &pntri, &imnoi); /* Get database, open it and get info Store the table descriptor into the MO variables */ stat = SCKGETC("TBLDB",1,60,&iav,tabdb); /* database table */ (void) TCTOPN(tabdb, F_IO_MODE, &tdb); (void) TCIGET(tdb, &ncoldb, &nrowdb, &nsortdb, &allcoldb, &allrowdb); TCCSER(tdb, x_colo, &colo[0]); TCCSER(tdb, y_colo, &colo[1]); TCCSER(tdb, i_colo, &colo[2]); MO_TBLRPAR(tdb, im_area, med_area); /* read the desciptor */ nimages = MO_NXSUB * MO_NYSUB; /* Get the output file */ stat = SCKGETC("ccdout",1,60,&iav,framec); /* final output file */ /* Get the alignment keyword */ stat = SCKGETC("ALIGN",1,10,&iav,output); /* get subtraction option */ CGN_UPCOPY(align,output,10); /* upper case -> */ if (align[0] == 'C') /* here for coordinates table*/ { stat = SCKGETC("TBLCO",1,60,&iav,tabco); /* coordinate input */ ialign = 0; } else if (align[0] == 'S') /* here for shift table */ { stat = SCKRDR("SHIFT",1,2,&iav,xyshift,&uni,&nulo); /* x and y shifts */ xshift = xyshift[0]; yshift = xyshift[1]; ialign = 1; } else if (align[0] == 'R') /* here for fixed shift */ { stat = SCKGETC("TBLCO",1,60,&iav,tabco); /* coordinate input */ ialign = 2; } /* get the image section */ stat = SCKRDI("TR_SEC",1,4,&iav,tr_area,&uni,&nulo); /* section to include */ /* Column and row number of reference raster */ stat = SCKRDI("NRSUB",1,2,&iav,nrsub,&uni,&nulo); MO_NXRSUB = nrsub[0]; MO_NYRSUB = nrsub[1]; if (MO_NXRSUB == 0 || MO_NXRSUB < 1 || MO_NXRSUB > MO_NXSUB) MO_NXRSUB = (MO_NXSUB + 1) / 2; if (MO_NYRSUB == 0 || MO_NYRSUB < 1 || MO_NYRSUB > MO_NYSUB) MO_NYRSUB = (MO_NYSUB + 1) / 2; /* Get the x and y offset of the reference subraster */ stat = SCKRDI("XYREF",1,2,&iav,xyref,&uni,&nulo); MO_XREF = xyref[0]; MO_YREF = xyref[1]; /* Get the size of the output */ stat = SCKRDI("OSIZE",1,2,&iav,imsize,&uni,&nulo); /* size of output */ nimcol = imsize[0]; nimrow = imsize[1]; if ( nimcol != 0 && nimcol > 0 && nimcol >= npix[0]) npixc[0] = nimcol; else npixc[0] = npix[0]; if ( nimrow != 0 && nimrow > 0 && nimrow >= npix[1]) npixc[1] = nimrow; else npixc[1] = npix[1]; /* Get Null value for the output frame */ 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; } /* get the interpolation */ stat = SCKGETC("INTERPOL",1,40,&iav,output); /* section to include */ CGN_UPSTR(output); /* convert to upper case */ if (strncmp(output,"NEA",3) == 0) /* set flag */ interp = MO_BINEAREST; else if (strncmp(output,"LIN",3) == 0) interp = MO_BILINEAR; else if (strncmp(output,"POLY3",5) == 0) interp = MO_BIPOLY3; else if (strncmp(output,"POLY5",5) == 0) interp = MO_BIPOLY5; else if (strncmp(output,"SPLINE3",7) == 0) interp = MO_BISPLINE3; else interp = 999; /* Get the output mode */ 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"); } /* Allocate memory for output frame and zero it */ sizec = npixc[0]*npixc[1]; stat = SCFCRE(framec,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnoc); stat = SCFMAP(imnoc, F_O_MODE, 1, sizec, &iav, &pntrc); if (stat != 0) /* get pointer output frame */ SCETER(66,"*** FATAL: Could not allocate virtual memory ..."); stat = SCDWRI(imnoc,"NAXIS",&naxis,1,1,&uni); stat = SCDWRI(imnoc,"NPIX",npix,1,naxis,&uni); stat = SCDWRD(imnoc,"START",start,1,naxis,&uni); stat = SCDWRD(imnoc,"STEP",step,1,naxis,&uni); stat = SCDCOP(imnoi, imnoc, 4, "CUNIT"); sprintf(ident,"Alignment of input frame %s", framei); stat = SCDWRC(imnoc,"IDENT",1,ident,1,72,&uni); MO_ZERO(pntrc, npixc, MO_BLANK); (void) SCFPUT(imnoc, 1, npixc[0]*npixc[1], pntrc); /* 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 frame: %s", framei); SCTPUT(line); sprintf(line,"Database table: %s", tabdb); SCTPUT(line); sprintf(line,"Output frame: %s",framec); SCTPUT(line); sprintf(line,"Number of subrasters (x,y): %d,%d", MO_NXSUB, MO_NYSUB); SCTPUT(line); if (align[0] == 'C') /* here for coordinates table */ sprintf(line,"Method is Coord; Coordinate table: %s",tabco); else if (align[0] == 'S') /* here for fixed shift */ sprintf(line,"Method is Shift; xshift = %f; yshift = %f", xshift, yshift); else if (align[0] == 'R') /* here for reference shift */ sprintf(line,"Method is Reference; Shift table: %s",tabco); SCTPUT(line); /* Now we have all parameter input available; time to go to work now */ switch(ialign) { case 0: if (strcmp(tabco,tabdb) != 0) { (void) TCTOPN(tabco, F_IO_MODE, &tco); TCCSER(tco, x_colc, &colc[0]); TCCSER(tco, y_colc, &colc[1]); TCCSER(tco, i_colc, &colc[2]); } else { tco = tdb; TCCSER(tco, x_colo, &colc[0]); TCCSER(tco, y_colo, &colc[1]); TCCSER(tco, i_colo, &colc[2]); } if (colc[0] < 0 || colc[1] < 0 ) SCETER(4,"*** FATAL: X and/or Y offset column(s) not found"); MO_LINKS(tco, colc, MO_XRSHIFTS, MO_YRSHIFTS, MO_XCSHIFTS, MO_YCSHIFTS, MO_NRSHIFTS, MO_NCSHIFTS, MO_NCOLS, MO_NROWS, MO_NXRSUB, MO_NYRSUB, MO_NXSUB, MO_NYSUB, MO_NXOVERLAP, MO_NYOVERLAP, MO_ORDER, &nshift); if (nshift > 0) MO_SHIFTS(imnoi,imnoc, MO_XRSHIFTS, MO_YRSHIFTS, MO_XCSHIFTS, MO_YCSHIFTS, MO_IC1, MO_IC2, MO_IL1, MO_IL2, MO_OC1, MO_OC2, MO_OL1, MO_OL2, MO_DELTAX, MO_DELTAY); else SCETER(4,"*** FATAL: Illegal shifts in the coordinate file"); break; case 1: MO_CLINKS(MO_XRSHIFTS, MO_YRSHIFTS, MO_XCSHIFTS, MO_YCSHIFTS, MO_NXRSUB, MO_NYRSUB, MO_NXSUB, MO_NYSUB, xshift, yshift, &nshift); if (nshift > 0) MO_SHIFTS(imnoi,imnoc, MO_XRSHIFTS, MO_YRSHIFTS, MO_XCSHIFTS, MO_YCSHIFTS, MO_IC1, MO_IC2, MO_IL1, MO_IL2, MO_OC1, MO_OC2, MO_OL1, MO_OL2, MO_DELTAX, MO_DELTAY); else SCETER(4,"*** FATAL: Illegal shifts in the coordinate file"); break; case 2: if (strcmp(tabco,tabdb) != 0) { (void) TCTOPN(tabco, F_IO_MODE, &tco); TCCSER(tco, x_colo, &colc[0]); TCCSER(tco, y_colo, &colc[1]); TCCSER(tco, i_colo, &colc[2]); } else { tco = tdb; TCCSER(tco, x_colo, &colc[0]); TCCSER(tco, y_colo, &colc[1]); TCCSER(tco, i_colo, &colc[2]); } if (colc[0] < 0 || colc[1] < 0 ) SCETER(4,"*** FATAL: X and/or Y offset column(s) not found"); MO_FLINKS(tco, colc, MO_DELTAX, MO_DELTAY, MO_DELTAI, nimages, &nshift); if (nshift < nimages) SCETER(4,"*** FATAL: Fewer shift than subrasters"); else MO_FSHIFTS(imnoi, imnoc, MO_DELTAX, MO_DELTAY, MO_IC1, MO_IC2, MO_IL1, MO_IL2, MO_OC1, MO_OC2, MO_OL1, MO_OL2); break; default: SCETER(4,"*** FATAL undefined alignment algorithm"); } /* Shift all subrasters */ for (i = 0; i < nimages; i++) MO_DELTAI[i] = 0.0; match = 0; MO_SUBALIGN(imnoi, pntri, imnoc, pntrc,tr_area, MO_IC1,MO_IC2,MO_IL1,MO_IL2, MO_OC1,MO_OC2,MO_OL1,MO_OL2, MO_DELTAX,MO_DELTAY,MO_DELTAI, match,interp,verbose); /* Store the data in the data base file and close files. Finally exit */ (void) SCFPUT(imnoc, 1, npixc[0]*npixc[1], pntrc); switch(ialign) { case 0: for (i = 0; i < nrowdb; i++) { data[0] = MO_DELTAX[i]; data[1] = MO_DELTAY[i]; data[2] = MO_DELTAI[i]; stat = TCRWRR(tdb, i+1, 2, colo, data); } (void) TCTCLO(tdb); (void) TCTCLO(tco); break; case 1: for (i = 0; i < nrowdb; i++) { data[0] = MO_DELTAX[i]; data[1] = MO_DELTAY[i]; data[2] = MO_DELTAI[i]; stat = TCRWRR(tdb, i+1, 2, colo, data); } (void) TCTCLO(tdb); break; case 2: for (i = 0; i < nrowdb; i++) { data[0] = MO_DELTAX[i]; data[1] = MO_DELTAY[i]; data[2] = MO_DELTAI[i]; stat = TCRWRR(tdb, i+1, 2, colo, data); } (void) TCTCLO(tdb); if (strcmp(tabco,tabdb) != 0) (void) TCTCLO(tco); break; default: SCETER(4,"*** FATAL undefined alignment algorithm"); } return SCSEPI(); }