IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 11, 2009, 3:54:00 PM (17 years ago)
Author:
watersc1
Message:

Bringing branch up-to-date at revision 25049

Location:
branches/czw_branch/cleanup
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/cleanup

  • branches/czw_branch/cleanup/psLib/src/fits/psFits.c

    r24285 r25051  
    150150            char errorBuf[MAX_STRING_LENGTH], *errorMsg;
    151151#if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE)
    152             errorMsg = strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH);
    153 #else
    154152            strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH);
    155153            errorMsg = errorBuf;
     154#else
     155            errorMsg = strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH);
    156156#endif
    157157            psError(PS_ERR_IO, true, "Directory (%s) for requested file is not accessible: %s",
    … …  
    209209                char errorBuf[64], *errorMsg;
    210210# if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE)
    211                 errorMsg = strerror_r (errno, errorBuf, 64);
    212 # else
    213211                strerror_r (errno, errorBuf, 64);
    214212                errorMsg = errorBuf;
     213# else
     214                errorMsg = strerror_r (errno, errorBuf, 64);
    215215# endif
    216                 psError (PS_ERR_IO, true, "Failed to delete a previously-existing file: %s", errorMsg);
    217                 fprintf (stderr, "errno: %d, %s, %s : %lx, %lx\n", thisErrno, errorMsg, errorBuf, (long int) errorMsg, (long int) errorBuf);
     216                psError(PS_ERR_IO, true, "Failed to delete a previously-existing file (%s), error %d: %s",
     217                        name, thisErrno, errorMsg);
     218                //fprintf(stderr, "errno: %d, %s, %s : %lx, %lx\n", thisErrno, errorMsg, errorBuf, (long int) errorMsg, (long int) errorBuf);
    218219                return NULL;
    219220            }

  • branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.c

    r23376 r25051  
    147147    return keyword;
    148148}
    149 
    150149
    151150bool psFitsCheckCompressedImagePHU(const psFits *fits, psMetadata *header)
    … …  
    497496static bool fitsWriteHeader(psFits *fits, // The FITS file handle
    498497                            const psMetadata *output, // Metadata that is to be output into the FITS file
    499                             bool keyStarts // Write out the key starts?
     498                            bool keyStarts, // Write out the key starts?
     499                            bool phuImage   // Are we writing a PHU image?
    500500                           )
    501501{
    502502    int status = 0;                     // Status of cfitsio calls
     503    int extnum = psFitsGetExtNum(fits); // Number of extension
    503504    bool simple = true;                 // If SIMPLE is T, then the file should conform to the FITS standard
    504505    psFitsCompressionType compress = psFitsCompressionGetType(fits); // Compression type
    505     if (psFitsGetExtNum(fits) == 0) {
     506    if (extnum == 0) {
     507
    506508        // We allow the user to write SIMPLE, but it must be boolean
    507509        psMetadataItem *simpleItem = psMetadataLookup(output, "SIMPLE"); // SIMPLE in the header
    … …  
    526528    bool compressing = ((!fits->options || fits->options->conventions.compression) &&
    527529                        compress != PS_FITS_COMPRESS_NONE) ? true : false; // Are we compressing?
    528 
    529530    if (compressing) {
    530531        psMetadataItem *simpleItem = psMetadataLookup(output, "SIMPLE"); // SIMPLE in the header
    … …  
    538539            }
    539540        }
    540         //        int value = simple;             // Temporary holder for boolean
    541         //        fits_update_key(fits->fd, TLOGICAL, "ZSIMPLE", &value,
    542         //                        "Uncompressed file's conformance to FITS standard", &status);
     541        if (simple && phuImage && extnum == 1) {
     542            // ZSIMPLE is required for decompression with funpack, etc.  For funpack to work, ZSIMPLE needs to
     543            // go early in the FITS header (otherwise we get "Extension doesn't start with SIMPLE or XTENSION
     544            // keyword.").  We put it after ZIMAGE by reading ZIMAGE (which sets the insertion pointer) and
     545            // then inserting ZSIMPLE.
     546            char comment[FLEN_CARD];    // Comment for ZIMAGE; unused
     547            int value;                  // Value for ZIMAGE; unused
     548            fits_read_key(fits->fd, TLOGICAL, "ZIMAGE", &value, comment, &status);
     549            fits_insert_key_log(fits->fd, "ZSIMPLE", simple, "Uncompressed file's conforms to FITS", &status);
     550        }
    543551    }
    544552
    … …  
    699707    PS_ASSERT_METADATA_NON_NULL(output, false);
    700708
    701     return fitsWriteHeader(fits, output, true);
     709    return fitsWriteHeader(fits, output, true, false);
     710}
     711
     712bool psFitsWriteHeaderImage(psFits *fits,
     713                            const psMetadata *output,
     714                            bool phuImage
     715                            )
     716{
     717    PS_ASSERT_FITS_NON_NULL(fits, false);
     718    PS_ASSERT_FITS_WRITABLE(fits, false);
     719    PS_ASSERT_METADATA_NON_NULL(output, false);
     720
     721    return fitsWriteHeader(fits, output, true, phuImage);
    702722}
    703723
    … …  
    733753    }
    734754
    735     if (output && !fitsWriteHeader(fits, output, false)) {
     755    if (output && !fitsWriteHeader(fits, output, false, false)) {
    736756        psError(PS_ERR_IO, false, "Unable to write FITS header.\n");
    737757        return false;

  • branches/czw_branch/cleanup/psLib/src/fits/psFitsHeader.h

    r22729 r25051  
    6666);
    6767
     68/** Write a header for an image
     69 *  Principal difference with psFitsWriteHeader is this allows writing ZSIMPLE for uncompressing a single image
     70 */
     71bool psFitsWriteHeaderImage(
     72    psFits *fits,                       ///< FITS file
     73    const psMetadata *output,           ///< Header to output
     74    bool phuImage                       ///< Is this image supposed to be in the PHU?
     75    );
     76
     77
    6878/** Writes a "blank" --- a header only, with no image or table.
    6979 *

  • branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c

    r24951 r25051  
    616616    }
    617617
     618    bool createPHU = false;             // Are we creating a PHU?
     619
    618620    // Create the image HDU
    619621    int hdus = psFitsGetSize(fits);     // Number of HDUs in file
    … …  
    621623        // We're creating the first image
    622624        fits_create_img(fits->fd, bitPix, naxis, naxes, &status);
     625        createPHU = true;
    623626    } else {
    624627        if (!after) {
    … …  
    627630                // Set status to signal fits_insert_img to insert a new primary HDU
    628631                status = PREPEND_PRIMARY;
     632                createPHU = true;
    629633            } else {
    630634                // Move back one to perform an insert after the previous HDU
    … …  
    641645
    642646    // write the header, if any.
    643     if (header && !psFitsWriteHeader(fits, header)) {
     647    if (header && !psFitsWriteHeaderImage(fits, header, createPHU)) {
    644648        psError(PS_ERR_IO, false, "Unable to write FITS header.\n");
    645649        success = false;

  • branches/czw_branch/cleanup/psLib/src/math/psStats.c

    r24951 r25051  
    129129
    130130# define COUNT_WARNING(LIMIT, INTERVAL, ...) { \
    131         static int nCalls = 1; \
    132         if (nCalls < LIMIT) { \
    133             psWarning(__VA_ARGS__); \
    134         } \
    135         if (!(nCalls % INTERVAL)) { \
    136             psWarning(__VA_ARGS__); \
     131        static int nCalls = 1; \
     132        if (nCalls < LIMIT) { \
     133            psWarning(__VA_ARGS__); \
     134        } \
     135        if (!(nCalls % INTERVAL)) { \
     136            psWarning(__VA_ARGS__); \
    137137            psWarning("(warning raised %d times)", nCalls); \
    138         } \
    139         nCalls ++; \
     138        } \
     139        nCalls ++; \
    140140}
    141  
     141
    142142
    143143/*****************************************************************************/
    … …  
    217217    for (long i = 0; i < numData; i++) {
    218218        // Check if the data is with the specified range
    219         if (!isfinite(data[i]))
    220             continue;
     219        if (!isfinite(data[i]))
     220            continue;
    221221        if (useRange && (data[i] < stats->min))
    222222            continue;
    … …  
    242242
    243243    if (!isnan(mean)) {
    244         stats->results |= PS_STAT_SAMPLE_MEAN;
     244        stats->results |= PS_STAT_SAMPLE_MEAN;
    245245    }
    246246    return true;
    … …  
    278278    for (long i = 0; i < num; i++) {
    279279        // Check if the data is with the specified range
    280         if (!isfinite(vector[i]))
    281             continue;
     280        if (!isfinite(vector[i]))
     281            continue;
    282282        if (useRange && (vector[i] < stats->min))
    283283            continue;
    … …  
    329329    // into the temporary vectors.
    330330    for (long i = 0; i < inVector->n; i++) {
    331         if (!isfinite(input[i]))
    332             continue;
     331        if (!isfinite(input[i]))
     332            continue;
    333333        if (useRange && (input[i] < stats->min))
    334334            continue;
    … …  
    353353    // Sort the temporary vector.
    354354    if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
    355         // an error in psVectorSort is a serious error
     355        // an error in psVectorSort is a serious error
    356356        psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
    357357        stats->sampleUQ = NAN;
    … …  
    407407    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    408408    if (isnan(stats->sampleMean)) {
    409         COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.");
     409        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.");
    410410        stats->sampleStdev = NAN;
    411411        return true;
    … …  
    425425    for (long i = 0; i < myVector->n; i++) {
    426426        // Check if the data is with the specified range
    427         if (!isfinite(data[i]))
    428             continue;
     427        if (!isfinite(data[i]))
     428            continue;
    429429        if (useRange && (data[i] < stats->min)) {
    430430            continue;
    … …  
    502502    for (long i = 0; i < myVector->n; i++) {
    503503        // Check if the data is with the specified range
    504         if (!isfinite(data[i]))
    505             continue;
     504        if (!isfinite(data[i]))
     505            continue;
    506506        if (useRange && (data[i] < stats->min)) {
    507507            continue;
    … …  
    771771            stats->robustLQ = min;
    772772            stats->robustN50 = numValid;
    773             // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits:
     773            // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits:
    774774            stats->results |= PS_STAT_ROBUST_MEDIAN;
    775775            stats->results |= PS_STAT_ROBUST_STDEV;
    … …  
    801801        // XXXXX we need to consider this step if errors -> variance
    802802        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    803             // if psVectorHistogram returns false, we have a programming error
     803            // if psVectorHistogram returns false, we have a programming error
    804804            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n");
    805             psFree(histogram);
    806             psFree(cumulative);
    807             psFree(statsMinMax);
    808             psFree(mask);
    809 
    810             return false;
     805            psFree(histogram);
     806            psFree(cumulative);
     807            psFree(statsMinMax);
     808            psFree(mask);
     809
     810            return false;
    811811        }
    812812        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    838838
    839839        // ADD step 3: Interpolate to the exact 50% position in bin units
    840         stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    841         // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    842         // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
    843 
    844         // convert bin to bin value: this is the robust histogram median.
     840        stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     841        // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     842        // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
     843
     844        // convert bin to bin value: this is the robust histogram median.
    845845        if (isnan(stats->robustMedian)) {
    846846            COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n");
    … …  
    10541054        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
    10551055            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
    1056             return false;
    1057         }
     1056            return false;
     1057        }
    10581058    }
    10591059
    … …  
    11071107        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    11081108        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    1109             // if psVectorHistogram returns false, we have a programming error
     1109            // if psVectorHistogram returns false, we have a programming error
    11101110            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n");
    11111111            psFree(histogram);
    … …  
    11721172            PS_VECTOR_PRINT_F32(y);
    11731173        }
    1174        
    1175         // psMinimizeLMChi2 can return false for bad data as well as for serious failures
     1174
     1175        // psMinimizeLMChi2 can return false for bad data as well as for serious failures
    11761176        if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
    11771177            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    … …  
    12351235        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
    12361236            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
    1237             return false;
    1238         }
     1237            return false;
     1238        }
    12391239    }
    12401240
    … …  
    13551355        // psVectorInit (fitMask, 0);
    13561356
    1357         // XXX not sure if these should result in errors or not...
     1357        // XXX not sure if these should result in errors or not...
    13581358        if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
    13591359            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    … …  
    14321432        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
    14331433            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
    1434             return false;
    1435         }
     1434            return false;
     1435        }
    14361436    }
    14371437
    … …  
    17291729        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
    17301730            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
    1731             return false;
    1732         }
     1731            return false;
     1732        }
    17331733    }
    17341734
    … …  
    17661766            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    17671767            psFree(statsMinMax);
    1768             goto escape;
     1768            goto escape;
    17691769        }
    17701770
    … …  
    17741774            stats->fittedMean = min;
    17751775            stats->fittedStdev = 0.0;
    1776             stats->results |= PS_STAT_FITTED_MEAN_V4;
    1777             stats->results |= PS_STAT_FITTED_STDEV_V4;
     1776            stats->results |= PS_STAT_FITTED_MEAN_V4;
     1777            stats->results |= PS_STAT_FITTED_STDEV_V4;
    17781778            return true;
    17791779        }
    … …  
    17931793            psFree(histogram);
    17941794            psFree(statsMinMax);
    1795             goto escape;
     1795            goto escape;
    17961796        }
    17971797        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    18251825            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    18261826            psFree(statsMinMax);
    1827             goto escape;
     1827            goto escape;
    18281828        }
    18291829
    … …  
    18861886
    18871887            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
    1888             // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
    1889             // we probably should be able to handle the data errors gracefully
     1888            // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
     1889            // we probably should be able to handle the data errors gracefully
    18901890            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    18911891            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
    … …  
    18981898                psFree(histogram);
    18991899                psFree(statsMinMax);
    1900                 goto escape;
     1900                goto escape;
    19011901            }
    19021902
    … …  
    19171917
    19181918                COUNT_WARNING(10, 100, "fit did not converge\n");
    1919                 goto escape;
     1919                goto escape;
    19201920            }
    19211921
    … …  
    19871987                psFree(histogram);
    19881988                psFree(statsMinMax);
    1989                 goto escape;
     1989                goto escape;
    19901990            }
    19911991
    … …  
    28982898*****************************************************************************/
    28992899static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec,
    2900                                                 psVector *yVec,
    2901                                                 psS32 binNum,
    2902                                                 psF32 yVal
     2900                                                psVector *yVec,
     2901                                                psS32 binNum,
     2902                                                psF32 yVal
    29032903    )
    29042904{
    … …  
    29822982        }
    29832983
    2984         // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
    2985         assert (binValue >= binNum - 1);
    2986         assert (binValue <= binNum + 1);
    2987 
    2988         int fitBin = binValue;
    2989         float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
    2990         float dY = binValue - fitBin;
    2991         tmpFloat = xVec->data.F32[fitBin] + dY * dX;
     2984        // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
     2985        assert (binValue >= binNum - 1);
     2986        assert (binValue <= binNum + 1);
     2987
     2988        int fitBin = binValue;
     2989        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
     2990        float dY = binValue - fitBin;
     2991        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
    29922992    } else {
    29932993        // These are special cases where the bin is at the beginning or end of the vector.
    29942994        if (binNum == 0) {
    29952995            // We have two points only at the beginning of the vectors x and y.
    2996             // X = (dX/dY)(Y - Yo) + Xo
    2997             float dX = xVec->data.F32[1] - xVec->data.F32[0];
    2998             float dY = yVec->data.F32[1] - yVec->data.F32[0];       
    2999             tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
     2996            // X = (dX/dY)(Y - Yo) + Xo
     2997            float dX = xVec->data.F32[1] - xVec->data.F32[0];
     2998            float dY = yVec->data.F32[1] - yVec->data.F32[0];
     2999            if (dY == 0.0) {
     3000                tmpFloat = xVec->data.F32[0];
     3001            } else {
     3002                tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
     3003            }
    30003004        } else if (binNum == (xVec->n - 1)) {
    30013005            // We have two points only at the end of the vectors x and y.
    3002             // X = (dX/dY)(Y - Yo) + Xo
    3003             float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
    3004             float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];           
    3005             tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
    3006         }
     3006            // X = (dX/dY)(Y - Yo) + Xo
     3007            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
     3008            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];
     3009            if (dY == 0.0) {
     3010                tmpFloat = xVec->data.F32[binNum-1];
     3011            } else {
     3012                tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
     3013            }
     3014        }
    30073015    }
    30083016
Note: See TracChangeset for help on using the changeset viewer.