diff options
Diffstat (limited to 'pkg/tbtables/cfitsio/histo.c')
| -rw-r--r-- | pkg/tbtables/cfitsio/histo.c | 1300 |
1 files changed, 1300 insertions, 0 deletions
diff --git a/pkg/tbtables/cfitsio/histo.c b/pkg/tbtables/cfitsio/histo.c new file mode 100644 index 00000000..525f50fd --- /dev/null +++ b/pkg/tbtables/cfitsio/histo.c @@ -0,0 +1,1300 @@ +/* Globally defined histogram parameters */ +#include <string.h> +#include <ctype.h> +#include <math.h> +#include <stdlib.h> +#include "fitsio2.h" + +typedef struct { /* Structure holding all the histogramming information */ + union { /* the iterator work functions (ffwritehist, ffcalchist) */ + char *b; /* need to do their job... passed via *userPointer. */ + short *i; + int *j; + float *r; + double *d; + } hist; + + fitsfile *tblptr; + + int haxis, hcolnum[4], himagetype; + long haxis1, haxis2, haxis3, haxis4; + float amin1, amin2, amin3, amin4; + float maxbin1, maxbin2, maxbin3, maxbin4; + float binsize1, binsize2, binsize3, binsize4; + int wtrecip, wtcolnum; + float weight; + char *rowselector; + +} histType; + +/*--------------------------------------------------------------------------*/ +int ffbins(char *binspec, /* I - binning specification */ + int *imagetype, /* O - image type, TINT or TSHORT */ + int *histaxis, /* O - no. of axes in the histogram */ + char colname[4][FLEN_VALUE], /* column name for axis */ + double *minin, /* minimum value for each axis */ + double *maxin, /* maximum value for each axis */ + double *binsizein, /* size of bins on each axis */ + char minname[4][FLEN_VALUE], /* keyword name for min */ + char maxname[4][FLEN_VALUE], /* keyword name for max */ + char binname[4][FLEN_VALUE], /* keyword name for binsize */ + double *wt, /* weighting factor */ + char *wtname, /* keyword or column name for weight */ + int *recip, /* the reciprocal of the weight? */ + int *status) +{ +/* + Parse the input binning specification string, returning the binning + parameters. Supports up to 4 dimensions. The binspec string has + one of these forms: + + bin binsize - 2D histogram with binsize on each axis + bin xcol - 1D histogram on column xcol + bin (xcol, ycol) = binsize - 2D histogram with binsize on each axis + bin x=min:max:size, y=min:max:size, z..., t... + bin x=:max, y=::size + bin x=size, y=min::size + + most other reasonable combinations are supported. +*/ + int ii, slen, defaulttype; + char *ptr, tmpname[30], *file_expr = NULL; + double dummy; + + if (*status > 0) + return(*status); + + /* set the default values */ + *histaxis = 2; + *imagetype = TINT; + defaulttype = 1; + *wt = 1.; + *recip = 0; + *wtname = '\0'; + + /* set default values */ + for (ii = 0; ii < 4; ii++) + { + *colname[ii] = '\0'; + *minname[ii] = '\0'; + *maxname[ii] = '\0'; + *binname[ii] = '\0'; + minin[ii] = DOUBLENULLVALUE; /* undefined values */ + maxin[ii] = DOUBLENULLVALUE; + binsizein[ii] = DOUBLENULLVALUE; + } + + ptr = binspec + 3; /* skip over 'bin' */ + + if (*ptr == 'i' ) /* bini */ + { + *imagetype = TSHORT; + defaulttype = 0; + ptr++; + } + else if (*ptr == 'j' ) /* binj; same as default */ + { + defaulttype = 0; + ptr ++; + } + else if (*ptr == 'r' ) /* binr */ + { + *imagetype = TFLOAT; + defaulttype = 0; + ptr ++; + } + else if (*ptr == 'd' ) /* bind */ + { + *imagetype = TDOUBLE; + defaulttype = 0; + ptr ++; + } + else if (*ptr == 'b' ) /* binb */ + { + *imagetype = TBYTE; + defaulttype = 0; + ptr ++; + } + + if (*ptr == '\0') /* use all defaults for other parameters */ + return(*status); + else if (*ptr != ' ') /* must be at least one blank */ + { + ffpmsg("binning specification syntax error:"); + ffpmsg(binspec); + return(*status = URL_PARSE_ERROR); + } + + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + if (*ptr == '\0') /* no other parameters; use defaults */ + return(*status); + + /* Check if need to import expression from a file */ + + if( *ptr=='@' ) { + if( ffimport_file( ptr+1, &file_expr, status ) ) return(*status); + ptr = file_expr; + while (*ptr == ' ') + ptr++; /* skip leading white space... again */ + } + + if (*ptr == '(' ) + { + /* this must be the opening parenthesis around a list of column */ + /* names, optionally followed by a '=' and the binning spec. */ + + for (ii = 0; ii < 4; ii++) + { + ptr++; /* skip over the '(', ',', or ' ') */ + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + slen = strcspn(ptr, " ,)"); + strncat(colname[ii], ptr, slen); /* copy 1st column name */ + + ptr += slen; + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + if (*ptr == ')' ) /* end of the list of names */ + { + *histaxis = ii + 1; + break; + } + } + + if (ii == 4) /* too many names in the list , or missing ')' */ + { + ffpmsg( + "binning specification has too many column names or is missing closing ')':"); + ffpmsg(binspec); + if( file_expr ) free( file_expr ); + return(*status = URL_PARSE_ERROR); + } + + ptr++; /* skip over the closing parenthesis */ + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + if (*ptr == '\0') { + if( file_expr ) free( file_expr ); + return(*status); /* parsed the entire string */ + } + + else if (*ptr != '=') /* must be an equals sign now*/ + { + ffpmsg("illegal binning specification in URL:"); + ffpmsg(" an equals sign '=' must follow the column names"); + ffpmsg(binspec); + if( file_expr ) free( file_expr ); + return(*status = URL_PARSE_ERROR); + } + + ptr++; /* skip over the equals sign */ + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + /* get the single range specification for all the columns */ + ffbinr(&ptr, tmpname, minin, + maxin, binsizein, minname[0], + maxname[0], binname[0], status); + if (*status > 0) + { + ffpmsg("illegal binning specification in URL:"); + ffpmsg(binspec); + if( file_expr ) free( file_expr ); + return(*status); + } + + for (ii = 1; ii < *histaxis; ii++) + { + minin[ii] = minin[0]; + maxin[ii] = maxin[0]; + binsizein[ii] = binsizein[0]; + strcpy(minname[ii], minname[0]); + strcpy(maxname[ii], maxname[0]); + strcpy(binname[ii], binname[0]); + } + + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + if (*ptr == ';') + goto getweight; /* a weighting factor is specified */ + + if (*ptr != '\0') /* must have reached end of string */ + { + ffpmsg("illegal syntax after binning range specification in URL:"); + ffpmsg(binspec); + if( file_expr ) free( file_expr ); + return(*status = URL_PARSE_ERROR); + } + + return(*status); + } /* end of case with list of column names in ( ) */ + + /* if we've reached this point, then the binning specification */ + /* must be of the form: XCOL = min:max:binsize, YCOL = ... */ + /* where the column name followed by '=' are optional. */ + /* If the column name is not specified, then use the default name */ + + for (ii = 0; ii < 4; ii++) /* allow up to 4 histogram dimensions */ + { + ffbinr(&ptr, colname[ii], &minin[ii], + &maxin[ii], &binsizein[ii], minname[ii], + maxname[ii], binname[ii], status); + + if (*status > 0) + { + ffpmsg("illegal syntax in binning range specification in URL:"); + ffpmsg(binspec); + if( file_expr ) free( file_expr ); + return(*status); + } + + if (*ptr == '\0' || *ptr == ';') + break; /* reached the end of the string */ + + if (*ptr == ' ') + { + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + if (*ptr == '\0' || *ptr == ';') + break; /* reached the end of the string */ + + if (*ptr == ',') + ptr++; /* comma separates the next column specification */ + } + else if (*ptr == ',') + { + ptr++; /* comma separates the next column specification */ + } + else + { + ffpmsg("illegal characters following binning specification in URL:"); + ffpmsg(binspec); + if( file_expr ) free( file_expr ); + return(*status = URL_PARSE_ERROR); + } + } + + if (ii == 4) + { + /* there are yet more characters in the string */ + ffpmsg("illegal binning specification in URL:"); + ffpmsg("apparently greater than 4 histogram dimensions"); + ffpmsg(binspec); + return(*status = URL_PARSE_ERROR); + } + else + *histaxis = ii + 1; + + /* special case: if a single number was entered it should be */ + /* interpreted as the binning factor for the default X and Y axes */ + + if (*histaxis == 1 && *colname[0] == '\0' && + minin[0] == DOUBLENULLVALUE && maxin[0] == DOUBLENULLVALUE) + { + *histaxis = 2; + binsizein[1] = binsizein[0]; + } + +getweight: + if (*ptr == ';') /* looks like a weighting factor is given */ + { + ptr++; + + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + recip = 0; + if (*ptr == '/') + { + *recip = 1; /* the reciprocal of the weight is entered */ + ptr++; + + while (*ptr == ' ') /* skip over blanks */ + ptr++; + } + + /* parse the weight as though it were a binrange. */ + /* either a column name or a numerical value will be returned */ + + ffbinr(&ptr, wtname, &dummy, &dummy, wt, tmpname, + tmpname, tmpname, status); + + if (*status > 0) + { + ffpmsg("illegal binning weight specification in URL:"); + ffpmsg(binspec); + if( file_expr ) free( file_expr ); + return(*status); + } + + /* creat a float datatype histogram by default, if weight */ + /* factor is not = 1.0 */ + + if (defaulttype && *wt != 1.0) + *imagetype = TFLOAT; + } + + while (*ptr == ' ') /* skip over blanks */ + ptr++; + + if (*ptr != '\0') /* should have reached the end of string */ + { + ffpmsg("illegal syntax after binning weight specification in URL:"); + ffpmsg(binspec); + *status = URL_PARSE_ERROR; + } + + if( file_expr ) free( file_expr ); + return(*status); +} +/*--------------------------------------------------------------------------*/ +int ffbinr(char **ptr, + char *colname, + double *minin, + double *maxin, + double *binsizein, + char *minname, + char *maxname, + char *binname, + int *status) +/* + Parse the input binning range specification string, returning + the column name, histogram min and max values, and bin size. +*/ +{ + int slen, isanumber; + char token[FLEN_VALUE]; + + if (*status > 0) + return(*status); + + slen = fits_get_token(ptr, " ,=:;", token, &isanumber); /* get 1st token */ + + if (slen == 0 && (**ptr == '\0' || **ptr == ',' || **ptr == ';') ) + return(*status); /* a null range string */ + + if (!isanumber && **ptr != ':') + { + /* this looks like the column name */ + + if (token[0] == '#' && isdigit((int) token[1]) ) + { + /* omit the leading '#' in the column number */ + strcpy(colname, token+1); + } + else + strcpy(colname, token); + + while (**ptr == ' ') /* skip over blanks */ + (*ptr)++; + + if (**ptr != '=') + return(*status); /* reached the end */ + + (*ptr)++; /* skip over the = sign */ + + while (**ptr == ' ') /* skip over blanks */ + (*ptr)++; + + slen = fits_get_token(ptr, " ,:;", token, &isanumber); /* get token */ + } + + if (**ptr != ':') + { + /* this is the first token, and since it is not followed by */ + /* a ':' this must be the binsize token */ + if (!isanumber) + strcpy(binname, token); + else + *binsizein = strtod(token, NULL); + + return(*status); /* reached the end */ + } + else + { + /* the token contains the min value */ + if (slen) + { + if (!isanumber) + strcpy(minname, token); + else + *minin = strtod(token, NULL); + } + } + + (*ptr)++; /* skip the colon between the min and max values */ + slen = fits_get_token(ptr, " ,:;", token, &isanumber); /* get token */ + + /* the token contains the max value */ + if (slen) + { + if (!isanumber) + strcpy(maxname, token); + else + *maxin = strtod(token, NULL); + } + + if (**ptr != ':') + return(*status); /* reached the end; no binsize token */ + + (*ptr)++; /* skip the colon between the max and binsize values */ + slen = fits_get_token(ptr, " ,:;", token, &isanumber); /* get token */ + + /* the token contains the binsize value */ + if (slen) + { + if (!isanumber) + strcpy(binname, token); + else + *binsizein = strtod(token, NULL); + } + + return(*status); +} +/*--------------------------------------------------------------------------*/ +int ffhist(fitsfile **fptr, /* IO - pointer to table with X and Y cols; */ + /* on output, points to histogram image */ + char *outfile, /* I - name for the output histogram file */ + int imagetype, /* I - datatype for image: TINT, TSHORT, etc */ + int naxis, /* I - number of axes in the histogram image */ + char colname[4][FLEN_VALUE], /* I - column names */ + double *minin, /* I - minimum histogram value, for each axis */ + double *maxin, /* I - maximum histogram value, for each axis */ + double *binsizein, /* I - bin size along each axis */ + char minname[4][FLEN_VALUE], /* I - optional keywords for min */ + char maxname[4][FLEN_VALUE], /* I - optional keywords for max */ + char binname[4][FLEN_VALUE], /* I - optional keywords for binsize */ + double weightin, /* I - binning weighting factor */ + char wtcol[FLEN_VALUE], /* I - optional keyword or col for weight*/ + int recip, /* I - use reciprocal of the weight? */ + char *selectrow, /* I - optional array (length = no. of */ + /* rows in the table). If the element is true */ + /* then the corresponding row of the table will*/ + /* be included in the histogram, otherwise the */ + /* row will be skipped. Ingnored if *selectrow*/ + /* is equal to NULL. */ + int *status) +{ + int ii, datatype, repeat, imin, imax, ibin, bitpix, tstatus, use_datamax = 0; + long haxes[4]; + fitsfile *histptr; + char errmsg[FLEN_ERRMSG], keyname[FLEN_KEYWORD], card[FLEN_CARD]; + tcolumn *colptr; + iteratorCol imagepars[1]; + int n_cols = 1, nkeys; + long offset = 0; + long n_per_loop = -1; /* force whole array to be passed at one time */ + histType histData; /* Structure holding histogram info for iterator */ + + float amin[4], amax[4], binsize[4], maxbin[4]; + float datamin = FLOATNULLVALUE, datamax = FLOATNULLVALUE; + char svalue[FLEN_VALUE]; + double dvalue; + char cpref[4][FLEN_VALUE]; + char *cptr; + + if (*status > 0) + return(*status); + + if (naxis > 4) + { + ffpmsg("histogram has more than 4 dimensions"); + return(*status = BAD_DIMEN); + } + + /* reset position to the correct HDU if necessary */ + if ((*fptr)->HDUposition != ((*fptr)->Fptr)->curhdu) + ffmahd(*fptr, ((*fptr)->HDUposition) + 1, NULL, status); + + histData.tblptr = *fptr; + histData.himagetype = imagetype; + histData.haxis = naxis; + histData.rowselector = selectrow; + + if (imagetype == TBYTE) + bitpix = BYTE_IMG; + else if (imagetype == TSHORT) + bitpix = SHORT_IMG; + else if (imagetype == TINT) + bitpix = LONG_IMG; + else if (imagetype == TFLOAT) + bitpix = FLOAT_IMG; + else if (imagetype == TDOUBLE) + bitpix = DOUBLE_IMG; + else + return(*status = BAD_DATATYPE); + + /* The CPREF keyword, if it exists, gives the preferred columns. */ + /* Otherwise, assume "X", "Y", "Z", and "T" */ + + tstatus = 0; + ffgky(*fptr, TSTRING, "CPREF", cpref[0], NULL, &tstatus); + + if (!tstatus) + { + /* Preferred column names are given; separate them */ + cptr = cpref[0]; + + /* the first preferred axis... */ + while (*cptr != ',' && *cptr != '\0') + cptr++; + + if (*cptr != '\0') + { + *cptr = '\0'; + cptr++; + while (*cptr == ' ') + cptr++; + + strcpy(cpref[1], cptr); + cptr = cpref[1]; + + /* the second preferred axis... */ + while (*cptr != ',' && *cptr != '\0') + cptr++; + + if (*cptr != '\0') + { + *cptr = '\0'; + cptr++; + while (*cptr == ' ') + cptr++; + + strcpy(cpref[2], cptr); + cptr = cpref[2]; + + /* the third preferred axis... */ + while (*cptr != ',' && *cptr != '\0') + cptr++; + + if (*cptr != '\0') + { + *cptr = '\0'; + cptr++; + while (*cptr == ' ') + cptr++; + + strcpy(cpref[3], cptr); + + } + } + } + } + + for (ii = 0; ii < naxis; ii++) + { + + /* get the min, max, and binsize values from keywords, if specified */ + + if (*minname[ii]) + { + if (ffgky(*fptr, TDOUBLE, minname[ii], &minin[ii], NULL, status) ) + { + ffpmsg("error reading histogramming minimum keyword"); + ffpmsg(minname[ii]); + return(*status); + } + } + + if (*maxname[ii]) + { + if (ffgky(*fptr, TDOUBLE, maxname[ii], &maxin[ii], NULL, status) ) + { + ffpmsg("error reading histogramming maximum keyword"); + ffpmsg(maxname[ii]); + return(*status); + } + } + + if (*binname[ii]) + { + if (ffgky(*fptr, TDOUBLE, binname[ii], &binsizein[ii], NULL, status) ) + { + ffpmsg("error reading histogramming binsize keyword"); + ffpmsg(binname[ii]); + return(*status); + } + } + + if (binsizein[ii] == 0.) + { + ffpmsg("error: histogram binsize = 0"); + return(*status = ZERO_SCALE); + } + + if (*colname[ii] == '\0') + { + strcpy(colname[ii], cpref[ii]); /* try using the preferred column */ + if (*colname[ii] == '\0') + { + if (ii == 0) + strcpy(colname[ii], "X"); + else if (ii == 1) + strcpy(colname[ii], "Y"); + else if (ii == 2) + strcpy(colname[ii], "Z"); + else if (ii == 3) + strcpy(colname[ii], "T"); + } + } + + /* get the column number in the table */ + if (ffgcno(*fptr, CASEINSEN, colname[ii], histData.hcolnum+ii, status) + > 0) + { + strcpy(errmsg, "column for histogram axis doesn't exist: "); + strcat(errmsg, colname[ii]); + ffpmsg(errmsg); + return(*status); + } + + colptr = ((*fptr)->Fptr)->tableptr; + colptr += (histData.hcolnum[ii] - 1); + + repeat = colptr->trepeat; /* vector repeat factor of the column */ + if (repeat > 1) + { + strcpy(errmsg, "Can't bin a vector column: "); + strcat(errmsg, colname[ii]); + ffpmsg(errmsg); + return(*status = BAD_DATATYPE); + } + + /* get the datatype of the column */ + fits_get_coltype(*fptr, histData.hcolnum[ii], &datatype, + NULL, NULL, status); + + if (datatype < 0 || datatype == TSTRING) + { + strcpy(errmsg, "Inappropriate datatype; can't bin this column: "); + strcat(errmsg, colname[ii]); + ffpmsg(errmsg); + return(*status = BAD_DATATYPE); + } + + /* use TLMINn and TLMAXn keyword values if min and max were not given */ + /* else use actual data min and max if TLMINn and TLMAXn don't exist */ + + if (minin[ii] == DOUBLENULLVALUE) + { + ffkeyn("TLMIN", histData.hcolnum[ii], keyname, status); + if (ffgky(*fptr, TFLOAT, keyname, amin+ii, NULL, status) > 0) + { + /* use actual data minimum value for the histogram minimum */ + *status = 0; + if (fits_get_col_minmax(*fptr, histData.hcolnum[ii], amin+ii, &datamax, status) > 0) + { + strcpy(errmsg, "Error calculating datamin and datamax for column: "); + strcat(errmsg, colname[ii]); + ffpmsg(errmsg); + return(*status); + } + } + } + else + { + amin[ii] = minin[ii]; + } + + if (maxin[ii] == DOUBLENULLVALUE) + { + ffkeyn("TLMAX", histData.hcolnum[ii], keyname, status); + if (ffgky(*fptr, TFLOAT, keyname, &amax[ii], NULL, status) > 0) + { + *status = 0; + if(datamax != FLOATNULLVALUE) /* already computed max value */ + { + amax[ii] = datamax; + } + else + { + /* use actual data maximum value for the histogram maximum */ + if (fits_get_col_minmax(*fptr, histData.hcolnum[ii], &datamin, &amax[ii], status) > 0) + { + strcpy(errmsg, "Error calculating datamin and datamax for column: "); + strcat(errmsg, colname[ii]); + ffpmsg(errmsg); + return(*status); + } + } + } + use_datamax = 1; /* flag that the max was determined by the data values */ + /* and not specifically set by the calling program */ + } + else + { + amax[ii] = maxin[ii]; + } + + /* use TDBINn keyword or else 1 if bin size is not given */ + if (binsizein[ii] == DOUBLENULLVALUE) + { + tstatus = 0; + ffkeyn("TDBIN", histData.hcolnum[ii], keyname, &tstatus); + + if (ffgky(*fptr, TDOUBLE, keyname, binsizein + ii, NULL, &tstatus) > 0) + { + /* make at least 10 bins */ + binsizein[ii] = (amax[ii] - amin[ii]) / 10. ; + if (binsizein[ii] > 1.) + binsizein[ii] = 1.; /* use default bin size */ + } + } + + if ( (amin[ii] > amax[ii] && binsizein[ii] > 0. ) || + (amin[ii] < amax[ii] && binsizein[ii] < 0. ) ) + binsize[ii] = -binsizein[ii]; /* reverse the sign of binsize */ + else + binsize[ii] = binsizein[ii]; /* binsize has the correct sign */ + + ibin = binsize[ii]; + imin = amin[ii]; + imax = amax[ii]; + + /* Determine the range and number of bins in the histogram. This */ + /* depends on whether the input columns are integer or floats, so */ + /* treat each case separately. */ + + if (datatype <= TLONG && (float) imin == amin[ii] && + (float) imax == amax[ii] && + (float) ibin == binsize[ii] ) + { + /* This is an integer column and integer limits were entered. */ + /* Shift the lower and upper histogramming limits by 0.5, so that */ + /* the values fall in the center of the bin, not on the edge. */ + + haxes[ii] = (imax - imin) / ibin + 1; /* last bin may only */ + /* be partially full */ + maxbin[ii] = haxes[ii] + 1.; /* add 1. instead of .5 to avoid roundoff */ + + if (amin[ii] < amax[ii]) + { + amin[ii] = amin[ii] - 0.5; + amax[ii] = amax[ii] + 0.5; + } + else + { + amin[ii] = amin[ii] + 0.5; + amax[ii] = amax[ii] - 0.5; + } + } + else if (use_datamax) + { + /* Either the column datatype and/or the limits are floating point, */ + /* and the histogram limits are being defined by the min and max */ + /* values of the array. Add 1 to the number of histogram bins to */ + /* make sure that pixels that are equal to the maximum or are */ + /* in the last partial bin are included. */ + + maxbin[ii] = (amax[ii] - amin[ii]) / binsize[ii]; + haxes[ii] = maxbin[ii] + 1; + } + else + { + /* float datatype column and/or limits, and the maximum value to */ + /* include in the histogram is specified by the calling program. */ + /* The lower limit is inclusive, but upper limit is exclusive */ + maxbin[ii] = (amax[ii] - amin[ii]) / binsize[ii]; + haxes[ii] = maxbin[ii]; + + if (amin[ii] < amax[ii]) + { + if (amin[ii] + (haxes[ii] * binsize[ii]) < amax[ii]) + haxes[ii]++; /* need to include another partial bin */ + } + else + { + if (amin[ii] + (haxes[ii] * binsize[ii]) > amax[ii]) + haxes[ii]++; /* need to include another partial bin */ + } + } + } + + /* get the histogramming weighting factor */ + if (*wtcol) + { + /* first, look for a keyword with the weight value */ + if (ffgky(*fptr, TFLOAT, wtcol, &histData.weight, NULL, status) ) + { + /* not a keyword, so look for column with this name */ + *status = 0; + + /* get the column number in the table */ + if (ffgcno(*fptr, CASEINSEN, wtcol, &histData.wtcolnum, status) > 0) + { + ffpmsg( + "keyword or column for histogram weights doesn't exist: "); + ffpmsg(wtcol); + return(*status); + } + + histData.weight = FLOATNULLVALUE; + } + } + else + histData.weight = (float) weightin; + + if (histData.weight <= 0. && histData.weight != FLOATNULLVALUE) + { + ffpmsg("Illegal histogramming weighting factor <= 0."); + return(*status = URL_PARSE_ERROR); + } + + if (recip && histData.weight != FLOATNULLVALUE) + /* take reciprocal of weight */ + histData.weight = 1.0 / histData.weight; + + histData.wtrecip = recip; + + /* size of histogram is now known, so create temp output file */ + if (ffinit(&histptr, outfile, status) > 0) + { + ffpmsg("failed to create temp output file for histogram"); + return(*status); + } + + if (ffcrim(histptr, bitpix, histData.haxis, haxes, status) > 0) + { + ffpmsg("failed to create primary array histogram in temp file"); + ffclos(histptr, status); + return(*status); + } + + /* copy all non-structural keywords from the table to the image */ + fits_get_hdrspace(*fptr, &nkeys, NULL, status); + for (ii = 1; ii <= nkeys; ii++) + { + fits_read_record(*fptr, ii, card, status); + if (fits_get_keyclass(card) >= 120) + fits_write_record(histptr, card, status); + } + + /* Set global variables with histogram parameter values. */ + /* Use separate scalar variables rather than arrays because */ + /* it is more efficient when computing the histogram. */ + + histData.amin1 = amin[0]; + histData.maxbin1 = maxbin[0]; + histData.binsize1 = binsize[0]; + histData.haxis1 = haxes[0]; + + if (histData.haxis > 1) + { + histData.amin2 = amin[1]; + histData.maxbin2 = maxbin[1]; + histData.binsize2 = binsize[1]; + histData.haxis2 = haxes[1]; + + if (histData.haxis > 2) + { + histData.amin3 = amin[2]; + histData.maxbin3 = maxbin[2]; + histData.binsize3 = binsize[2]; + histData.haxis3 = haxes[2]; + + if (histData.haxis > 3) + { + histData.amin4 = amin[3]; + histData.maxbin4 = maxbin[3]; + histData.binsize4 = binsize[3]; + histData.haxis4 = haxes[3]; + } + } + } + + /* define parameters of image for the iterator function */ + fits_iter_set_file(imagepars, histptr); /* pointer to image */ + fits_iter_set_datatype(imagepars, imagetype); /* image datatype */ + fits_iter_set_iotype(imagepars, OutputCol); /* image is output */ + + /* call the iterator function to write out the histogram image */ + if (fits_iterate_data(n_cols, imagepars, offset, n_per_loop, + ffwritehisto, (void*)&histData, status) ) + return(*status); + + /* write the World Coordinate System (WCS) keywords */ + /* create default values if WCS keywords are not present in the table */ + for (ii = 0; ii < histData.haxis; ii++) + { + /* CTYPEn */ + tstatus = 0; + ffkeyn("TCTYP", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus); + if (tstatus) + { /* just use column name as the type */ + tstatus = 0; + ffkeyn("TTYPE", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus); + } + + if (!tstatus) + { + ffkeyn("CTYPE", ii + 1, keyname, &tstatus); + ffpky(histptr, TSTRING, keyname, svalue, "Coordinate Type", &tstatus); + } + else + tstatus = 0; + + /* CUNITn */ + ffkeyn("TCUNI", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus); + if (tstatus) + { /* use the column units */ + tstatus = 0; + ffkeyn("TUNIT", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus); + } + + if (!tstatus) + { + ffkeyn("CUNIT", ii + 1, keyname, &tstatus); + ffpky(histptr, TSTRING, keyname, svalue, "Coordinate Units", &tstatus); + } + else + tstatus = 0; + + /* CRPIXn - Reference Pixel */ + ffkeyn("TCRPX", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); + if (tstatus) + { + dvalue = 1.0; /* choose first pixel in new image as ref. pix. */ + tstatus = 0; + } + else + { + /* calculate locate of the ref. pix. in the new image */ + dvalue = (dvalue - amin[ii]) / binsize[ii] + .5; + } + + ffkeyn("CRPIX", ii + 1, keyname, &tstatus); + ffpky(histptr, TDOUBLE, keyname, &dvalue, "Reference Pixel", &tstatus); + + /* CRVALn - Value at the location of the reference pixel */ + ffkeyn("TCRVL", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); + if (tstatus) + { + /* calculate value at ref. pix. location (at center of 1st pixel) */ + dvalue = amin[ii] + binsize[ii]/2.; + tstatus = 0; + } + + ffkeyn("CRVAL", ii + 1, keyname, &tstatus); + ffpky(histptr, TDOUBLE, keyname, &dvalue, "Reference Value", &tstatus); + + /* CDELTn - unit size of pixels */ + ffkeyn("TCDLT", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); + if (tstatus) + { + dvalue = 1.0; /* use default pixel size */ + tstatus = 0; + } + + dvalue = dvalue * binsize[ii]; + ffkeyn("CDELT", ii + 1, keyname, &tstatus); + ffpky(histptr, TDOUBLE, keyname, &dvalue, "Pixel size", &tstatus); + + /* CROTAn - Rotation angle (degrees CCW) */ + /* There should only be a CROTA2 keyword, and only for 2+ D images */ + if (ii == 1) + { + ffkeyn("TCROT", histData.hcolnum[ii], keyname, &tstatus); + ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); + if (!tstatus && dvalue != 0.) /* only write keyword if angle != 0 */ + { + ffkeyn("CROTA", ii + 1, keyname, &tstatus); + ffpky(histptr, TDOUBLE, keyname, &dvalue, + "Rotation angle", &tstatus); + } + else + { + /* didn't find CROTA for the 2nd axis, so look for one */ + /* on the first axis */ + tstatus = 0; + ffkeyn("TCROT", histData.hcolnum[0], keyname, &tstatus); + ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); + if (!tstatus && dvalue != 0.) /* only write keyword if angle != 0 */ + { + dvalue *= -1.; /* negate the value, because mirror image */ + ffkeyn("CROTA", ii + 1, keyname, &tstatus); + ffpky(histptr, TDOUBLE, keyname, &dvalue, + "Rotation angle", &tstatus); + } + } + } + } + + /* finally, close the original file and return ptr to the new image */ + ffclos(*fptr, status); + *fptr = histptr; + + return(*status); +} +/*--------------------------------------------------------------------------*/ +int fits_get_col_minmax(fitsfile *fptr, int colnum, float *datamin, + float *datamax, int *status) +/* + Simple utility routine to compute the min and max value in a column +*/ +{ + int anynul; + long nrows, ntodo, firstrow, ii; + float array[1000], nulval; + + ffgky(fptr, TLONG, "NAXIS2", &nrows, NULL, status); /* no. of rows */ + + firstrow = 1; + nulval = FLOATNULLVALUE; + *datamin = 9.0E36F; + *datamax = -9.0E36F; + + while(nrows) + { + ntodo = minvalue(nrows, 100); + ffgcv(fptr, TFLOAT, colnum, firstrow, 1, ntodo, &nulval, array, + &anynul, status); + + for (ii = 0; ii < ntodo; ii++) + { + if (array[ii] != nulval) + { + *datamin = minvalue(*datamin, array[ii]); + *datamax = maxvalue(*datamax, array[ii]); + } + } + + nrows -= ntodo; + firstrow += ntodo; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +int ffwritehisto(long totaln, long pixoffset, long firstn, long nvalues, + int narrays, iteratorCol *imagepars, void *userPointer) +/* + Interator work function that writes out the histogram. + The histogram values are calculated by another work function, ffcalchisto. + This work function only gets called once, and totaln = nvalues. +*/ +{ + iteratorCol colpars[5]; + int ii, status = 0, ncols; + long rows_per_loop = 0, offset = 0; + histType *histData; + + histData = (histType *)userPointer; + + /* store pointer to the histogram array, and initialize to zero */ + + switch( histData->himagetype ) { + case TBYTE: + histData->hist.b = (char * ) fits_iter_get_array(imagepars); + break; + case TSHORT: + histData->hist.i = (short * ) fits_iter_get_array(imagepars); + break; + case TINT: + histData->hist.j = (int * ) fits_iter_get_array(imagepars); + break; + case TFLOAT: + histData->hist.r = (float * ) fits_iter_get_array(imagepars); + break; + case TDOUBLE: + histData->hist.d = (double *) fits_iter_get_array(imagepars); + break; + } + + /* set the column parameters for the iterator function */ + for (ii = 0; ii < histData->haxis; ii++) + { + fits_iter_set_by_num(&colpars[ii], histData->tblptr, + histData->hcolnum[ii], TFLOAT, InputCol); + } + ncols = histData->haxis; + + if (histData->weight == FLOATNULLVALUE) + { + fits_iter_set_by_num(&colpars[histData->haxis], histData->tblptr, + histData->wtcolnum, TFLOAT, InputCol); + ncols = histData->haxis + 1; + } + + /* call iterator function to calc the histogram pixel values */ + fits_iterate_data(ncols, colpars, offset, rows_per_loop, + ffcalchist, (void*)histData, &status); + + return(status); +} +/*--------------------------------------------------------------------------*/ +int ffcalchist(long totalrows, long offset, long firstrow, long nrows, + int ncols, iteratorCol *colpars, void *userPointer) +/* + Interator work function that calculates values for the 2D histogram. +*/ +{ + long ii, ipix, iaxisbin; + float pix, axisbin; + static float *col1, *col2, *col3, *col4; /* static to preserve values */ + static float *wtcol; + static long incr2, incr3, incr4; + static histType histData; + static char *rowselect; + + /* Initialization procedures: execute on the first call */ + if (firstrow == 1) + { + + /* Copy input histogram data to static local variable so we */ + /* don't have to constantly dereference it. */ + + histData = *(histType*)userPointer; + rowselect = histData.rowselector; + + /* assign the input array pointers to local pointers */ + col1 = (float *) fits_iter_get_array(&colpars[0]); + if (histData.haxis > 1) + { + col2 = (float *) fits_iter_get_array(&colpars[1]); + incr2 = histData.haxis1; + + if (histData.haxis > 2) + { + col3 = (float *) fits_iter_get_array(&colpars[2]); + incr3 = incr2 * histData.haxis2; + + if (histData.haxis > 3) + { + col4 = (float *) fits_iter_get_array(&colpars[3]); + incr4 = incr3 * histData.haxis3; + } + } + } + + if (ncols > histData.haxis) /* then weights are give in a column */ + { + wtcol = (float *) fits_iter_get_array(&colpars[histData.haxis]); + } + } /* end of Initialization procedures */ + + /* Main loop: increment the histogram at position of each event */ + for (ii = 1; ii <= nrows; ii++) + { + if (rowselect) /* if a row selector array is supplied... */ + { + if (*rowselect) + { + rowselect++; /* this row is included in the histogram */ + } + else + { + rowselect++; /* this row is excluded from the histogram */ + continue; + } + } + + if (col1[ii] == FLOATNULLVALUE) /* test for null value */ + continue; + + pix = (col1[ii] - histData.amin1) / histData.binsize1; + ipix = (long) (pix + 1.); /* add 1 because the 1st pixel is the null value */ + + /* test if bin is within range */ + if (ipix < 1 || ipix > histData.haxis1 || pix > histData.maxbin1) + continue; + + if (histData.haxis > 1) + { + if (col2[ii] == FLOATNULLVALUE) + continue; + + axisbin = (col2[ii] - histData.amin2) / histData.binsize2; + iaxisbin = axisbin; + + if (axisbin < 0. || iaxisbin >= histData.haxis2 || axisbin > histData.maxbin2) + continue; + + ipix += (iaxisbin * incr2); + + if (histData.haxis > 2) + { + if (col3[ii] == FLOATNULLVALUE) + continue; + + axisbin = (col3[ii] - histData.amin3) / histData.binsize3; + iaxisbin = axisbin; + if (axisbin < 0. || iaxisbin >= histData.haxis3 || axisbin > histData.maxbin3) + continue; + + ipix += (iaxisbin * incr3); + + if (histData.haxis > 3) + { + if (col4[ii] == FLOATNULLVALUE) + continue; + + axisbin = (col4[ii] - histData.amin4) / histData.binsize4; + iaxisbin = axisbin; + if (axisbin < 0. || iaxisbin >= histData.haxis4 || axisbin > histData.maxbin4) + continue; + + ipix += (iaxisbin * incr4); + + } /* end of haxis > 3 case */ + } /* end of haxis > 2 case */ + } /* end of haxis > 1 case */ + + /* increment the histogram pixel */ + if (histData.weight != FLOATNULLVALUE) /* constant weight factor */ + { + if (histData.himagetype == TINT) + histData.hist.j[ipix] += histData.weight; + else if (histData.himagetype == TSHORT) + histData.hist.i[ipix] += histData.weight; + else if (histData.himagetype == TFLOAT) + histData.hist.r[ipix] += histData.weight; + else if (histData.himagetype == TDOUBLE) + histData.hist.d[ipix] += histData.weight; + else if (histData.himagetype == TBYTE) + histData.hist.b[ipix] += histData.weight; + } + else if (histData.wtrecip) /* use reciprocal of the weight */ + { + if (histData.himagetype == TINT) + histData.hist.j[ipix] += 1./wtcol[ii]; + else if (histData.himagetype == TSHORT) + histData.hist.i[ipix] += 1./wtcol[ii]; + else if (histData.himagetype == TFLOAT) + histData.hist.r[ipix] += 1./wtcol[ii]; + else if (histData.himagetype == TDOUBLE) + histData.hist.d[ipix] += 1./wtcol[ii]; + else if (histData.himagetype == TBYTE) + histData.hist.b[ipix] += 1./wtcol[ii]; + } + else /* no weights */ + { + if (histData.himagetype == TINT) + histData.hist.j[ipix] += wtcol[ii]; + else if (histData.himagetype == TSHORT) + histData.hist.i[ipix] += wtcol[ii]; + else if (histData.himagetype == TFLOAT) + histData.hist.r[ipix] += wtcol[ii]; + else if (histData.himagetype == TDOUBLE) + histData.hist.d[ipix] += wtcol[ii]; + else if (histData.himagetype == TBYTE) + histData.hist.b[ipix] += wtcol[ii]; + } + + } /* end of main loop over all rows */ + + return(0); +} + |