Index: /branches/eam_branches/ipp-20130509/psModules/src/imcombine/pmSubtractionSimple.c =================================================================== --- /branches/eam_branches/ipp-20130509/psModules/src/imcombine/pmSubtractionSimple.c (revision 35750) +++ /branches/eam_branches/ipp-20130509/psModules/src/imcombine/pmSubtractionSimple.c (revision 35750) @@ -0,0 +1,316 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmHDUUtils.h" +#include "pmSubtractionTypes.h" +#include "pmSubtraction.h" +#include "pmSubtractionParams.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionStamps.h" +#include "pmSubtractionEquation.h" +#include "pmSubtractionAnalysis.h" +#include "pmSubtractionMask.h" +#include "pmSubtractionThreads.h" +#include "pmSubtractionVisual.h" +#include "pmErrorCodes.h" + +#include "pmSubtractionSimple.h" + +bool pmSubtractionSimpleMatch(pmReadout *conv1, + pmReadout *conv2, + const pmReadout *ro1, + const pmReadout *ro2, + const psArray *sources, + int size, + psImageMaskType maskVal, + psImageMaskType maskBad, + psImageMaskType maskPoor + ) { + // + // We've already validated the input values at this level + float sig2fwhm = 2.0 * sqrt(2.0 * log(2.0)); + float fwhm1,fwhm2; + float sigma1,sigma2,sigmaKern; + int convolution_direction = 0; + psImage *image1; + psImage *mask1; + psImage *var1; + + psImage *image2; + psImage *mask2; + psImage *var2; + + psImage *imageC1; + psImage *maskC1; + psImage *varC1; + + psImage *imageC2; + psImage *maskC2; + psImage *varC2; + + // Allocate images, as this is usually done by subtractionMatchAlloc after this function is called. + int numCols = ro1->image->numCols; + int numRows = ro1->image->numRows; + + if (conv1) { + if (!conv1->image) { + conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv1->image, NAN); + if (ro1->variance) { + if (!conv1->variance) { + conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv1->variance, NAN); + } + if (!conv1->mask) { + conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } + psImageInit(conv1->mask, maskBad); + } + if (conv2) { + if (!conv2->image) { + conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv2->image, NAN); + if (ro2->variance) { + if (!conv2->variance) { + conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv2->variance, NAN); + } + if (!conv2->mask) { + conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } + psImageInit(conv2->mask, maskBad); + } + + // Assign local aliases to images + image1 = ro1->image; + mask1 = ro1->mask; + var1 = ro1->variance; + + image2 = ro2->image; + mask2 = ro2->mask; + var2 = ro2->variance; + + if (conv1) { + imageC1 = conv1->image; + maskC1 = conv1->mask; + varC1 = conv1->variance; + } + + if (conv2) { + imageC2 = conv2->image; + maskC2 = conv2->mask; + varC2 = conv2->variance; + } + + // + // Determine Gaussian widths + pmSubtractionGetFWHMs(&fwhm1,&fwhm2); + sigma1 = fwhm1 / sig2fwhm; + sigma2 = fwhm2 / sig2fwhm; + if (sigma1 > sigma2) { + convolution_direction = 2; + sigmaKern = sqrt(PS_SQR(sigma1) - PS_SQR(sigma2)); + } + else if (sigma1 < sigma2) { + convolution_direction = 1; + sigmaKern = sqrt(PS_SQR(sigma2) - PS_SQR(sigma1)); + } + if (!conv1) { + convolution_direction = 2; + } + if (!conv2) { + convolution_direction = 1; + } + + // + // Determine Normalization scaling + psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity. + + + // + // Do convolutions + psImage *temp = psImageAlloc(numCols,numRows,PS_TYPE_F32); + psImage *tempMask = psImageAlloc(numCols,numRows,PS_TYPE_IMAGE_MASK); + psImage *tempVar = psImageAlloc(numCols,numRows,PS_TYPE_F32); + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + temp->data.F32[y][x] = 0.0; + tempVar->data.F32[y][x] = 0.0; + tempMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + + // Copy the image we're not convolving into the convolved output + if (convolution_direction == 1) { + if (conv1) { + imageC1->data.F32[y][x] = 0.0; + maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + varC1->data.F32[y][x] = 0.0; + } + if (conv2) { + imageC2->data.F32[y][x] = image2->data.F32[y][x]; + maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask2->data.PS_TYPE_IMAGE_MASK_DATA[y][x]; + varC2->data.F32[y][x] = var2->data.F32[y][x]; + } + } + else if (convolution_direction == 2) { + if (conv2) { + imageC2->data.F32[y][x] = 0.0; + maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + varC2->data.F32[y][x] = 0.0; + } + if (conv1) { + imageC1->data.F32[y][x] = image1->data.F32[y][x]; + maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask1->data.PS_TYPE_IMAGE_MASK_DATA[y][x]; + varC1->data.F32[y][x] = var1->data.F32[y][x]; + } + } + + // Do y-direction convolution + for (int v = -size; v <= size; v++) { + if ((y-v >= 0)&&(y-v < numRows)) { + if (convolution_direction == 1) { + // Insert mask/finite check here. + temp->data.F32[y][x] += kernelVec->data.F32[v + size] * image1->data.F32[y - v][x]; + tempVar->data.F32[y][x] += kernelVec->data.F32[v + size] * var1->data.F32[y - v][x]; + } + else if (convolution_direction == 2) { + // And here. + temp->data.F32[y][x] += kernelVec->data.F32[v + size] * image2->data.F32[y - v][x]; + tempVar->data.F32[y][x] += kernelVec->data.F32[v + size] * var2->data.F32[y - v][x]; + } + } + } + } + } + + // Do x-direction convolution + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + for (int u = -size; u <= size; u++) { + if ((x-u >= 0)&&(x-u < numCols)) { + if (convolution_direction == 1) { + // And here + imageC1->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u]; + varC1->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u]; + } + else if (convolution_direction == 2) { + // And here + imageC2->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u]; + varC2->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u]; + } + } + } + } + } + psFree(temp); + psFree(kernelVec); + + // + // Do normalization + float normalization = 1.0; + + // Something with the source list here + + // Apply normalization + if (normalization != 1.0) { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if ((conv1)&&((convolution_direction == 1))) { + imageC1->data.F32[y][x] *= normalization; + } + else if ((conv2)&&(convolution_direction == 2)) { + imageC2->data.F32[y][x] *= normalization; + } + } + } + } + + + // + // Make a fake pmSubtractionKernels element so we can add it appropriately. + // I call it fake because we've successfully done everything at this point + // without having to define these things. + psVector *fwhms = psVectorAlloc(1,PS_TYPE_F32); + fwhms->data.F32[0] = sigmaKern * sig2fwhm; + psVector *orders = psVectorAlloc(1,PS_TYPE_S32); + orders->data.S32[0] = 0; + psRegion bounds; + bounds.x0 = 0; + bounds.y0 = 0; + bounds.x1 = numCols; + bounds.y1 = numRows; + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(1,PM_SUBTRACTION_KERNEL_SIMPLE, + size,fwhms,orders,0,0.0,bounds, + convolution_direction); + pmSubtractionKernelsMakeDescription(kernels); // Need this defined + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_SIMPLE,0,0,size,sigmaKern); + kernels->widths->data.F32[0] = sigmaKern * sig2fwhm; + kernels->u->data.S32[0] = 0; + kernels->v->data.S32[0] = 0; + if (kernels->preCalc->data[0]) { + psFree(kernels->preCalc->data[0]); + } + kernels->preCalc->data[0] = preCalc; + kernels->solution1 = psVectorAlloc(3,PS_TYPE_F64); + kernels->solution1->data.F32[0] = 1.0; + kernels->solution1->data.F32[1] = 0.0; + kernels->solution1->data.F32[2] = 0.0; + kernels->solution1err = psVectorAlloc(3,PS_TYPE_F32); + kernels->solution1err->data.F32[0] = 0.0; + kernels->solution1err->data.F32[1] = 0.0; + kernels->solution1err->data.F32[2] = 0.0; + + // + // Actually add it to the headers + psMetadata *outAnalysis = psMetadataAlloc(); + psMetadata *outHeader = psMetadataAlloc(); + if (!pmSubtractionAnalysis(outAnalysis,outHeader,kernels,NULL,numCols,numRows)) { + psError(psErrorCodeLast(),false,"Unable to generate QA data"); + psFree(fwhms); + psFree(orders); + psFree(preCalc); + psFree(kernels); + } + psFree(fwhms); + psFree(orders); + psFree(preCalc); + psFree(kernels); + + if (conv1) { + conv1->analysis = psMetadataCopy(conv1->analysis, outAnalysis); + } + if (conv2) { + conv2->analysis = psMetadataCopy(conv2->analysis, outAnalysis); + } + + if (conv1 && conv1->parent) { + pmHDU *hdu = pmHDUFromCell(conv1->parent); + if (hdu) { + hdu->header = psMetadataCopy(hdu->header, outHeader); + } + } + if (conv2 && conv2->parent) { + pmHDU *hdu = pmHDUFromCell(conv2->parent); + if (hdu) { + hdu->header = psMetadataCopy(hdu->header, outHeader); + } + } + psFree(outAnalysis); + psFree(outHeader); + + // + // Return + return(true); +} Index: /branches/eam_branches/ipp-20130509/psModules/src/imcombine/pmSubtractionSimple.h =================================================================== --- /branches/eam_branches/ipp-20130509/psModules/src/imcombine/pmSubtractionSimple.h (revision 35750) +++ /branches/eam_branches/ipp-20130509/psModules/src/imcombine/pmSubtractionSimple.h (revision 35750) @@ -0,0 +1,19 @@ +#ifndef PM_SUBTRACTION_SIMPLE_H +#define PM_SUBTRACTION_SIMPLE_H + +#include + +#include +#include +#include + +bool pmSubtractionSimpleMatch(pmReadout *conv1, + pmReadout *conv2, + const pmReadout *ro1, + const pmReadout *ro2, + const psArray *sources, + int size + ); + + +#endif Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO.c (revision 35749) +++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO.c (revision 35750) @@ -58,6 +58,6 @@ #define BLANK_HEADERS "BLANK.HEADERS" // Name of metadata in camera configuration containing header names // for putting values into a blank PHU -static bool pmReadoutReadXSRC(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex); -static bool pmReadoutReadXFIT(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); +static bool pmReadoutReadXSRC(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex); +static bool pmReadoutReadXFIT(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); @@ -1093,5 +1093,5 @@ if (XSRC_OUTPUT && xsrcname) { // a cmf file may have an XSRC extension, but it is not required - if (!pmReadoutReadXSRC(file, exttype, hdu->header, xsrcname, sources, sourceIndex)) { + if (!pmReadoutReadXSRC(file, readout, exttype, hdu->header, xsrcname, sources, sourceIndex)) { // do anything? } @@ -1100,5 +1100,5 @@ if (XFIT_OUTPUT && xfitname) { // a cmf file may have an XFIT extension, but it is not required - if (!pmReadoutReadXFIT(file, exttype, hdu->header, xfitname, sources, sourceIndex)) { + if (!pmReadoutReadXFIT(file, readout, exttype, hdu->header, xfitname, sources, sourceIndex)) { // do anything? } @@ -1243,5 +1243,5 @@ // XXX: We might be able to macroize this and reuse for the other types -static bool pmReadoutReadXSRC(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex) +static bool pmReadoutReadXSRC(pmFPAfile *file, pmReadout *readout, char *exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex) { if (!psFitsMoveExtNameClean (file->fits, xsrcname)) { @@ -1263,5 +1263,5 @@ # define PM_SOURCES_READ_XSRC(NAME,TYPE) \ if (!strcmp (exttype, NAME)) { \ - status = pmSourcesRead_##TYPE##_XSRC(file->fits, hduHeader, sources, sourceIndex); \ + status = pmSourcesRead_##TYPE##_XSRC(file->fits, readout, hduHeader, tableHeader, sources, sourceIndex); \ } @@ -1282,5 +1282,5 @@ } -static bool pmReadoutReadXFIT(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) +static bool pmReadoutReadXFIT(pmFPAfile *file, pmReadout *readout, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) { if (!psFitsMoveExtNameClean (file->fits, extname)) { @@ -1302,5 +1302,5 @@ # define PM_SOURCES_READ_XFIT(NAME,TYPE) \ if (!strcmp (exttype, NAME)) { \ - status = pmSourcesRead_##TYPE##_XFIT(file->fits, hduHeader, sources, sourceIndex); \ + status = pmSourcesRead_##TYPE##_XFIT(file->fits, readout, hduHeader, tableHeader, sources, sourceIndex); \ } @@ -1340,5 +1340,5 @@ # define PM_SOURCES_READ_XRAD(NAME,TYPE) \ if (!strcmp (exttype, NAME)) { \ - status = pmSourcesRead_##TYPE##_XRAD(file->fits, readout, hduHeader, sources, sourceIndex); \ + status = pmSourcesRead_##TYPE##_XRAD(file->fits, readout, hduHeader, tableHeader, sources, sourceIndex); \ } Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO.h =================================================================== --- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO.h (revision 35749) +++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO.h (revision 35750) @@ -22,7 +22,7 @@ bool pmSourcesWrite_##TYPE##_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); \ psArray *pmSourcesRead_##TYPE (psFits *fits, psMetadata *header); \ - bool pmSourcesRead_##TYPE##_XSRC (psFits *fits, psMetadata *header, psArray *sources, long *index); \ - bool pmSourcesRead_##TYPE##_XFIT (psFits *fits, psMetadata *header, psArray *sources, long *index); \ - bool pmSourcesRead_##TYPE##_XRAD (psFits *fits, pmReadout *readout, psMetadata *header, psArray *sources, long *index);\ + bool pmSourcesRead_##TYPE##_XSRC (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index); \ + bool pmSourcesRead_##TYPE##_XFIT (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index); \ + bool pmSourcesRead_##TYPE##_XRAD (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index);\ // All of these functions need to use the same API, even if not all elements are used in a specific case Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO_CMF.c.in (revision 35749) +++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceIO_CMF.c.in (revision 35750) @@ -457,6 +457,14 @@ bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + // First look for radial bin definition in readout->analysis ... + psVector *radMin = psMetadataLookupPtr (&status, readout->analysis, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, readout->analysis, "RADIAL.ANNULAR.BINS.UPPER"); + if (!radMin) { + // .. if not found use the recipe values + radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + } + psAssert(radMin != NULL, "unable to find RADIAL.ANNULAR.BINS.LOWER"); + psAssert(radMax != NULL, "unable to find RADIAL.ANNULAR.BINS.LOWER"); psAssert (radMin->n == radMax->n, "inconsistent annular bins"); @@ -610,5 +618,66 @@ } -bool pmSourcesRead_CMF_@CMFMODE@_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +static bool setRadialBinsInAnalysis(pmReadout *readout, psMetadata *tableHeader) +{ + if (!readout->analysis) { + readout->analysis = psMetadataAlloc(); + } + + bool status; + psVector *oldRMin = psMetadataLookupVector(&status, readout->analysis, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *oldRMax = psMetadataLookupVector(&status, readout->analysis, "RADIAL.ANNULAR.BINS.UPPER"); + + psVector *rMin = psVectorAllocEmpty(20, PS_TYPE_F32); + psVector *rMax = psVectorAllocEmpty(20, PS_TYPE_F32); + + for (int i = 0; ; i++) { + char key [24]; + sprintf(key, "RMIN_%02d", i); + psF32 rMinVal = psMetadataLookupF32(&status, tableHeader, key); + if (!status) { + break; + } + psVectorAppend(rMin, rMinVal); + + sprintf(key, "RMAX_%02d", i); + psF32 rMaxVal = psMetadataLookupF32(&status, tableHeader, key); + if (!status) { + break; + } + psVectorAppend(rMax, rMaxVal); + } + + if (rMin->n != rMax->n) { + psError(PS_ERR_UNKNOWN, true, "number of RMIN entries %ld does not equal number of RMAX entries %ld", + rMin->n, rMax->n); + return false; + } + + if (oldRMin) { + if (oldRMin->n != rMin->n) { + psError(PS_ERR_UNKNOWN, true, "number of RMIN entries in header: %ld does not equal number of RMAX entries in analysis: %ld", + rMin->n, oldRMin->n); + return false; + } + } else { + psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.LOWER", PS_META_REPLACE, "", rMin); + } + psFree(rMin); + + if (oldRMax) { + if (oldRMax->n != rMax->n) { + psError(PS_ERR_UNKNOWN, true, "number of RMIN entries in header: %ld does not equal number of RMAX entries in analysis: %ld", + rMax->n, oldRMax->n); + return false; + } + } else { + psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.UPPER", PS_META_REPLACE, "", rMax); + } + psFree(rMax); + + return true; +} + +bool pmSourcesRead_CMF_@CMFMODE@_XSRC(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psMetadata *tableHeader, psArray *sources, long *sourceIndex) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -626,4 +695,9 @@ float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME"); float magOffset = zeropt + 2.5*log10(exptime); + + if (!setRadialBinsInAnalysis(readout, tableHeader)) { + psError(PS_ERR_UNKNOWN, false, "Failed to save radial bins in analysis"); + return false; + } for (long i = 0; i < numSources; i++) { @@ -920,5 +994,5 @@ } -bool pmSourcesRead_CMF_@CMFMODE@_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +bool pmSourcesRead_CMF_@CMFMODE@_XFIT(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psMetadata *tableHeader, psArray *sources, long *sourceIndex) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -1043,8 +1117,14 @@ // we use this just to define the output vectors (which must be present for all objects) - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + // First look for radial bin definition in readout->analysis ... + psVector *radMin = psMetadataLookupPtr (&status, readout->analysis, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, readout->analysis, "RADIAL.ANNULAR.BINS.UPPER"); + if (!radMin) { + // .. if not found use the recipe values + radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + } + psAssert (radMin, "this must have been defined and tested earlier!"); psAssert (radMax, "this must have been defined and tested earlier!"); - psAssert (radMax->n, "this must have been defined and tested earlier!"); psAssert (radMin->n == radMax->n, "inconsistent annular bins"); @@ -1167,5 +1247,5 @@ } -bool pmSourcesRead_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +bool pmSourcesRead_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psMetadata *tableHeader, psArray *sources, long *sourceIndex) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -1176,4 +1256,9 @@ if (numSources == 0) { psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n"); + return false; + } + + if (!setRadialBinsInAnalysis(readout, tableHeader)) { + psError(PS_ERR_UNKNOWN, false, "Failed to save radial bins in analysis"); return false; }