Index: /trunk/psModules/src/detrend/pmBias.c =================================================================== --- /trunk/psModules/src/detrend/pmBias.c (revision 18959) +++ /trunk/psModules/src/detrend/pmBias.c (revision 18960) @@ -19,5 +19,7 @@ #include "pmBias.h" -bool pmBiasSubtractScan_Threaded (psThreadJob *job) { +bool pmBiasSubtractScan_Threaded(const psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); pmReadout *in = job->args->data[0]; pmReadout *sub = job->args->data[1]; @@ -63,6 +65,6 @@ // image. The bias image is subtracted in place from the input image. bool pmBiasSubtractFrame(pmReadout *in, // Input readout - pmReadout *sub, // Readout to be subtracted from input - float scale // Scale to apply before subtracting + pmReadout *sub, // Readout to be subtracted from input + float scale // Scale to apply before subtracting ) { @@ -121,7 +123,7 @@ bool threaded = true; int scanRows = pmDetrendGetScanRows(); - if (scanRows == 0) { - threaded = false; - scanRows = inImage->numRows; + if (scanRows == 0) { + threaded = false; + scanRows = inImage->numRows; } @@ -129,43 +131,40 @@ int rowStop = PS_MIN (rowStart + scanRows, inImage->numRows); -# define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \ - psScalar *scalar = psScalarAlloc(VALUE, TYPE); \ - psArrayAdd(ARRAY, 1, scalar); \ - psFree (scalar); } - if (threaded) { - // allocate a job, construct the arguments for this job - psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_BIAS"); - psArrayAdd (job->args, 1, in); - psArrayAdd (job->args, 1, sub); - PS_ARRAY_ADD_SCALAR (job->args, scale, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); - - // ppImageDetrendReadout(config, options, view) - if (!psThreadJobAddPending (job)) { - return false; - } + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_BIAS"); + psArrayAdd (job->args, 1, in); + psArrayAdd (job->args, 1, sub); + PS_ARRAY_ADD_SCALAR (job->args, scale, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); + + // ppImageDetrendReadout(config, options, view) + if (!psThreadJobAddPending (job)) { + psFree(job); + return false; + } + psFree(job); } else { - pmBiasSubtractScan (in, sub, scale, xOffset, yOffset, rowStart, rowStop); + pmBiasSubtractScan (in, sub, scale, xOffset, yOffset, rowStart, rowStop); } } if (threaded) { - // wait here for the threaded jobs to finish - if (!psThreadPoolWait ()) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - return false; - } - fprintf (stderr, "success for threaded jobs\n"); - - // each job records its own goodPixel values; sum them here - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - psFree (job); - } + // wait here for the threaded jobs to finish + if (!psThreadPoolWait(false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + fprintf (stderr, "success for threaded jobs\n"); + + // each job records its own goodPixel values; sum them here + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + psFree (job); + } } Index: /trunk/psModules/src/detrend/pmBias.h =================================================================== --- /trunk/psModules/src/detrend/pmBias.h (revision 18959) +++ /trunk/psModules/src/detrend/pmBias.h (revision 18960) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-05 01:24:47 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-08 18:09:07 $ * Copyright 2004--2006 Institute for Astronomy, University of Hawaii */ @@ -30,9 +30,9 @@ // image. The bias image is subtracted in place from the input image. bool pmBiasSubtractFrame(pmReadout *in, // Input readout - pmReadout *sub, // Readout to be subtracted from input - float scale // Scale to apply before subtracting + pmReadout *sub, // Readout to be subtracted from input + float scale // Scale to apply before subtracting ); -bool pmBiasSubtractScan_Threaded (psThreadJob *job); +bool pmBiasSubtractScan_Threaded(const psThreadJob *job); bool pmBiasSubtractScan (pmReadout *in, pmReadout *sub, float scan, int xOffset, int yOffset, int rowStart, int rowStop); Index: /trunk/psModules/src/detrend/pmDark.c =================================================================== --- /trunk/psModules/src/detrend/pmDark.c (revision 18959) +++ /trunk/psModules/src/detrend/pmDark.c (revision 18960) @@ -102,5 +102,5 @@ bool pmDarkCombinePrepare(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept) { - psArray *values = psArrayAlloc(inputs->n); + psArray *values = psArrayAlloc(inputs->n); psVector *roMask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for bad readouts psVector *norm = normConcept ? psVectorAlloc(inputs->n, PS_TYPE_F32) : NULL; // Normalizations for each @@ -117,21 +117,21 @@ if (!norm) continue; - pmReadout *readout = inputs->data[i]; // Readout of interest - float normValue; // Normalisation value - if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) { - psWarning("Unable to find value of %s for readout %d", normConcept, i); - roMask->data.U8[i] = 0xff; - norm->data.F32[i] = NAN; - numBadInputs++; - continue; - } - if (normValue == 0.0) { - psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i); - roMask->data.U8[i] = 0xff; - norm->data.F32[i] = NAN; - numBadInputs++; - continue; - } - norm->data.F32[i] = 1.0 / normValue; + pmReadout *readout = inputs->data[i]; // Readout of interest + float normValue; // Normalisation value + if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) { + psWarning("Unable to find value of %s for readout %d", normConcept, i); + roMask->data.U8[i] = 0xff; + norm->data.F32[i] = NAN; + numBadInputs++; + continue; + } + if (normValue == 0.0) { + psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i); + roMask->data.U8[i] = 0xff; + norm->data.F32[i] = NAN; + numBadInputs++; + continue; + } + norm->data.F32[i] = 1.0 / normValue; } @@ -144,5 +144,5 @@ psFree(roMask); psFree(orders); - psFree(norm); + psFree(norm); return false; } @@ -215,6 +215,6 @@ } - pmReadoutStackDefineOutput(readout, col0, row0, numCols, numRows, false, false, 0); - psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", col0, row0); + pmReadoutStackDefineOutput(readout, col0, row0, numCols, numRows, false, false, 0); + psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", col0, row0); } @@ -269,11 +269,11 @@ bool mdok = false; - // retrieve the required parameter vectors - psArray *values = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES"); psAssert (values, "values not supplied"); - psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK"); psAssert (roMask, "roMask not supplied"); + // retrieve the required parameter vectors + psArray *values = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES"); psAssert (values, "values not supplied"); + psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK"); psAssert (roMask, "roMask not supplied"); psPolynomialMD *poly = psMetadataLookupPtr(&mdok, output->analysis, "DARK.POLY"); psAssert (poly, "orders not supplied"); // retrieve the norm vector, if supplied - psVector *norm = psMetadataLookupPtr(&mdok, output->analysis, "DARK.NORM"); + psVector *norm = psMetadataLookupPtr(&mdok, output->analysis, "DARK.NORM"); // retrieve the 'counts' and 'sigma' images @@ -367,5 +367,8 @@ } -bool pmDarkApplyScan_Threaded (psThreadJob *job) { +bool pmDarkApplyScan_Threaded(const psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + pmReadout *readout = job->args->data[0]; pmCell *dark = job->args->data[1]; @@ -478,7 +481,7 @@ bool threaded = true; int scanRows = pmDetrendGetScanRows(); - if (scanRows == 0) { - threaded = false; - scanRows = readout->image->numRows; + if (scanRows == 0) { + threaded = false; + scanRows = readout->image->numRows; } @@ -486,44 +489,41 @@ int rowStop = PS_MIN (rowStart + scanRows, readout->image->numRows); -# define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \ - psScalar *scalar = psScalarAlloc(VALUE, TYPE); \ - psArrayAdd(ARRAY, 1, scalar); \ - psFree (scalar); } - if (threaded) { - // allocate a job, construct the arguments for this job - psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_DARK"); - psArrayAdd (job->args, 1, readout); - psArrayAdd (job->args, 1, dark); - psArrayAdd (job->args, 1, poly); - psArrayAdd (job->args, 1, values); - PS_ARRAY_ADD_SCALAR (job->args, bad, PS_TYPE_MASK); - PS_ARRAY_ADD_SCALAR (job->args, doNorm, PS_TYPE_U8); - PS_ARRAY_ADD_SCALAR (job->args, norm, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); - - // ppImageDetrendReadout(config, options, view) - if (!psThreadJobAddPending (job)) { - return false; - } + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_DARK"); + psArrayAdd (job->args, 1, readout); + psArrayAdd (job->args, 1, dark); + psArrayAdd (job->args, 1, poly); + psArrayAdd (job->args, 1, values); + PS_ARRAY_ADD_SCALAR (job->args, bad, PS_TYPE_MASK); + PS_ARRAY_ADD_SCALAR (job->args, doNorm, PS_TYPE_U8); + PS_ARRAY_ADD_SCALAR (job->args, norm, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); + + // ppImageDetrendReadout(config, options, view) + if (!psThreadJobAddPending (job)) { + psFree(job); + return false; + } + psFree(job); } else { - pmDarkApplyScan (readout, dark, poly, values, bad, doNorm, norm, rowStart, rowStop); + pmDarkApplyScan (readout, dark, poly, values, bad, doNorm, norm, rowStart, rowStop); } } if (threaded) { - // wait here for the threaded jobs to finish - if (!psThreadPoolWait ()) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - return false; - } - fprintf (stderr, "success for threaded jobs\n"); - - // free the done jobs - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - psFree (job); - } + // wait here for the threaded jobs to finish + if (!psThreadPoolWait(false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + fprintf (stderr, "success for threaded jobs\n"); + + // free the done jobs + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + psFree (job); + } } Index: /trunk/psModules/src/detrend/pmDark.h =================================================================== --- /trunk/psModules/src/detrend/pmDark.h (revision 18959) +++ /trunk/psModules/src/detrend/pmDark.h (revision 18960) @@ -27,7 +27,7 @@ // Combine darks -- preparation step bool pmDarkCombinePrepare(pmCell *output, // Output cell; readouts will be attached - const psArray *inputs, // Input readouts for combination - psArray *ordinates, // Ordinates for fitting - const char *normConcept // Concept name to use to divide input pixel values + const psArray *inputs, // Input readouts for combination + psArray *ordinates, // Ordinates for fitting + const char *normConcept // Concept name to use to divide input pixel values ); @@ -40,5 +40,5 @@ ); -bool pmDarkApplyScan_Threaded (psThreadJob *job); +bool pmDarkApplyScan_Threaded(const psThreadJob *job); bool pmDarkApplyScan (pmReadout *readout, const pmCell *dark, psPolynomialMD *poly, psVector *values, psMaskType bad, bool doNorm, float norm, int rowStart, int rowStop); Index: /trunk/psModules/src/detrend/pmDetrendThreads.c =================================================================== --- /trunk/psModules/src/detrend/pmDetrendThreads.c (revision 18959) +++ /trunk/psModules/src/detrend/pmDetrendThreads.c (revision 18960) @@ -18,32 +18,45 @@ #include "pmDetrendThreads.h" -static int scanRows = 0; +static int scanRows = 0; // Number of rows to work on at once -int pmDetrendGetScanRows () { +int pmDetrendGetScanRows(void) +{ return scanRows; -} +} -bool pmDetrendSetThreadTasks (int newScanRows) { +bool pmDetrendSetThreadTasks (int newScanRows) +{ + psAssert(scanRows == 0, "programming error: program called pmDetrendSetThreadTasks twice"); - psAssert (scanRows == 0, "programming error: program called pmDetrendSetThreadTasks twice"); + PS_ASSERT_INT_POSITIVE(newScanRows, false); scanRows = newScanRows; - psThreadTask *task = NULL; + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_BIAS", 7); + task->function = &pmBiasSubtractScan_Threaded; + psThreadTaskAdd(task); + psFree(task); + } - task = psThreadTaskAlloc ("PSMODULES_DETREND_BIAS", 7); - task->function = &pmBiasSubtractScan_Threaded; - psThreadTaskAdd (task); + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_DARK", 9); + task->function = &pmDarkApplyScan_Threaded; + psThreadTaskAdd(task); + psFree(task); + } - task = psThreadTaskAlloc ("PSMODULES_DETREND_DARK", 9); - task->function = &pmDarkApplyScan_Threaded; - psThreadTaskAdd (task); + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_SHUTTER", 7); + task->function = &pmShutterCorrectionApplyScan_Threaded; + psThreadTaskAdd(task); + psFree(task); + } - task = psThreadTaskAlloc ("PSMODULES_DETREND_SHUTTER", 7); - task->function = &pmShutterCorrectionApplyScan_Threaded; - psThreadTaskAdd (task); - - task = psThreadTaskAlloc ("PSMODULES_DETREND_FLAT", 9); - task->function = &pmFlatFieldScan_Threaded; - psThreadTaskAdd (task); + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_FLAT", 9); + task->function = &pmFlatFieldScan_Threaded; + psThreadTaskAdd(task); + psFree(task); + } return true; Index: /trunk/psModules/src/detrend/pmFlatField.c =================================================================== --- /trunk/psModules/src/detrend/pmFlatField.c (revision 18959) +++ /trunk/psModules/src/detrend/pmFlatField.c (revision 18960) @@ -13,5 +13,8 @@ #include "pmDetrendThreads.h" -bool pmFlatFieldScan_Threaded (psThreadJob *job) { +bool pmFlatFieldScan_Threaded(const psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psImage *inImage = job->args->data[0]; psImage *inMask = job->args->data[1]; @@ -29,6 +32,6 @@ // Macro for all PS types -#define FLAT_DIVISION_CASE(TYPE, SPECIAL) \ - case PS_TYPE_##TYPE: \ +#define FLAT_DIVISION_CASE(TYPE, SPECIAL) \ + case PS_TYPE_##TYPE: \ for (int j = rowStart; j < rowStop; j++) { \ for (int i = 0; i < inImage->numCols; i++) { \ @@ -129,7 +132,7 @@ bool threaded = true; int scanRows = pmDetrendGetScanRows(); - if (scanRows == 0) { - threaded = false; - scanRows = inImage->numRows; + if (scanRows == 0) { + threaded = false; + scanRows = inImage->numRows; } @@ -137,44 +140,41 @@ int rowStop = PS_MIN (rowStart + scanRows, inImage->numRows); -# define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \ - psScalar *scalar = psScalarAlloc(VALUE, TYPE); \ - psArrayAdd(ARRAY, 1, scalar); \ - psFree (scalar); } - if (threaded) { - // allocate a job, construct the arguments for this job - psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_FLAT"); - psArrayAdd (job->args, 1, inImage); - psArrayAdd (job->args, 1, inMask); - psArrayAdd (job->args, 1, flatImage); - psArrayAdd (job->args, 1, flatMask); - PS_ARRAY_ADD_SCALAR (job->args, badFlat, PS_TYPE_U8); - PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_FLAT"); + psArrayAdd (job->args, 1, inImage); + psArrayAdd (job->args, 1, inMask); + psArrayAdd (job->args, 1, flatImage); + psArrayAdd (job->args, 1, flatMask); + PS_ARRAY_ADD_SCALAR (job->args, badFlat, PS_TYPE_U8); + PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); - // ppImageDetrendReadout(config, options, view) - if (!psThreadJobAddPending (job)) { - return false; - } + // ppImageDetrendReadout(config, options, view) + if (!psThreadJobAddPending (job)) { + psFree(job); + return false; + } + psFree(job); } else { - pmFlatFieldScan (inImage, inMask, flatImage, flatMask, badFlat, xOffset, yOffset, rowStart, rowStop); + pmFlatFieldScan (inImage, inMask, flatImage, flatMask, badFlat, xOffset, yOffset, rowStart, rowStop); } } if (threaded) { - // wait here for the threaded jobs to finish - if (!psThreadPoolWait ()) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - return false; - } - fprintf (stderr, "success for threaded jobs\n"); + // wait here for the threaded jobs to finish + if (!psThreadPoolWait(false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + fprintf (stderr, "success for threaded jobs\n"); - // free done jobs - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - psFree (job); - } + // free done jobs + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + psFree (job); + } } Index: /trunk/psModules/src/detrend/pmFlatField.h =================================================================== --- /trunk/psModules/src/detrend/pmFlatField.h (revision 18959) +++ /trunk/psModules/src/detrend/pmFlatField.h (revision 18960) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-05 01:24:47 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-08 18:09:07 $ * Copyright 2004-2006 Institute for Astronomy, University of Hawaii */ @@ -28,5 +28,5 @@ ); -bool pmFlatFieldScan_Threaded (psThreadJob *job); +bool pmFlatFieldScan_Threaded(const psThreadJob *job); bool pmFlatFieldScan (psImage *inImage, psImage *inMask, psImage *flatImage, psImage *flatMask, psMaskType badFlag, int xOffset, int yOffset, int rowStart, int rowStop); Index: /trunk/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 18959) +++ /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 18960) @@ -122,5 +122,5 @@ // Iterate only - for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++); + for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++); if (index == N - 1) { @@ -344,8 +344,8 @@ psVector *resid = psVectorAlloc (exptime->n, PS_TYPE_F32); for (int i = 0; i < exptime->n; i++) { - float fitCounts = corr->scale * (exptime->data.F32[i] + corr->offset) / (exptime->data.F32[i] + corr->offref); - resid->data.F32[i] = counts->data.F32[i] - fitCounts; - } - + float fitCounts = corr->scale * (exptime->data.F32[i] + corr->offset) / (exptime->data.F32[i] + corr->offref); + resid->data.F32[i] = counts->data.F32[i] - fitCounts; + } + psStats *rawStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); psStats *resStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); @@ -353,5 +353,5 @@ psVectorStats (resStats, resid, NULL, NULL, 0); - // XXX temporary hard-wired minimum stdev improvement factor + // XXX temporary hard-wired minimum stdev improvement factor psTrace("psModules.detrend", 3, "raw scatter %f vs res scatter %f\n", rawStats->sampleStdev, resStats->sampleStdev); if (rawStats->sampleStdev / resStats->sampleStdev < 1.5) corr->valid = false; @@ -433,5 +433,5 @@ numRows = image->numRows; numCols = image->numCols; - // define the reference region : a box of size 'size' at the center + // define the reference region : a box of size 'size' at the center refRegion = psRegionForSquare(0.5 * numCols, 0.5 * numRows, size); // Set up the sample regions : boxes of size 'size' at the 4 image corners @@ -656,5 +656,8 @@ -bool pmShutterCorrectionApplyScan_Threaded (psThreadJob *job) { +bool pmShutterCorrectionApplyScan_Threaded(const psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psImage *image = job->args->data[0]; psImage *shutterImage = job->args->data[1]; @@ -672,12 +675,12 @@ { for (int y = rowStart; y < rowStop; y++) { - for (int x = 0; x < image->numCols; x++) { - if (mask && !isfinite(shutterImage->data.F32[y][x])) { - mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; - image->data.F32[y][x] = NAN; - continue; - } - image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]); - } + for (int x = 0; x < image->numCols; x++) { + if (mask && !isfinite(shutterImage->data.F32[y][x])) { + mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + image->data.F32[y][x] = NAN; + continue; + } + image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]); + } } return true; @@ -732,7 +735,7 @@ bool threaded = true; int scanRows = pmDetrendGetScanRows(); - if (scanRows == 0) { - threaded = false; - scanRows = image->numRows; + if (scanRows == 0) { + threaded = false; + scanRows = image->numRows; } @@ -756,45 +759,42 @@ psFree (line); } else { - for (int rowStart = 0; rowStart < image->numRows; rowStart += scanRows) { - int rowStop = PS_MIN (rowStart + scanRows, image->numRows); - -# define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \ - psScalar *scalar = psScalarAlloc(VALUE, TYPE); \ - psArrayAdd(ARRAY, 1, scalar); \ - psFree (scalar); } - - if (threaded) { - // allocate a job, construct the arguments for this job - psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_SHUTTER"); - psArrayAdd (job->args, 1, image); - psArrayAdd (job->args, 1, shutterImage); - psArrayAdd (job->args, 1, mask); - PS_ARRAY_ADD_SCALAR (job->args, exptime, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR (job->args, blank, PS_TYPE_MASK); - PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); - - // ppImageDetrendReadout(config, options, view) - if (!psThreadJobAddPending (job)) { - return false; - } - } else { - pmShutterCorrectionApplyScan (image, shutterImage, mask, exptime, blank, rowStart, rowStop); - } - } - if (threaded) { - // wait here for the threaded jobs to finish - if (!psThreadPoolWait ()) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - return false; - } - fprintf (stderr, "success for threaded jobs\n"); - - // free the done jobs - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - psFree (job); - } - } + for (int rowStart = 0; rowStart < image->numRows; rowStart += scanRows) { + int rowStop = PS_MIN (rowStart + scanRows, image->numRows); + + if (threaded) { + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_SHUTTER"); + psArrayAdd (job->args, 1, image); + psArrayAdd (job->args, 1, shutterImage); + psArrayAdd (job->args, 1, mask); + PS_ARRAY_ADD_SCALAR (job->args, exptime, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR (job->args, blank, PS_TYPE_MASK); + PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32); + + // ppImageDetrendReadout(config, options, view) + if (!psThreadJobAddPending (job)) { + psFree(job); + return false; + } + psFree(job); + } else { + pmShutterCorrectionApplyScan (image, shutterImage, mask, exptime, blank, rowStart, rowStop); + } + } + if (threaded) { + // wait here for the threaded jobs to finish + if (!psThreadPoolWait(false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + fprintf (stderr, "success for threaded jobs\n"); + + // free the done jobs + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + psFree (job); + } + } } psFree(shutterImage); @@ -952,7 +952,7 @@ stdev->data.F32[stdev->n] = psStatsGetValue(stats, stdevStat) * refValue; - psTrace("psModules.detrend", 5, "input shutter image sample value %d: %f +/- %f -> %f +/- %f\n", j, - psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat), - mean->data.F32[mean->n], stdev->data.F32[stdev->n]); + psTrace("psModules.detrend", 5, "input shutter image sample value %d: %f +/- %f -> %f +/- %f\n", j, + psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat), + mean->data.F32[mean->n], stdev->data.F32[stdev->n]); data->mean->data[j] = psVectorExtend(mean, IMAGES_BUFFER, 1); @@ -978,5 +978,5 @@ // reshuffle exptimes to new sequence (this is only a local value) for (int j = 0; j < newtimes->n; j++) { - newtimes->data.F32[j] = data->exptimes->data.F32[index->data.S32[j]]; + newtimes->data.F32[j] = data->exptimes->data.F32[index->data.S32[j]]; } @@ -984,20 +984,20 @@ int numGood = 0; // Number of good measurements for (int i = 0; i < MEASURE_SAMPLES; i++) { - psVector *newcounts = psVectorAlloc (data->exptimes->n, PS_TYPE_F32); - psVector *newerrors = psVectorAlloc (data->exptimes->n, PS_TYPE_F32); - psVector *counts = data->mean->data[i]; - psVector *errors = data->stdev->data[i]; - - for (int j = 0; j < newcounts->n; j++) { - newcounts->data.F32[j] = counts->data.F32[index->data.S32[j]]; - newerrors->data.F32[j] = errors->data.F32[index->data.S32[j]]; - } - - // use the sorted exptime, counts, and errors for the measurements + psVector *newcounts = psVectorAlloc (data->exptimes->n, PS_TYPE_F32); + psVector *newerrors = psVectorAlloc (data->exptimes->n, PS_TYPE_F32); + psVector *counts = data->mean->data[i]; + psVector *errors = data->stdev->data[i]; + + for (int j = 0; j < newcounts->n; j++) { + newcounts->data.F32[j] = counts->data.F32[index->data.S32[j]]; + newerrors->data.F32[j] = errors->data.F32[index->data.S32[j]]; + } + + // use the sorted exptime, counts, and errors for the measurements pmShutterCorrection *guess = pmShutterCorrectionGuess(newtimes, newcounts); // Guess at correction - psTrace("psModules.detrend", 5, "Shutter correction guess: scale: %f, offset: %f, offref: %f\n", guess->scale, guess->offset, guess->offref); + psTrace("psModules.detrend", 5, "Shutter correction guess: scale: %f, offset: %f, offref: %f\n", guess->scale, guess->offset, guess->offref); pmShutterCorrection *corr = pmShutterCorrectionFullFit(newtimes, newcounts, newerrors, guess); // The actual correction - psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref); + psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref); if (corr && isfinite(corr->offref) && corr->valid) { @@ -1009,6 +1009,6 @@ psFree(corr); psFree(guess); - psFree (newcounts); - psFree (newerrors); + psFree (newcounts); + psFree (newerrors); } psFree (newtimes); @@ -1040,5 +1040,5 @@ pmReadoutStackDefineOutput(shutter, col0, row0, numCols, numRows, false, false, maskVal); if (pattern) { - pmReadoutStackDefineOutput(pattern, col0, row0, numCols, numRows, false, false, maskVal); + pmReadoutStackDefineOutput(pattern, col0, row0, numCols, numRows, false, false, maskVal); } @@ -1072,9 +1072,9 @@ pattern->data_exists = true; - if (pattern->parent) { - pattern->parent->data_exists = true; - if (pattern->parent->parent) { - pattern->parent->parent->data_exists = true; - } + if (pattern->parent) { + pattern->parent->data_exists = true; + if (pattern->parent->parent) { + pattern->parent->parent->data_exists = true; + } } @@ -1137,5 +1137,5 @@ errors->data.F32[r] = sqrtf(readout->weight->data.F32[yIn][xIn]) * ref; } else { - // XXX guess that the input data is Poisson distributed; if we go negative, force high + // XXX guess that the input data is Poisson distributed; if we go negative, force high errors->data.F32[r] = sqrtf(fabs(image->data.F32[yIn][xIn])) * ref; } Index: /trunk/psModules/src/detrend/pmShutterCorrection.h =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 18959) +++ /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 18960) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-05 01:24:47 $ + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-08 18:09:07 $ * Copyright 2006 Institute for Astronomy, University of Hawaii */ @@ -60,10 +60,10 @@ /// Shutter correction parameters, applicable for a single pixel typedef struct { - double scale; ///< The normalisation for an exposure, A(k) or f'(x,y) + double scale; ///< The normalisation for an exposure, A(k) or f'(x,y) double offset; ///< The time offset, dTk double offref; ///< The reference time offset, dTo int num; ///< Number of points used float stdev; ///< Standard deviation - bool valid; // is the fitted shutter correction valid (produce a significant improvement?) + bool valid; // is the fitted shutter correction valid (produce a significant improvement?) } pmShutterCorrection; @@ -124,5 +124,5 @@ ); -bool pmShutterCorrectionApplyScan_Threaded (psThreadJob *job); +bool pmShutterCorrectionApplyScan_Threaded(const psThreadJob *job); bool pmShutterCorrectionApplyScan(psImage *image, psImage *shutterImage, psImage *mask, float exptime, psMaskType blank, int rowStart, int rowStop);