Index: /branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsImage.c =================================================================== --- /branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsImage.c (revision 28794) +++ /branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsImage.c (revision 28795) @@ -415,5 +415,5 @@ if (fits_read_subset(fits->fd, info->fitsDatatype, info->firstPixel, info->lastPixel, info->increment, nullValue, output->data.V[0], &anynull, &status) != 0) { - psFitsError(status, true, "Reading FITS file failed."); + psFitsError(status, true, "Reading FITS file %s failed.", fits->fd->Fptr->filename); return false; } Index: /branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsTable.c =================================================================== --- /branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsTable.c (revision 28794) +++ /branches/eam_branches/ipp-20100621/psLib/src/fits/psFitsTable.c (revision 28795) @@ -39,4 +39,25 @@ int status = 0; // CFITSIO status + + // Check for empty table, which looks like an image + int hdutype; // Type of HDU + fits_get_hdu_type(fits->fd, &hdutype, &status); + if (psFitsError(status, true, "Could not determine the HDU type.")) { + return -1; + } + if (hdutype == IMAGE_HDU) { + // It could be an empty table + int naxis = 0; // Dimensions of image + if (fits_get_img_dim(fits->fd, &naxis, &status) != 0) { + psFitsError(status, true, "Unable to determine dimension for table."); + return -1; + } + if (naxis != 0) { + psFitsError(PS_ERR_BAD_FITS, true, "Current FITS HDU is not a table."); + return -1; + } + return 0; + } + long numRows; // Number of rows if (fits_get_num_rows(fits->fd, &numRows, &status)) { @@ -48,9 +69,44 @@ } +// Check if the FITS file is an empty table +static bool emptyTableCheck(const psFits *fits // FITS file + ) +{ + PS_ASSERT_FITS_NON_NULL(fits, false); + + if (psFitsGetExtNum(fits) == 0 && !psFitsMoveExtNum(fits, 1, false)) { + psError(PS_ERR_IO, false, "Unable to move to first extension to read table."); + return false; + } + + int status = 0; // CFITSIO status + int hdutype; // Type of HDU + fits_get_hdu_type(fits->fd, &hdutype, &status); + if (status) { + psFitsError(status, true, "Could not determine the HDU type."); + return false; + } + if (hdutype == IMAGE_HDU) { + // It could be an empty table + int naxis = 0; // Dimensions of image + if (fits_get_img_dim(fits->fd, &naxis, &status) != 0) { + psFitsError(status, true, "Unable to determine dimension for table."); + return false; + } + if (naxis != 0) { + psFitsError(PS_ERR_BAD_FITS, true, "Current FITS HDU is not a table."); + return false; + } + return true; + } + + return false; +} + // Check the FITS file in preparation for reading a table static bool readTableCheck(const psFits *fits // FITS file ) { - PS_ASSERT_FITS_NON_NULL(fits, NULL); + PS_ASSERT_FITS_NON_NULL(fits, false); if (psFitsGetExtNum(fits) == 0 && !psFitsMoveExtNum(fits, 1, false)) { @@ -58,5 +114,4 @@ return false; } - // check that we are positioned on a table HDU @@ -64,9 +119,9 @@ int hdutype; // Type of HDU fits_get_hdu_type(fits->fd, &hdutype, &status); - if (psFitsError(status, true, "Could not determine the HDU type.")) { + if (status) { + psFitsError(status, true, "Could not determine the HDU type."); return false; } if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) { - psError(PS_ERR_IO, true, _("Current FITS HDU is not a table.")); return false; } @@ -82,4 +137,8 @@ if (!readTableCheck(fits)) { + if (emptyTableCheck(fits)) { + return psMetadataAlloc(); + } + psError(PS_ERR_BAD_FITS, true, "Extension is not a table"); return NULL; } @@ -91,5 +150,5 @@ fits_get_num_rows(fits->fd, &numRows, &status); fits_get_num_cols(fits->fd, &numCols, &status); - if (status != 0) { + if (status) { psFitsError(status, true, "Failed to determine the size of the current HDU table."); return NULL; @@ -204,4 +263,8 @@ if (!readTableCheck(fits)) { + if (emptyTableCheck(fits)) { + return psArrayAlloc(0); + } + psError(PS_ERR_BAD_FITS, true, "Extension is not a table"); return NULL; } @@ -251,4 +314,5 @@ if (!readTableCheck(fits)) { + psError(PS_ERR_BAD_FITS, true, "Extension is not a table"); return NULL; } @@ -295,4 +359,8 @@ if (!readTableCheck(fits)) { + if (emptyTableCheck(fits)) { + return psArrayAlloc(0); + } + psError(PS_ERR_BAD_FITS, true, "Extension is not a table"); return NULL; } @@ -704,4 +772,5 @@ if (!readTableCheck(fits)) { + psError(PS_ERR_BAD_FITS, true, "Extension is not a table"); return NULL; } @@ -753,4 +822,8 @@ if (!readTableCheck(fits)) { + if (emptyTableCheck(fits)) { + return psMetadataAlloc(); + } + psError(PS_ERR_BAD_FITS, true, "Extension is not a table"); return NULL; } @@ -762,5 +835,6 @@ fits_get_num_rows(fits->fd, &numRows, &status); fits_get_num_cols(fits->fd, &numCols, &status); - if (psFitsError(status, true, "Failed to determine the size of the current HDU table.")) { + if (status) { + psFitsError(status, true, "Failed to determine the size of the current HDU table."); return NULL; } @@ -768,5 +842,6 @@ int hdutype; // Type of HDU: need to distinguish ASCII and binary tables fits_get_hdu_type(fits->fd, &hdutype, &status); - if (psFitsError(status, true, "Could not determine the HDU type.")) { + if (status) { + psFitsError(status, true, "Could not determine the HDU type."); return false; } @@ -786,5 +861,6 @@ long repeat; // Number of repeats fits_get_eqcoltype(fits->fd, col, &cfitsioType, &repeat, NULL, &status); - if (psFitsError(status, true, "Could not determine the column data for %s", name)) { + if (status) { + psFitsError(status, true, "Could not determine the column data for %s", name); psFree(table); return false; @@ -805,5 +881,6 @@ } fits_read_col_str(fits->fd, col, 1, 1, numRows, "", (char**)array->data, NULL, &status); - if (psFitsError(status, true, "Failed to read column %s", name)) { + if (status) { + psFitsError(status, true, "Failed to read column %s", name); psFree(array); psFree(table); @@ -822,5 +899,6 @@ fits_read_col(fits->fd, cfitsioType, col, 1, 1, numRows, NULL, vector->data.U8, NULL, &status); - if (psFitsError(status, true, "Failed to read column %s", name)) { + if (status) { + psFitsError(status, true, "Failed to read column %s", name); psFree(vector); psFree(table); @@ -933,6 +1011,7 @@ psFree(names); psFree(types); - if (psFitsError(status, true, "Unable to create FITS table with %d columns and %ld rows", - numCols, numRows)) { + if (status) { + psFitsError(status, true, "Unable to create FITS table with %d columns and %ld rows", + numCols, numRows); psFree(columns); return false; @@ -986,5 +1065,6 @@ } // Check error status from writing column - if (psFitsError(status, true, "Unable to write column %ld of FITS table", col)) { + if (status) { + psFitsError(status, true, "Unable to write column %ld of FITS table", col); psFree(columns); return false; @@ -997,5 +1077,6 @@ // hurt. ffrdef(fits->fd, &status); - if (psFitsError(status, true, "Could not re-scan HDU.")) { + if (status) { + psFitsError(status, true, "Could not re-scan HDU."); return false; } Index: /branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageCovariance.c =================================================================== --- /branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageCovariance.c (revision 28794) +++ /branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageCovariance.c (revision 28795) @@ -34,5 +34,6 @@ static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix const psKernel *kernel, // Convolution kernel - int x, int y // Coordinates in output covariance matrix + int x, int y, // Coordinates in output covariance matrix + float scale // Scale to apply ) { @@ -83,5 +84,5 @@ } - return sum; + return scale * sum; } @@ -91,5 +92,5 @@ PS_ASSERT_THREAD_JOB_NON_NULL(job, false); psAssert(job->args, "No job arguments"); - psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n); + psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n); psKernel *out = job->args->data[0]; // Output covariance matrix @@ -98,6 +99,7 @@ int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix - - out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y); + float scale = PS_SCALAR_VALUE(job->args->data[5], F32); // Scaling to apply + + out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale); return true; @@ -127,4 +129,5 @@ // Check for non-finite elements + double sumKernel = 0.0, sumKernel2 = 0.0; // Sum of the kernel for (int y = kernel->yMin; y <= kernel->yMax; y++) { for (int x = kernel->xMin; x <= kernel->xMax; x++) { @@ -135,4 +138,6 @@ return NULL; } + sumKernel += kernel->kernel[y][x]; + sumKernel2 += PS_SQR(kernel->kernel[y][x]); } } @@ -155,4 +160,5 @@ int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax; psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output + float scale = 1.0 / sumKernel2; // Scaling to apply for (int y = yMin; y <= yMax; y++) { @@ -165,4 +171,5 @@ PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32); PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, scale, PS_TYPE_F32); if (!psThreadJobAddPending(job)) { psFree(covar); @@ -170,9 +177,8 @@ } } else { - out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y); - } - } - } - psFree(covar); + out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale); + } + } + } if (threaded && !psThreadPoolWait(true)) { @@ -180,4 +186,6 @@ return false; } + + psFree(covar); return out; @@ -194,4 +202,5 @@ // Check for non-finite elements + double sumKernel2 = 0.0; // Sum of the squared kernel for (int y = kernel->yMin; y <= kernel->yMax; y++) { for (int x = kernel->xMin; x <= kernel->xMax; x++) { @@ -202,4 +211,5 @@ return NAN; } + sumKernel2 += PS_SQR(kernel->kernel[y][x]); } } @@ -215,5 +225,6 @@ } - float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor + float scale = 1.0 / sumKernel2; // Scale to apply + float factor = imageCovarianceCalculate(covar, kernel, 0, 0, scale); // Covariance factor psFree(covar); return factor; @@ -338,10 +349,9 @@ } } - psFree(covar); - if (threaded && !psThreadPoolWait(true)) { psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); return false; } + psFree(covar); return out; @@ -616,5 +626,5 @@ if (set && !threaded) { { - psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5); + psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 6); task->function = &imageCovarianceCalculateThread; psThreadTaskAdd(task); Index: /branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageInterpolate.c (revision 28794) +++ /branches/eam_branches/ipp-20100621/psLib/src/imageops/psImageInterpolate.c (revision 28795) @@ -427,32 +427,42 @@ // Determine the result of the interpolation after all the math has been done -#define INTERPOLATE_RESULT() \ - psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; /* Status of interpolation */ \ - *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage; \ - if (wantVariance) { \ - *varianceValue = sumVariance / (sumKernel2 - sumBad); \ - } \ - if (sumKernel == 0.0) { \ - /* No kernel contributions */ \ - if (haveMask && maskValue) { \ - *maskValue |= interp->badMask; \ - } \ - status = PS_INTERPOLATE_STATUS_BAD; \ - } else if (sumBad == 0) { \ - /* Completely good pixel */ \ - status = PS_INTERPOLATE_STATUS_GOOD; \ - } else if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { \ - /* Some pixels masked: poor pixel */ \ - if (haveMask && maskValue) { \ - *maskValue |= interp->poorMask; \ - } \ - status = PS_INTERPOLATE_STATUS_POOR; \ - } else { \ - /* Many pixels (or a few important pixels) masked: bad pixel */ \ - if (haveMask && maskValue) { \ - *maskValue |= interp->badMask; \ - } \ - status = PS_INTERPOLATE_STATUS_BAD; \ - } +static psImageInterpolateStatus interpolateResult(const psImageInterpolation *interp, + double *imageValue, double *varianceValue, + psImageMaskType *maskValue, + double sumImage, double sumVariance, double sumBad, + double sumKernel, double sumKernel2, + bool wantVariance, bool haveMask) +{ + *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage; + if (wantVariance) { + if (sumBad > 0) { + sumVariance *= sumKernel2 / (sumKernel2 - sumBad); + } + *varianceValue = sumVariance / PS_SQR(sumKernel); + } + if (sumKernel == 0.0) { + // No kernel contributions at all + if (haveMask && maskValue) { + *maskValue |= interp->badMask; + } + return PS_INTERPOLATE_STATUS_BAD; + } + if (sumBad == 0) { + // Completely good pixel + return PS_INTERPOLATE_STATUS_GOOD; + } + if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { + // Some pixels masked: poor pixel + if (haveMask && maskValue) { + *maskValue |= interp->poorMask; + } + return PS_INTERPOLATE_STATUS_POOR; + } + // Many pixels (or a few important pixels) masked: bad pixel + if (haveMask && maskValue) { + *maskValue |= interp->badMask; + } + return PS_INTERPOLATE_STATUS_BAD; +} // Interpolation engine for separable interpolation kernels @@ -703,6 +713,4 @@ } - INTERPOLATE_RESULT(); - psFree(xKernelNew); psFree(yKernelNew); @@ -710,5 +718,6 @@ psFree(yKernel2New); - return status; + return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad, + sumKernel, sumKernel2, wantVariance, haveMask); } @@ -861,7 +870,6 @@ } - INTERPOLATE_RESULT(); - - return status; + return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad, + sumKernel, sumKernel2, wantVariance, haveMask); } Index: /branches/eam_branches/ipp-20100621/psLib/src/math/psMinimizePowell.c =================================================================== --- /branches/eam_branches/ipp-20100621/psLib/src/math/psMinimizePowell.c (revision 28794) +++ /branches/eam_branches/ipp-20100621/psLib/src/math/psMinimizePowell.c (revision 28795) @@ -83,12 +83,12 @@ f(param + b * line) < f(param + c * line) a < b < c - + Algorithm: - + XXX completely ad hoc: start with the user-supplied starting parameter and call that b. Calculate a/c as a fractional amount smaller/larger than b. Repeat this process until a local minimum is found. - + XXX: new algorithm: start at x=0, expand in one direction until the function @@ -96,5 +96,5 @@ increases, or x is too large. If thst does not work, expand in the other direction. - + XXX: output bracket vector should be an input as well. *****************************************************************************/ @@ -335,5 +335,5 @@ along that vector and returns the offset along that vector at which the minimum is determined. - + XXX: This routine is not very efficient in terms of total evaluations of the function. @@ -491,5 +491,5 @@ points at which the function is varied are in the argument "coords" which is a psArray of psVectors: each vector represents a different coordinate. - + XXX: We do not use Brent's method. *****************************************************************************/ @@ -709,5 +709,5 @@ This routine is to be used with the psMinimizeChi2Powell() function below. and the psMinimizePowell() function above. - + The basic idea is calculate chi-squared for a set of params/coords/errors. This functions uses global variables to receive the function pointer, the @@ -764,5 +764,5 @@ This routine must minimize the chi-squared match of a set of data points and values for a possibly multi-dimensional function. - + The basic idea is to use the psMinimizePowell() function defined above. In order to do so, we defined above a function myPowellChi2Func() which takes @@ -780,4 +780,11 @@ psMinimizeChi2PowellFunc model) { + PS_ASSERT_VECTOR_NON_NULL(params, false); + PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); + PS_ASSERT_ARRAY_NON_NULL(coords, false); + PS_ASSERT_VECTOR_NON_NULL(value, false); + PS_ASSERT_VECTOR_TYPE(value, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(coords, value, false); + // Generate extended version of coords array, so we can pass in extra data to the chi^2 function psArray *newCoords = psArrayAlloc(coords->n + 3); @@ -791,5 +798,6 @@ newCoords->data[coords->n + 2] = &model; - bool success = psMinimizePowell(min, params, constraint->paramMask, newCoords, myPowellChi2Func); + bool success = psMinimizePowell(min, params, constraint ? constraint->paramMask : NULL, + newCoords, myPowellChi2Func); newCoords->data[coords->n - 1] = NULL; // We can't free the array with a function pointer on it