Index: branches/czw_branch/cleanup/psLib/src/fits/psFits.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/fits/psFits.c (revision 24951) +++ branches/czw_branch/cleanup/psLib/src/fits/psFits.c (revision 25051) @@ -150,8 +150,8 @@ char errorBuf[MAX_STRING_LENGTH], *errorMsg; #if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) - errorMsg = strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH); -#else strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH); errorMsg = errorBuf; +#else + errorMsg = strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH); #endif psError(PS_ERR_IO, true, "Directory (%s) for requested file is not accessible: %s", @@ -209,11 +209,12 @@ char errorBuf[64], *errorMsg; # if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) - errorMsg = strerror_r (errno, errorBuf, 64); -# else strerror_r (errno, errorBuf, 64); errorMsg = errorBuf; +# else + errorMsg = strerror_r (errno, errorBuf, 64); # endif - psError (PS_ERR_IO, true, "Failed to delete a previously-existing file: %s", errorMsg); - fprintf (stderr, "errno: %d, %s, %s : %lx, %lx\n", thisErrno, errorMsg, errorBuf, (long int) errorMsg, (long int) errorBuf); + psError(PS_ERR_IO, true, "Failed to delete a previously-existing file (%s), error %d: %s", + name, thisErrno, errorMsg); + //fprintf(stderr, "errno: %d, %s, %s : %lx, %lx\n", thisErrno, errorMsg, errorBuf, (long int) errorMsg, (long int) errorBuf); return NULL; } Index: branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.c (revision 24951) +++ branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.c (revision 25051) @@ -147,5 +147,4 @@ return keyword; } - bool psFitsCheckCompressedImagePHU(const psFits *fits, psMetadata *header) @@ -497,11 +496,14 @@ static bool fitsWriteHeader(psFits *fits, // The FITS file handle const psMetadata *output, // Metadata that is to be output into the FITS file - bool keyStarts // Write out the key starts? + bool keyStarts, // Write out the key starts? + bool phuImage // Are we writing a PHU image? ) { int status = 0; // Status of cfitsio calls + int extnum = psFitsGetExtNum(fits); // Number of extension bool simple = true; // If SIMPLE is T, then the file should conform to the FITS standard psFitsCompressionType compress = psFitsCompressionGetType(fits); // Compression type - if (psFitsGetExtNum(fits) == 0) { + if (extnum == 0) { + // We allow the user to write SIMPLE, but it must be boolean psMetadataItem *simpleItem = psMetadataLookup(output, "SIMPLE"); // SIMPLE in the header @@ -526,5 +528,4 @@ bool compressing = ((!fits->options || fits->options->conventions.compression) && compress != PS_FITS_COMPRESS_NONE) ? true : false; // Are we compressing? - if (compressing) { psMetadataItem *simpleItem = psMetadataLookup(output, "SIMPLE"); // SIMPLE in the header @@ -538,7 +539,14 @@ } } - // int value = simple; // Temporary holder for boolean - // fits_update_key(fits->fd, TLOGICAL, "ZSIMPLE", &value, - // "Uncompressed file's conformance to FITS standard", &status); + if (simple && phuImage && extnum == 1) { + // ZSIMPLE is required for decompression with funpack, etc. For funpack to work, ZSIMPLE needs to + // go early in the FITS header (otherwise we get "Extension doesn't start with SIMPLE or XTENSION + // keyword."). We put it after ZIMAGE by reading ZIMAGE (which sets the insertion pointer) and + // then inserting ZSIMPLE. + char comment[FLEN_CARD]; // Comment for ZIMAGE; unused + int value; // Value for ZIMAGE; unused + fits_read_key(fits->fd, TLOGICAL, "ZIMAGE", &value, comment, &status); + fits_insert_key_log(fits->fd, "ZSIMPLE", simple, "Uncompressed file's conforms to FITS", &status); + } } @@ -699,5 +707,17 @@ PS_ASSERT_METADATA_NON_NULL(output, false); - return fitsWriteHeader(fits, output, true); + return fitsWriteHeader(fits, output, true, false); +} + +bool psFitsWriteHeaderImage(psFits *fits, + const psMetadata *output, + bool phuImage + ) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + PS_ASSERT_FITS_WRITABLE(fits, false); + PS_ASSERT_METADATA_NON_NULL(output, false); + + return fitsWriteHeader(fits, output, true, phuImage); } @@ -733,5 +753,5 @@ } - if (output && !fitsWriteHeader(fits, output, false)) { + if (output && !fitsWriteHeader(fits, output, false, false)) { psError(PS_ERR_IO, false, "Unable to write FITS header.\n"); return false; Index: branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.h =================================================================== --- branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.h (revision 24951) +++ branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.h (revision 25051) @@ -66,4 +66,14 @@ ); +/** Write a header for an image + * Principal difference with psFitsWriteHeader is this allows writing ZSIMPLE for uncompressing a single image + */ +bool psFitsWriteHeaderImage( + psFits *fits, ///< FITS file + const psMetadata *output, ///< Header to output + bool phuImage ///< Is this image supposed to be in the PHU? + ); + + /** Writes a "blank" --- a header only, with no image or table. * Index: branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c (revision 24951) +++ branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c (revision 25051) @@ -616,4 +616,6 @@ } + bool createPHU = false; // Are we creating a PHU? + // Create the image HDU int hdus = psFitsGetSize(fits); // Number of HDUs in file @@ -621,4 +623,5 @@ // We're creating the first image fits_create_img(fits->fd, bitPix, naxis, naxes, &status); + createPHU = true; } else { if (!after) { @@ -627,4 +630,5 @@ // Set status to signal fits_insert_img to insert a new primary HDU status = PREPEND_PRIMARY; + createPHU = true; } else { // Move back one to perform an insert after the previous HDU @@ -641,5 +645,5 @@ // write the header, if any. - if (header && !psFitsWriteHeader(fits, header)) { + if (header && !psFitsWriteHeaderImage(fits, header, createPHU)) { psError(PS_ERR_IO, false, "Unable to write FITS header.\n"); success = false; Index: branches/czw_branch/cleanup/psLib/src/math/psStats.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/math/psStats.c (revision 24951) +++ branches/czw_branch/cleanup/psLib/src/math/psStats.c (revision 25051) @@ -129,15 +129,15 @@ # define COUNT_WARNING(LIMIT, INTERVAL, ...) { \ - static int nCalls = 1; \ - if (nCalls < LIMIT) { \ - psWarning(__VA_ARGS__); \ - } \ - if (!(nCalls % INTERVAL)) { \ - psWarning(__VA_ARGS__); \ + static int nCalls = 1; \ + if (nCalls < LIMIT) { \ + psWarning(__VA_ARGS__); \ + } \ + if (!(nCalls % INTERVAL)) { \ + psWarning(__VA_ARGS__); \ psWarning("(warning raised %d times)", nCalls); \ - } \ - nCalls ++; \ + } \ + nCalls ++; \ } - + /*****************************************************************************/ @@ -217,6 +217,6 @@ for (long i = 0; i < numData; i++) { // Check if the data is with the specified range - if (!isfinite(data[i])) - continue; + if (!isfinite(data[i])) + continue; if (useRange && (data[i] < stats->min)) continue; @@ -242,5 +242,5 @@ if (!isnan(mean)) { - stats->results |= PS_STAT_SAMPLE_MEAN; + stats->results |= PS_STAT_SAMPLE_MEAN; } return true; @@ -278,6 +278,6 @@ for (long i = 0; i < num; i++) { // Check if the data is with the specified range - if (!isfinite(vector[i])) - continue; + if (!isfinite(vector[i])) + continue; if (useRange && (vector[i] < stats->min)) continue; @@ -329,6 +329,6 @@ // into the temporary vectors. for (long i = 0; i < inVector->n; i++) { - if (!isfinite(input[i])) - continue; + if (!isfinite(input[i])) + continue; if (useRange && (input[i] < stats->min)) continue; @@ -353,5 +353,5 @@ // Sort the temporary vector. if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK) - // an error in psVectorSort is a serious error + // an error in psVectorSort is a serious error psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data.")); stats->sampleUQ = NAN; @@ -407,5 +407,5 @@ // If the mean is NAN, then generate a warning and set the stdev to NAN. if (isnan(stats->sampleMean)) { - COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN."); + COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN."); stats->sampleStdev = NAN; return true; @@ -425,6 +425,6 @@ for (long i = 0; i < myVector->n; i++) { // Check if the data is with the specified range - if (!isfinite(data[i])) - continue; + if (!isfinite(data[i])) + continue; if (useRange && (data[i] < stats->min)) { continue; @@ -502,6 +502,6 @@ for (long i = 0; i < myVector->n; i++) { // Check if the data is with the specified range - if (!isfinite(data[i])) - continue; + if (!isfinite(data[i])) + continue; if (useRange && (data[i] < stats->min)) { continue; @@ -771,5 +771,5 @@ stats->robustLQ = min; stats->robustN50 = numValid; - // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits: + // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits: stats->results |= PS_STAT_ROBUST_MEDIAN; stats->results |= PS_STAT_ROBUST_STDEV; @@ -801,12 +801,12 @@ // XXXXX we need to consider this step if errors -> variance if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - // if psVectorHistogram returns false, we have a programming error + // if psVectorHistogram returns false, we have a programming error psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n"); - psFree(histogram); - psFree(cumulative); - psFree(statsMinMax); - psFree(mask); - - return false; + psFree(histogram); + psFree(cumulative); + psFree(statsMinMax); + psFree(mask); + + return false; } if (psTraceGetLevel("psLib.math") >= 8) { @@ -838,9 +838,9 @@ // ADD step 3: Interpolate to the exact 50% position in bin units - stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); - // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); - // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian); - - // convert bin to bin value: this is the robust histogram median. + stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); + // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); + // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian); + + // convert bin to bin value: this is the robust histogram median. if (isnan(stats->robustMedian)) { COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n"); @@ -1054,6 +1054,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1107,5 +1107,5 @@ psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - // if psVectorHistogram returns false, we have a programming error + // if psVectorHistogram returns false, we have a programming error psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n"); psFree(histogram); @@ -1172,6 +1172,6 @@ PS_VECTOR_PRINT_F32(y); } - - // psMinimizeLMChi2 can return false for bad data as well as for serious failures + + // psMinimizeLMChi2 can return false for bad data as well as for serious failures if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); @@ -1235,6 +1235,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1355,5 +1355,5 @@ // psVectorInit (fitMask, 0); - // XXX not sure if these should result in errors or not... + // XXX not sure if these should result in errors or not... if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) { psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); @@ -1432,6 +1432,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1729,6 +1729,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1766,5 +1766,5 @@ COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - goto escape; + goto escape; } @@ -1774,6 +1774,6 @@ stats->fittedMean = min; stats->fittedStdev = 0.0; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN_V4; + stats->results |= PS_STAT_FITTED_STDEV_V4; return true; } @@ -1793,5 +1793,5 @@ psFree(histogram); psFree(statsMinMax); - goto escape; + goto escape; } if (psTraceGetLevel("psLib.math") >= 8) { @@ -1825,5 +1825,5 @@ COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - goto escape; + goto escape; } @@ -1886,6 +1886,6 @@ // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 - // XXX this fit may fail with an error for an ill-conditioned matrix (bad data) - // we probably should be able to handle the data errors gracefully + // XXX this fit may fail with an error for an ill-conditioned matrix (bad data) + // we probably should be able to handle the data errors gracefully psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); @@ -1898,5 +1898,5 @@ psFree(histogram); psFree(statsMinMax); - goto escape; + goto escape; } @@ -1917,5 +1917,5 @@ COUNT_WARNING(10, 100, "fit did not converge\n"); - goto escape; + goto escape; } @@ -1987,5 +1987,5 @@ psFree(histogram); psFree(statsMinMax); - goto escape; + goto escape; } @@ -2898,7 +2898,7 @@ *****************************************************************************/ static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal + psVector *yVec, + psS32 binNum, + psF32 yVal ) { @@ -2982,27 +2982,35 @@ } - // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1 - assert (binValue >= binNum - 1); - assert (binValue <= binNum + 1); - - int fitBin = binValue; - float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin]; - float dY = binValue - fitBin; - tmpFloat = xVec->data.F32[fitBin] + dY * dX; + // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1 + assert (binValue >= binNum - 1); + assert (binValue <= binNum + 1); + + int fitBin = binValue; + float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin]; + float dY = binValue - fitBin; + tmpFloat = xVec->data.F32[fitBin] + dY * dX; } else { // These are special cases where the bin is at the beginning or end of the vector. if (binNum == 0) { // We have two points only at the beginning of the vectors x and y. - // X = (dX/dY)(Y - Yo) + Xo - float dX = xVec->data.F32[1] - xVec->data.F32[0]; - float dY = yVec->data.F32[1] - yVec->data.F32[0]; - tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; + // X = (dX/dY)(Y - Yo) + Xo + float dX = xVec->data.F32[1] - xVec->data.F32[0]; + float dY = yVec->data.F32[1] - yVec->data.F32[0]; + if (dY == 0.0) { + tmpFloat = xVec->data.F32[0]; + } else { + tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; + } } else if (binNum == (xVec->n - 1)) { // We have two points only at the end of the vectors x and y. - // X = (dX/dY)(Y - Yo) + Xo - float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1]; - float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1]; - tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1]; - } + // X = (dX/dY)(Y - Yo) + Xo + float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1]; + float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1]; + if (dY == 0.0) { + tmpFloat = xVec->data.F32[binNum-1]; + } else { + tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1]; + } + } }