/* COPYRIGHT (c) 1990 Kapteyn Astronomical Institute University of Groningen, The Netherlands All Rights Reserved. #> flat.dc1 Program: FLAT Purpose: Apply a flat field correction to a map by fitting a 2-D polynomial to the background. Category: MODELS, CALCULATION, FITTING, CALCULATION File: flat.c Author: M. Vogelaar Keywords: INSET= Give set, subsets: Maximum number of subsets is 2048. Dimension of the subsets must be 2. BOX= Area which needs to be flat fielded [whole map] OUTSET= Output Set/Subset(s) ORDER= Order of polynomial to be fitted to the map [1] highest allowed order is 21 ** RANGE= Give range of levels to include in fit: [all levels] Input consists of two values. The first value may be greater than the second. See the description for the correct use of this keyword. Description: This program can be used for example to extract a background by fitting a polynomial to it. The program uses a set (INSET=) which usually is blotted first. The subsets must be two dimensional. The area on which FLAT works is given in BOX= It writes the results to OUTSET= in the same box. Beware if OUTSET= already exists. It only overwrites data in the given box. A polynomial plane will be fitted to the data using the least squares method. The order of this polynomial is given in ORDER= This number cannot exceed 21. Values outside the range defined in RANGE= are treated as blanks, i.e. they are not included in the fit and are interpolated after a plane is fitted to the data. Examples for the use of the RANGE= keyword: RANGE=2, 5 (IF 2blank) RANGE=5, 2 (IF 2blank ELSE include in fit) At the RANGE= keyword, the values -INF and INF can be input also. These values represent the minimum and maximum values of the current system. RANGE=5, INF (IF value>5 THEN include ELSE value==>blank) The fitted background can be subtracted from the original set with the program COMBIN. Example: Notes: Updates: Oct 14, 1991: MV, Document created. #< */ #include "stdio.h" #include "stdlib.h" #include "string.h" #include "math.h" #include "cmain.h" #include "gipsyc.h" #include "init.h" #include "finis.h" #include "float.h" #include "gdsinp.h" #include "gdsasn.h" #include "gdsout.h" #include "setfblank.h" #include "myname.h" #include "anyout.h" #include "nelc.h" #include "presetr.h" #include "gdsbox.h" #include "gdsc_range.h" #include "gdsc_grid.h" #include "gdsc_fill.h" #include "gdsc_name.h" #include "gdsi_read.h" #include "gdsi_write.h" #include "gdsc_ndims.h" #include "userchar.h" #include "userreal.h" #include "userint.h" #include "error.h" #include "reject.h" #include "getrange.h" #include "minmax3.h" #include "wminmax.h" #include "status.h" #define AXESMAX 10 /* Max. allowed number of axes in a set */ #define SUBSMAX 2048 /* Max. number of substructures to be specified */ #define MAXELEM 256 /* Max. # coefficients in polynomial */ #define BIGSTORE 80 /* Length of a string */ #define VERSION "1.0" /* Version number of this program */ #define NONE 0 /* Default values for use in userxxx routines */ #define REQUEST 1 #define HIDDEN 2 #define EXACT 4 #define false 0 #define true 1 #define PI 3.141592653589793 /* Keywords and messages */ #define KEY_INSET tofchar("INSET=") #define MES_INSET tofchar("Give set (, subsets): ") #define KEY_BOX tofchar("BOX=") #define MES_BOX tofchar("Area which needs to be flat fielded: [whole map]") #define KEY_OUTSET tofchar("OUTSET=") #define MES_OUTSET tofchar("Output set:") #define KEY_ORDER tofchar("ORDER=") #define MES_ORDER tofchar("Order of polynomial to be fitted: [1]") #define KEY_RANGE tofchar("RANGE=") #define MES_RANGE tofchar("Range of levels to include in fit: [all levels]") /* Initialize string with macro */ #define fmake(fchr,size) { \ static char buff[size+1]; \ int i; \ for (i = 0; i < size; buff[i++] = ' '); \ buff[i] = 0; \ fchr.a = buff; \ fchr.l = size; \ } /* Malloc version of 'fmake' */ #define finit( fc , len ) { fc.a = malloc( ( len + 1 ) * sizeof( char ) ) ; \ fc.a[ len ] = '\0' ; \ fc.l = len ; } #define MYMAX(a,b) ((a) > (b) ? (a) : (b)) #define MYMIN(a,b) ((a) > (b) ? (b) : (a)) #define NINT(a) ( (a)<0 ? (int)((a)-.5) : (int)((a)+.5) ) /* Input of set, subsets: */ static fchar Setin; /* Name of the set */ static fint Subin[SUBSMAX]; /* Array for the subset coordinate words */ static fint Nsubsin; /* Number of input subsets */ static fint Dfault; /* Default option for input etc */ static fint Axnum[AXESMAX]; /* GDSINP axis numbers array */ static fint Axcount[AXESMAX]; /* GDSINP axis lengths array */ static fint Class = 1; /* Axis is operation axis */ static fint Setdim; /* Dimension of the set */ static fint Subdim; /* Dimension of the subset */ static fint Scrnum = 8; /* Destination of log output */ static fint Maxaxes = AXESMAX; /* Convert parameters to variables */ static fint Maxsubs = SUBSMAX; /* Max. input subsets */ static int subnrI; /* Index of current input subset */ static int i,j, m; /* Counters */ static fint Setlevel = 0; /* Indicate set level */ /* Output set related */ static fchar Setout; /* Name of the set */ static fint Subout[SUBSMAX]; /* Array for the subset coordinate words */ static fint Nsubsout; /* Number of input subsets */ static fint Axnumout[AXESMAX]; /* GDSOUT axis numbers array */ static fint Axcountout[AXESMAX]; /* GDSOUT axis lengths array */ /* Input of area etc.:*/ static fint CwloI; /* Coordinate words */ static fint CwhiI; static fint CwloO; static fint CwhiO; static fint GridloI[AXESMAX]; /* Coordinates for input-frame */ static fint GridhiI[AXESMAX]; static fint BgridloI[AXESMAX]; /* Coordinates for input-box */ static fint BgridhiI[AXESMAX]; static fint Boxopt; /* Default for gdsbox */ /* Data transfer: */ static fint Totpixels; /* Total # pixels in one subset */ static fint TidI; /* Tranfer id for input */ static fint TidO; /* Transfer id for output */ static float *ImageI = NULL; /* Multiple buffer for all subsets */ static fint Buflen; /* Max length of read buffer */ static int Bufheight; /* To determine max. buffer size */ static fint LenX, LenY; /* Size of user given box */ static float Zero; /* Variable containing 0.0 */ static fint Ylo, Yhi; /* Lines in read buffer */ static fint Readbuf, Afterread; /* Read buffer administration */ static fint Writepixels; /* Num. of pixels to write */ /* Calculations etc.; */ static fint Arraylen; /* Length used in presetting arrays */ static fint Order; /* User given order of polynomial */ static fint Coefficients; /* # coefficients in polynomial */ static float matrix[MAXELEM][MAXELEM]; /* least squares matrix */ static float B[MAXELEM]; /* Data vector */ static float SOL[MAXELEM]; /* The coefficients */ static int Xpos, Ypos; /* Two dimensional positions */ static int Xoffset, Yoffset; /* Used to determine 2-dim positions */ /* Related to update of header etc: */ static float Datamin[SUBSMAX]; /* Min, max of new set */ static float Datamax[SUBSMAX]; static fint Nblanks[SUBSMAX]; /* Number of blanks in new set */ static fint Remove; /* 0 means that minimum and maximum */ /* have not changed, anything else */ /* means that minimum and maximum */ /* have changed and that the MINMAX */ /* descriptors at intersecting levels */ /* will be removed by WMINMAX. */ /* Miscellaneous: */ static fint Numitems; /* Max. number of input items in userxxx routines */ static fint R1, R2; /* Results of userxxx routines */ static float Blank; /* Value of system blank */ static char messbuf[BIGSTORE]; /* Buffer for text message */ static fint Mcount; /* Counter for 'minmax3' routine */ static int agreed; /* Loop control */ static float Range[2]; /* User given data in/exclude range */ void anyoutC( char *anystr ) /*-------------------------------------------------------------------*/ /* The C version of 'anyout' needs a C type string as argument only. */ /* The value of 'Scrnum' is global. */ /*-------------------------------------------------------------------*/ { anyout_c( &Scrnum, tofchar( anystr ) ); } void polfill( int Xpos, int Ypos, float *B, fint Order ) /*---------------------------------------------------------------*/ /* Xpos, Ypos are the pixel positions, B is an array with */ /* solution parameters and Order is the order of the polynomial. */ /*---------------------------------------------------------------*/ { static int index; /* counter in polynomial vector */ static int i, j, k; index = 0; for (i = 0; i <= (int) Order; i++) { /* Loop over order */ for (j = 0; j <= i; j++) { /* loop over terms of that order */ B[index] = 1.0; for (k = 1; k <= (i-j); k++) { /* loop over X-powers */ B[index] *= (float) Xpos; } for (k = 1; k <= j; k++) { /* loop over Y-powers */ B[index] *= (float) Ypos; } index++; } } } void fillup( int Xpos, int Ypos, float Value, fint Coefficients ) /*-------------------------------------------------------------*/ /* Fill matrix used for lsq fit */ /*-------------------------------------------------------------*/ { static int i, j; static float C[MAXELEM]; (void) polfill( Xpos, Ypos, C, Order ); /* make from the row an array */ for(i = 0; i < Coefficients; i++) { /* loop over rows */ for (j = 0; j < Coefficients; j++) { /* loop over collums */ matrix[ j ][ i ] += C[ j ] * C[ i ]; } } /* Generate data vector */ for (i = 0; i < Coefficients; i++) { /* loop over rows */ B[ i ] += Value * C[ i ]; } } static fint invmat( fint nfree ) /*-------------------------------------------------------------------------*/ /* invmat calculates the inverse of matrix. The algorithm used is the */ /* Gauss-Jordan algorithm described in Stoer, Numerische matematik, 1 Teil.*/ /* The routine returns -6 if the determinant was zero. If successful, 0 is */ /* returned. !!!matrix is global!!! */ /*-------------------------------------------------------------------------*/ { #define MAXPAR 32 /* number of free parameters */ double even; double hv[MAXPAR]; double mjk; double rowmax; fint evin; fint i; fint j; fint k; fint per[MAXPAR]; fint row; for (i = 0; i < nfree; i++) per[i] = i; /* set permutation array */ for (j = 0; j < nfree; j++) { /* in j-th column, ... */ rowmax = fabs( matrix[j][j] ); /* determine row with ... */ row = j; /* largest element. */ for (i = j + 1; i < nfree; i++) { if (fabs( matrix[i][j] ) > rowmax) { rowmax = fabs( matrix[i][j] ); row = i; } } if (matrix[row][j] == 0.0) return( -6 ); /* determinant is zero! */ if (row > j) { /* if largest element not ... */ for (k = 0; k < nfree; k++) { /* on diagonal, then ... */ even = matrix[j][k]; /* permute rows. */ matrix[j][k] = matrix[row][k]; matrix[row][k] = even; } evin = per[j]; /* keep track of permutation */ per[j] = per[row]; per[row] = evin; } even = 1.0 / matrix[j][j]; /* modify column */ for (i = 0; i < nfree; i++) matrix[i][j] *= even; matrix[j][j] = even; for (k = 0; k < j; k++) { mjk = matrix[j][k]; for (i = 0; i < j; i++) matrix[i][k] -= matrix[i][j] * mjk; for (i = j + 1; i < nfree; i++) matrix[i][k] -= matrix[i][j] * mjk; matrix[j][k] = -even * mjk; } for (k = j + 1; k < nfree; k++) { mjk = matrix[j][k]; for (i = 0; i < j; i++) matrix[i][k] -= matrix[i][j] * mjk; for (i = j + 1; i < nfree; i++) matrix[i][k] -= matrix[i][j] * mjk; matrix[j][k] = -even * mjk; } } for (i = 0; i < nfree; i++) { /* finally, repermute the ... */ for (k = 0; k < nfree; k++) { /* columns. */ hv[per[k]] = matrix[i][k]; } for (k = 0; k < nfree; k++) { matrix[i][k] = hv[k]; } } return( 0 ); /* all is well */ } void matver( float *X, float *B, fint N ) /*----------------------------------------------------------------*/ /* This subroutine multiplies an array with a vector of length N */ /* The data array is given by matrix. The data vector by X and */ /* the result vector by B. matrix is global! */ /*----------------------------------------------------------------*/ { static int i, j; for (i = 0; i < (int) N; i++) { B[i] = 0.0; /* Reset result */ for (j = 0; j < (int) N; j++) { /* Loop over columns */ B[i] = B[i] + matrix[i][j] * X[j]; } } } void writeparms( float *SOL, fint Order ) /*---------------------------------------------*/ /* Output parameters to screen and Log file. */ /*---------------------------------------------*/ { static int index; static char longbuf[100]; static int len1, len2; index = 1; anyoutC("Result of polynomial plane fit"); anyoutC("------------------------------"); sprintf( messbuf, "Pol = %g", SOL[0] ); strcpy( longbuf, messbuf ); for (i = 1; i <= (int) Order; i++) { /* loop over orders from 1 */ for (j = 0; j <= i; j++) { /* loop over terms of that order */ if (SOL[index] < 0.0) { sprintf( messbuf, " - %g X**%d Y**%d", fabs(SOL[index]), i-j, j ); } else { sprintf( messbuf, " + %g X**%d Y**%d", SOL[index], i-j, j ); } len2 = strlen( messbuf ); len1 = strlen( longbuf ); if ( (len2+len1) <= 80 ) { strcat( longbuf, messbuf ); } else { anyoutC( longbuf ); strcpy( longbuf, messbuf ); } index++; } } anyoutC( longbuf ); anyoutC("------------------------------"); } float polval( int Xpos, int Ypos, float *SOL, fint Coefficients, fint Order ) /*-------------------------------------------------------------------*/ /* This function calculates the value of a two dimensional */ /* polynomial for position (X,Y). The order of the polynomial */ /* is given by Order. */ /*-------------------------------------------------------------------*/ { static float B[MAXELEM]; static float pol; static int i; (void) polfill( Xpos, Ypos, B, Order ); pol = 0.0; for (i = 0; i < Coefficients; i++) { pol += SOL[i] * B[i]; } return pol; } int inrange( float *value, float *Range ) /*----------------------------------------------------------*/ /* Check whether value of pixel is within user given range. */ /*----------------------------------------------------------*/ { if (Range[0] < Range[1]) { if (*value >= Range[0] && *value <= Range[1]) return true; else return false; } else { if (*value < Range[1] || *value > Range[0]) return true; else return false; } } MAIN_PROGRAM_ENTRY /*--------------------------------------------------------------------------*/ /* Because Fortran passes all arguments by reference, all C functions with */ /* a Fortran equivalent must do this also (GIPSY programmers guide, */ /* Chapter 9). Variables that can be interchanged between Fortran and C */ /* all start with a capital. */ /*--------------------------------------------------------------------------*/ { init_c(); /* contact Hermes */ /* Task identification */ { static fchar Task; /* Name of current task */ fmake( Task, 20 ); /* Macro 'fmake' must be available */ myname_c( Task ); /* Get task name */ Task.a[nelc_c(Task)] = '\0'; /* Terminate task name with null char. */ IDENTIFICATION( Task.a, VERSION ); /* Show task and version */ } setfblank_c( &Blank ); fmake(Setin, 80); Dfault = NONE; Scrnum = 8; /* Output device, 8 is experienced mode */ Subdim = 2; /* Dimension of subsets must be 2 */ Class = 1; Nsubsin = gdsinp_c( Setin, Subin, &Maxsubs, &Dfault, KEY_INSET, MES_INSET, &Scrnum, Axnum, Axcount, &Maxaxes, &Class, &Subdim ); Setdim = gdsc_ndims_c( Setin, &Setlevel ); /* *---------------------------------------------------------------------------- * Determine the edges of this its frame (GridloI/hiI) *---------------------------------------------------------------------------- */ R1 = 0; (void) gdsc_range_c( Setin, &Setlevel, &CwloI, &CwhiI, &R1 ); for (m = 0; m < (int) Setdim; m++) { R1 = R2 = 0; GridloI[m] = gdsc_grid_c( Setin, &Axnum[m], &CwloI, &R1 ); GridhiI[m] = gdsc_grid_c( Setin, &Axnum[m], &CwhiI, &R2 ); } /*--------------------------------------------------------*/ /* Prepare a box for INSET. Default is the entire subset. */ /*--------------------------------------------------------*/ Dfault = REQUEST; Boxopt = 0; Scrnum = 8; (void) gdsbox_c( BgridloI, BgridhiI, Setin, Subin, &Dfault, KEY_BOX, MES_BOX, &Scrnum, &Boxopt ); /* Count number of pixels in this substructure */ Totpixels = 1; /* For one subset */ for(m = 0; m < (int) Subdim; m++) Totpixels *= Axcount[m]; /*---------------------------------------------------------*/ /* Assign input characteristics to output and ask name for */ /* output. */ /*---------------------------------------------------------*/ (void) gdsasn_c( KEY_INSET, KEY_OUTSET, &Class ); fmake(Setout, BIGSTORE); Dfault = NONE; do { Nsubsout = gdsout_c( Setout, Subout, &Nsubsin, &Dfault, KEY_OUTSET, MES_OUTSET, &Scrnum, Axnumout, Axcountout, &Maxaxes ); agreed = (Nsubsout == Nsubsin); if (!agreed) { (void) reject_c( KEY_OUTSET, tofchar("Subset(s) error") ); } } while (!agreed); /*------------------------------------*/ /* Ask order of polynomial to be used */ /*------------------------------------*/ Order = 1; Dfault = REQUEST; Numitems = 1; do { R1 = userint_c( &Order, &Numitems, &Dfault, KEY_ORDER, MES_ORDER ); agreed = ((Order > 0) && (Order < 21)); if (!agreed) (void) reject_c( KEY_ORDER, tofchar("Wrong order!") ); } while (!agreed); /*------------------------------------------------------------------*/ /* Determine range of levels to include in the fit. The default */ /* values are the max. and min. floats of the system, so that all */ /* levels are included if the user pressed carriage return. */ /*------------------------------------------------------------------*/ Dfault = HIDDEN; Range[0] = -1.0*FLT_MAX; /* Defined in float.h */ Range[1] = FLT_MAX; getrange_c( Range, &Dfault, KEY_RANGE, MES_RANGE ); /*----------------------------------------------*/ /* Determine number of coefficients of function */ /*----------------------------------------------*/ Coefficients = 0; /* reset element counter */ for ( i = 0; i <= (int) Order; i++) { /* loop over all orders */ Coefficients += (i + 1); /* increase the # of elements by the order */ } /*------------------------------------------------------------*/ /* Try to find the maximum buffer size for the current system */ /* The total number of elements is given in 'Buflen' and the */ /* maximum number of lines is given in 'Bufheight' */ /*------------------------------------------------------------*/ LenX = (BgridhiI[0] - BgridloI[0] + 1); LenY = (BgridhiI[1] - BgridloI[1] + 1); Bufheight = LenY; Buflen = LenX * Bufheight; do { ImageI = (float *) calloc( (int) Buflen, sizeof(float) ); if (ImageI == NULL) { agreed = false; status_c(tofchar("Decreasing buffer size")); Bufheight /= 2; Buflen = LenX * Bufheight; } else agreed = true; } while( !agreed || (Bufheight < 1) ); if (Bufheight < 1) { static fint Errlev = 4; /* Generate fatal error */ anyoutC("Cannot allocate space for at least 1 line of input!"); error_c( &Errlev, tofchar("...Memory Allocation problems!...") ); } for (subnrI = 0; subnrI < Nsubsin; subnrI++) { /* Preset used arrays to zero */ Arraylen = MAXELEM; Zero = 0.0; presetr_c( &Zero, B, &Arraylen ); presetr_c( &Zero, SOL, &Arraylen ); Arraylen *= Arraylen; presetr_c( &Zero, &matrix[0][0], &Arraylen ); CwloI = gdsc_fill_c( Setin, &Subin[subnrI], BgridloI ); CwhiI = gdsc_fill_c( Setin, &Subin[subnrI], BgridhiI ); TidI = 0; Xoffset = Yoffset = 0; /*---------------------------------------------------------------*/ /* Create buffers with maximum size. If this size cannot contain */ /* The complete box, split it into smaller parts. Take as much */ /* parts as possible with size Bufheight*length and take care of */ /* leftovers. */ /*---------------------------------------------------------------*/ Ylo = BgridloI[1]; Yhi = BgridloI[1] + Bufheight - 1; while ( (Yhi-Ylo) > 0 ) { /* Read data in buffer of max size 'Bufheight' * length of box */ Readbuf = (Yhi - Ylo + 1) * LenX; (void) gdsi_read_c( Setin, &CwloI, &CwhiI, ImageI, &Readbuf, &Afterread, &TidI ); for (i = 0; i < (int) Afterread; i++) { if (Xoffset == LenX) { Yoffset++; /* Increase Y position */ Xoffset = 0; /* Reset Xoffset */ } Xpos = (int) BgridloI[0] + Xoffset; Ypos = (int) BgridloI[1] + Yoffset; Xoffset++; if (ImageI[i] != Blank) { if (inrange( &ImageI[i], Range )) { (void) fillup( Xpos, Ypos, ImageI[i], Coefficients ); } } } Ylo += Bufheight; Yhi += Bufheight; Yhi = MYMIN( BgridhiI[1], Yhi ); /* Yhi never exceeds Ymax in box */ } /*-----------------------------------------------------*/ /* We have: A * SOL = B ==> SOL = INV( A ) * B */ /* First invert matrix A and than multiply by vector B */ /*-----------------------------------------------------*/ R1 = invmat( Coefficients ); /* Matrix is global */ if (R1 != 0) { static fint Errlev = 4; /* Generate fatal error */ anyoutC("Cannot invert matrix, probably determinant = 0!"); error_c( &Errlev, tofchar("...Cannot invert matrix!...") ); } (void) matver( B, SOL, Coefficients ); /* Matrix is global */ (void) writeparms( SOL, Order ); /*------------------------------------------------------------*/ /* It is possible to put an arbitrary subset in the data set */ /* into just one output subset. The coordinate words will not */ /* correspond in this case, therefore we have to calculate */ /* the coordinate words of the output set separately. */ /*------------------------------------------------------------*/ CwloO = gdsc_fill_c( Setout, &Subout[subnrI], BgridloI ); CwhiO = gdsc_fill_c( Setout, &Subout[subnrI], BgridhiI ); TidO = 0; Mcount = 0; Xoffset = Yoffset = 0; m = 0; for (i = 0; i < (int) Totpixels; i++) { if (Xoffset == LenX) { Yoffset++; /* Increase Y position */ Xoffset = 0; /* Reset Xoffset */ } Xpos = (int) BgridloI[0] + Xoffset; Ypos = (int) BgridloI[1] + Yoffset; Xoffset++; /* Use the same buffer as for the input */ ImageI[m] = polval( Xpos, Ypos, SOL, Coefficients, Order ); if ( ((m+1) == Buflen) || ((i+1) == Totpixels) ) { /* Write to output if m is equal to Buflen or if this */ /* was the last pixel */ Writepixels = (fint) m + 1; (void) minmax3_c( ImageI, &Writepixels, &Datamin[subnrI], &Datamax[subnrI], &Nblanks[subnrI], &Mcount ); (void) gdsi_write_c( Setout, &CwloO, &CwhiO, ImageI, &Writepixels, &Writepixels, &TidO ); m = 0; /* Reset buflen. counter */ } else { m++; } } /* Give proceedings message */ sprintf( messbuf, "%d %% completed", 100*(subnrI+1)/Nsubsin ); status_c( tofchar(messbuf) ); } /* End of all subsets */ /* Remove in 'WMINMAX' old minmax descriptors because program */ /* changes values! */ Remove = true; (void) wminmax_c( Setout, Subout, Datamin, Datamax, Nblanks, &Nsubsout, &Remove ); finis_c(); /* Quit Hermes */ return 0; /* Dummy return for main */ }