IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28795


Ignore:
Timestamp:
Jul 30, 2010, 9:37:39 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/eam_branches/ipp-20100621/psLib/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsImage.c

    r27313 r28795  
    415415    if (fits_read_subset(fits->fd, info->fitsDatatype, info->firstPixel, info->lastPixel,
    416416                         info->increment, nullValue, output->data.V[0], &anynull, &status) != 0) {
    417         psFitsError(status, true, "Reading FITS file failed.");
     417        psFitsError(status, true, "Reading FITS file %s failed.", fits->fd->Fptr->filename);
    418418        return false;
    419419    }

  • branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsTable.c

    r28216 r28795  
    3939
    4040    int status = 0;                     // CFITSIO status
     41
     42    // Check for empty table, which looks like an image
     43    int hdutype;                        // Type of HDU
     44    fits_get_hdu_type(fits->fd, &hdutype, &status);
     45    if (psFitsError(status, true, "Could not determine the HDU type.")) {
     46        return -1;
     47    }
     48    if (hdutype == IMAGE_HDU) {
     49        // It could be an empty table
     50        int naxis = 0;                  // Dimensions of image
     51        if (fits_get_img_dim(fits->fd, &naxis, &status) != 0) {
     52            psFitsError(status, true, "Unable to determine dimension for table.");
     53            return -1;
     54        }
     55        if (naxis != 0) {
     56            psFitsError(PS_ERR_BAD_FITS, true, "Current FITS HDU is not a table.");
     57            return -1;
     58        }
     59        return 0;
     60    }
     61
    4162    long numRows;                       // Number of rows
    4263    if (fits_get_num_rows(fits->fd, &numRows, &status)) {
    … …  
    4869}
    4970
     71// Check if the FITS file is an empty table
     72static bool emptyTableCheck(const psFits *fits // FITS file
     73                            )
     74{
     75    PS_ASSERT_FITS_NON_NULL(fits, false);
     76
     77    if (psFitsGetExtNum(fits) == 0 && !psFitsMoveExtNum(fits, 1, false)) {
     78        psError(PS_ERR_IO, false, "Unable to move to first extension to read table.");
     79        return false;
     80    }
     81
     82    int status = 0;                     // CFITSIO status
     83    int hdutype;                        // Type of HDU
     84    fits_get_hdu_type(fits->fd, &hdutype, &status);
     85    if (status) {
     86        psFitsError(status, true, "Could not determine the HDU type.");
     87        return false;
     88    }
     89    if (hdutype == IMAGE_HDU) {
     90        // It could be an empty table
     91        int naxis = 0;                  // Dimensions of image
     92        if (fits_get_img_dim(fits->fd, &naxis, &status) != 0) {
     93            psFitsError(status, true, "Unable to determine dimension for table.");
     94            return false;
     95        }
     96        if (naxis != 0) {
     97            psFitsError(PS_ERR_BAD_FITS, true, "Current FITS HDU is not a table.");
     98            return false;
     99        }
     100        return true;
     101    }
     102
     103    return false;
     104}
     105
    50106// Check the FITS file in preparation for reading a table
    51107static bool readTableCheck(const psFits *fits // FITS file
    52108                           )
    53109{
    54     PS_ASSERT_FITS_NON_NULL(fits, NULL);
     110    PS_ASSERT_FITS_NON_NULL(fits, false);
    55111
    56112    if (psFitsGetExtNum(fits) == 0 && !psFitsMoveExtNum(fits, 1, false)) {
    … …  
    58114        return false;
    59115    }
    60 
    61116
    62117    // check that we are positioned on a table HDU
    … …  
    64119    int hdutype;                        // Type of HDU
    65120    fits_get_hdu_type(fits->fd, &hdutype, &status);
    66     if (psFitsError(status, true, "Could not determine the HDU type.")) {
     121    if (status) {
     122        psFitsError(status, true, "Could not determine the HDU type.");
    67123        return false;
    68124    }
    69125    if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
    70         psError(PS_ERR_IO, true, _("Current FITS HDU is not a table."));
    71126        return false;
    72127    }
    … …  
    82137
    83138    if (!readTableCheck(fits)) {
     139        if (emptyTableCheck(fits)) {
     140            return psMetadataAlloc();
     141        }
     142        psError(PS_ERR_BAD_FITS, true, "Extension is not a table");
    84143        return NULL;
    85144    }
    … …  
    91150    fits_get_num_rows(fits->fd, &numRows, &status);
    92151    fits_get_num_cols(fits->fd, &numCols, &status);
    93     if (status != 0) {
     152    if (status) {
    94153        psFitsError(status, true, "Failed to determine the size of the current HDU table.");
    95154        return NULL;
    … …  
    204263
    205264    if (!readTableCheck(fits)) {
     265        if (emptyTableCheck(fits)) {
     266            return psArrayAlloc(0);
     267        }
     268        psError(PS_ERR_BAD_FITS, true, "Extension is not a table");
    206269        return NULL;
    207270    }
    … …  
    251314
    252315    if (!readTableCheck(fits)) {
     316        psError(PS_ERR_BAD_FITS, true, "Extension is not a table");
    253317        return NULL;
    254318    }
    … …  
    295359
    296360    if (!readTableCheck(fits)) {
     361        if (emptyTableCheck(fits)) {
     362            return psArrayAlloc(0);
     363        }
     364        psError(PS_ERR_BAD_FITS, true, "Extension is not a table");
    297365        return NULL;
    298366    }
    … …  
    704772
    705773    if (!readTableCheck(fits)) {
     774        psError(PS_ERR_BAD_FITS, true, "Extension is not a table");
    706775        return NULL;
    707776    }
    … …  
    753822
    754823    if (!readTableCheck(fits)) {
     824        if (emptyTableCheck(fits)) {
     825            return psMetadataAlloc();
     826        }
     827        psError(PS_ERR_BAD_FITS, true, "Extension is not a table");
    755828        return NULL;
    756829    }
    … …  
    762835    fits_get_num_rows(fits->fd, &numRows, &status);
    763836    fits_get_num_cols(fits->fd, &numCols, &status);
    764     if (psFitsError(status, true, "Failed to determine the size of the current HDU table.")) {
     837    if (status) {
     838        psFitsError(status, true, "Failed to determine the size of the current HDU table.");
    765839        return NULL;
    766840    }
    … …  
    768842    int hdutype;                        // Type of HDU: need to distinguish ASCII and binary tables
    769843    fits_get_hdu_type(fits->fd, &hdutype, &status);
    770     if (psFitsError(status, true, "Could not determine the HDU type.")) {
     844    if (status) {
     845        psFitsError(status, true, "Could not determine the HDU type.");
    771846        return false;
    772847    }
    … …  
    786861        long repeat;                     // Number of repeats
    787862        fits_get_eqcoltype(fits->fd, col, &cfitsioType, &repeat, NULL, &status);
    788         if (psFitsError(status, true, "Could not determine the column data for %s", name)) {
     863        if (status) {
     864            psFitsError(status, true, "Could not determine the column data for %s", name);
    789865            psFree(table);
    790866            return false;
    … …  
    805881            }
    806882            fits_read_col_str(fits->fd, col, 1, 1, numRows, "", (char**)array->data, NULL, &status);
    807             if (psFitsError(status, true, "Failed to read column %s", name)) {
     883            if (status) {
     884                psFitsError(status, true, "Failed to read column %s", name);
    808885                psFree(array);
    809886                psFree(table);
    … …  
    822899            fits_read_col(fits->fd, cfitsioType, col, 1, 1, numRows, NULL,
    823900                          vector->data.U8, NULL, &status);
    824             if (psFitsError(status, true, "Failed to read column %s", name)) {
     901            if (status) {
     902                psFitsError(status, true, "Failed to read column %s", name);
    825903                psFree(vector);
    826904                psFree(table);
    … …  
    9331011    psFree(names);
    9341012    psFree(types);
    935     if (psFitsError(status, true, "Unable to create FITS table with %d columns and %ld rows",
    936                     numCols, numRows)) {
     1013    if (status) {
     1014        psFitsError(status, true, "Unable to create FITS table with %d columns and %ld rows",
     1015                    numCols, numRows);
    9371016        psFree(columns);
    9381017        return false;
    … …  
    9861065        }
    9871066        // Check error status from writing column
    988         if (psFitsError(status, true, "Unable to write column %ld of FITS table", col)) {
     1067        if (status) {
     1068            psFitsError(status, true, "Unable to write column %ld of FITS table", col);
    9891069            psFree(columns);
    9901070            return false;
    … …  
    9971077    // hurt.
    9981078    ffrdef(fits->fd, &status);
    999     if (psFitsError(status, true, "Could not re-scan HDU.")) {
     1079    if (status) {
     1080        psFitsError(status, true, "Could not re-scan HDU.");
    10001081        return false;
    10011082    }

  • branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageCovariance.c

    r28405 r28795  
    3434static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix
    3535                                      const psKernel *kernel, // Convolution kernel
    36                                       int x, int y            // Coordinates in output covariance matrix
     36                                      int x, int y,           // Coordinates in output covariance matrix
     37                                      float scale             // Scale to apply
    3738                                      )
    3839{
    … …  
    8384    }
    8485
    85     return sum;
     86    return scale * sum;
    8687}
    8788
    … …  
    9192    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
    9293    psAssert(job->args, "No job arguments");
    93     psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n);
     94    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
    9495
    9596    psKernel *out = job->args->data[0]; // Output covariance matrix
    … …  
    9899    int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix
    99100    int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix
    100 
    101     out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
     101    float scale = PS_SCALAR_VALUE(job->args->data[5], F32); // Scaling to apply
     102
     103    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
    102104
    103105    return true;
    … …  
    127129
    128130    // Check for non-finite elements
     131    double sumKernel = 0.0, sumKernel2 = 0.0; // Sum of the kernel
    129132    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    130133        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    … …  
    135138                return NULL;
    136139            }
     140            sumKernel += kernel->kernel[y][x];
     141            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
    137142        }
    138143    }
    … …  
    155160    int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
    156161    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
     162    float scale = 1.0 / sumKernel2;          // Scaling to apply
    157163
    158164    for (int y = yMin; y <= yMax; y++) {
    … …  
    165171                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
    166172                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
     173                PS_ARRAY_ADD_SCALAR(job->args, scale, PS_TYPE_F32);
    167174                if (!psThreadJobAddPending(job)) {
    168175                    psFree(covar);
    … …  
    170177                }
    171178            } else {
    172                 out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
    173             }
    174         }
    175     }
    176     psFree(covar);
     179                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
     180            }
     181        }
     182    }
    177183
    178184    if (threaded && !psThreadPoolWait(true)) {
    … …  
    180186        return false;
    181187    }
     188
     189    psFree(covar);
    182190
    183191    return out;
    … …  
    194202
    195203    // Check for non-finite elements
     204    double sumKernel2 = 0.0; // Sum of the squared kernel
    196205    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    197206        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    … …  
    202211                return NAN;
    203212            }
     213            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
    204214        }
    205215    }
    … …  
    215225    }
    216226
    217     float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
     227    float scale = 1.0 / sumKernel2;     // Scale to apply
     228    float factor = imageCovarianceCalculate(covar, kernel, 0, 0, scale); // Covariance factor
    218229    psFree(covar);
    219230    return factor;
    … …  
    338349        }
    339350    }
    340     psFree(covar);
    341 
    342351    if (threaded && !psThreadPoolWait(true)) {
    343352        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
    344353        return false;
    345354    }
     355    psFree(covar);
    346356
    347357    return out;
    … …  
    616626    if (set && !threaded) {
    617627        {
    618             psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
     628            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 6);
    619629            task->function = &imageCovarianceCalculateThread;
    620630            psThreadTaskAdd(task);

  • branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageInterpolate.c

    r26892 r28795  
    427427
    428428// Determine the result of the interpolation after all the math has been done
    429 #define INTERPOLATE_RESULT() \
    430     psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; /* Status of interpolation */ \
    431     *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage; \
    432     if (wantVariance) { \
    433         *varianceValue = sumVariance / (sumKernel2 - sumBad); \
    434     } \
    435     if (sumKernel == 0.0) { \
    436         /* No kernel contributions */ \
    437         if (haveMask && maskValue) { \
    438             *maskValue |= interp->badMask; \
    439         } \
    440         status = PS_INTERPOLATE_STATUS_BAD; \
    441     } else if (sumBad == 0) { \
    442         /* Completely good pixel */ \
    443         status = PS_INTERPOLATE_STATUS_GOOD; \
    444     } else if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { \
    445         /* Some pixels masked: poor pixel */ \
    446         if (haveMask && maskValue) { \
    447             *maskValue |= interp->poorMask; \
    448         } \
    449         status = PS_INTERPOLATE_STATUS_POOR; \
    450     } else { \
    451         /* Many pixels (or a few important pixels) masked: bad pixel */ \
    452         if (haveMask && maskValue) { \
    453             *maskValue |= interp->badMask; \
    454         } \
    455         status = PS_INTERPOLATE_STATUS_BAD; \
    456     }
     429static psImageInterpolateStatus interpolateResult(const psImageInterpolation *interp,
     430                                                  double *imageValue, double *varianceValue,
     431                                                  psImageMaskType *maskValue,
     432                                                  double sumImage, double sumVariance, double sumBad,
     433                                                  double sumKernel, double sumKernel2,
     434                                                  bool wantVariance, bool haveMask)
     435{
     436    *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage;
     437    if (wantVariance) {
     438        if (sumBad > 0) {
     439            sumVariance *= sumKernel2 / (sumKernel2 - sumBad);
     440        }
     441        *varianceValue = sumVariance / PS_SQR(sumKernel);
     442    }
     443    if (sumKernel == 0.0) {
     444        // No kernel contributions at all
     445        if (haveMask && maskValue) {
     446            *maskValue |= interp->badMask;
     447        }
     448        return PS_INTERPOLATE_STATUS_BAD;
     449    }
     450    if (sumBad == 0) {
     451        // Completely good pixel
     452        return PS_INTERPOLATE_STATUS_GOOD;
     453    }
     454    if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) {
     455        // Some pixels masked: poor pixel
     456        if (haveMask && maskValue) {
     457            *maskValue |= interp->poorMask;
     458        }
     459        return PS_INTERPOLATE_STATUS_POOR;
     460    }
     461    // Many pixels (or a few important pixels) masked: bad pixel
     462    if (haveMask && maskValue) {
     463        *maskValue |= interp->badMask;
     464    }
     465    return PS_INTERPOLATE_STATUS_BAD;
     466}
    457467
    458468// Interpolation engine for separable interpolation kernels
    … …  
    703713    }
    704714
    705     INTERPOLATE_RESULT();
    706 
    707715    psFree(xKernelNew);
    708716    psFree(yKernelNew);
    … …  
    710718    psFree(yKernel2New);
    711719
    712     return status;
     720    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
     721                             sumKernel, sumKernel2, wantVariance, haveMask);
    713722}
    714723
    … …  
    861870    }
    862871
    863     INTERPOLATE_RESULT();
    864 
    865     return status;
     872    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
     873                             sumKernel, sumKernel2, wantVariance, haveMask);
    866874}
    867875

  • branches/eam_branches/ipp-20100621/psLib/src/math/psMinimizePowell.c

    r28655 r28795  
    8383f(param + b * line) < f(param + c * line)
    8484a < b < c
    85  
     85
    8686Algorithm:
    87  
     87
    8888XXX completely ad hoc:
    8989start with the user-supplied starting parameter and
    9090call that b.  Calculate a/c as a fractional amount smaller/larger than b.
    9191Repeat this process until a local minimum is found.
    92  
     92
    9393XXX: new algorithm:
    9494start at x=0, expand in one direction until the function
    … …  
    9696increases, or x is too large.  If thst does not work, expand in the other
    9797direction.
    98  
     98
    9999XXX: output bracket vector should be an input as well.
    100100*****************************************************************************/
    … …  
    335335along that vector and returns the offset along that vector at which the
    336336minimum is determined.
    337  
     337
    338338XXX: This routine is not very efficient in terms of total evaluations of the
    339339function.
    … …  
    491491points at which the function is varied are in the argument "coords" which is
    492492a psArray of psVectors: each vector represents a different coordinate.
    493  
     493
    494494XXX: We do not use Brent's method.
    495495 *****************************************************************************/
    … …  
    709709This routine is to be used with the psMinimizeChi2Powell() function below.
    710710and the psMinimizePowell() function above.
    711  
     711
    712712The basic idea is calculate chi-squared for a set of params/coords/errors.
    713713This functions uses global variables to receive the function pointer, the
    … …  
    764764This routine must minimize the chi-squared match of a set of data points and
    765765values for a possibly multi-dimensional function.
    766  
     766
    767767The basic idea is to use the psMinimizePowell() function defined above.  In
    768768order to do so, we defined above a function myPowellChi2Func() which takes
    … …  
    780780    psMinimizeChi2PowellFunc model)
    781781{
     782    PS_ASSERT_VECTOR_NON_NULL(params, false);
     783    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
     784    PS_ASSERT_ARRAY_NON_NULL(coords, false);
     785    PS_ASSERT_VECTOR_NON_NULL(value, false);
     786    PS_ASSERT_VECTOR_TYPE(value, PS_TYPE_F32, false);
     787    PS_ASSERT_VECTORS_SIZE_EQUAL(coords, value, false);
     788
    782789    // Generate extended version of coords array, so we can pass in extra data to the chi^2 function
    783790    psArray *newCoords = psArrayAlloc(coords->n + 3);
    … …  
    791798    newCoords->data[coords->n + 2] = &model;
    792799
    793     bool success = psMinimizePowell(min, params, constraint->paramMask, newCoords, myPowellChi2Func);
     800    bool success = psMinimizePowell(min, params, constraint ? constraint->paramMask : NULL,
     801                                    newCoords, myPowellChi2Func);
    794802
    795803    newCoords->data[coords->n - 1] = NULL; // We can't free the array with a function pointer on it
Note: See TracChangeset for help on using the changeset viewer.