IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18960


Ignore:
Timestamp:
Aug 8, 2008, 8:09:07 AM (18 years ago)
Author:
Paul Price
Message:

Fixing to match slight differences in psThread

Location:
trunk/psModules/src/detrend
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/detrend/pmBias.c

    r18893 r18960  
    1919#include "pmBias.h"
    2020
    21 bool pmBiasSubtractScan_Threaded (psThreadJob *job) {
     21bool pmBiasSubtractScan_Threaded(const psThreadJob *job)
     22{
     23    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
    2224    pmReadout *in = job->args->data[0];
    2325    pmReadout *sub = job->args->data[1];
    … …  
    6365// image.  The bias image is subtracted in place from the input image.
    6466bool pmBiasSubtractFrame(pmReadout *in, // Input readout
    65                         pmReadout *sub, // Readout to be subtracted from input
    66                         float scale   // Scale to apply before subtracting
     67                        pmReadout *sub, // Readout to be subtracted from input
     68                        float scale   // Scale to apply before subtracting
    6769    )
    6870{
    … …  
    121123    bool threaded = true;
    122124    int scanRows = pmDetrendGetScanRows();
    123     if (scanRows == 0) { 
    124         threaded = false;
    125         scanRows = inImage->numRows;
     125    if (scanRows == 0) {
     126        threaded = false;
     127        scanRows = inImage->numRows;
    126128    }
    127129
    … …  
    129131      int rowStop = PS_MIN (rowStart + scanRows, inImage->numRows);
    130132
    131 # define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \
    132       psScalar *scalar = psScalarAlloc(VALUE, TYPE); \
    133       psArrayAdd(ARRAY, 1, scalar); \
    134       psFree (scalar); }
    135      
    136133      if (threaded) {
    137           // allocate a job, construct the arguments for this job
    138           psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_BIAS");
    139           psArrayAdd (job->args, 1, in);
    140           psArrayAdd (job->args, 1, sub);
    141           PS_ARRAY_ADD_SCALAR (job->args, scale, PS_TYPE_F32);
    142           PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32);
    143           PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32);
    144           PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
    145           PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
    146 
    147           // ppImageDetrendReadout(config, options, view)
    148           if (!psThreadJobAddPending (job)) {
    149               return false;
    150           }
     134          // allocate a job, construct the arguments for this job
     135          psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_BIAS");
     136          psArrayAdd (job->args, 1, in);
     137          psArrayAdd (job->args, 1, sub);
     138          PS_ARRAY_ADD_SCALAR (job->args, scale, PS_TYPE_F32);
     139          PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32);
     140          PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32);
     141          PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
     142          PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
     143
     144          // ppImageDetrendReadout(config, options, view)
     145          if (!psThreadJobAddPending (job)) {
     146              psFree(job);
     147              return false;
     148          }
     149          psFree(job);
    151150      } else {
    152           pmBiasSubtractScan (in, sub, scale, xOffset, yOffset, rowStart, rowStop);
     151          pmBiasSubtractScan (in, sub, scale, xOffset, yOffset, rowStart, rowStop);
    153152      }
    154153    }
    155154
    156155    if (threaded) {
    157         // wait here for the threaded jobs to finish
    158         if (!psThreadPoolWait ()) {
    159             psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    160             return false;
    161         }
    162         fprintf (stderr, "success for threaded jobs\n");
    163 
    164         // each job records its own goodPixel values; sum them here
    165         // we have only supplied one type of job, so we can assume the types here
    166         psThreadJob *job = NULL;
    167         while ((job = psThreadJobGetDone()) != NULL) {
    168             psFree (job);
    169         }
     156        // wait here for the threaded jobs to finish
     157        if (!psThreadPoolWait(false)) {
     158            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     159            return false;
     160        }
     161        fprintf (stderr, "success for threaded jobs\n");
     162
     163        // each job records its own goodPixel values; sum them here
     164        // we have only supplied one type of job, so we can assume the types here
     165        psThreadJob *job = NULL;
     166        while ((job = psThreadJobGetDone()) != NULL) {
     167            psFree (job);
     168        }
    170169    }
    171170

  • trunk/psModules/src/detrend/pmBias.h

    r18893 r18960  
    55 * @author Paul Price, IfA
    66 *
    7  * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2008-08-05 01:24:47 $
     7 * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2008-08-08 18:09:07 $
    99 * Copyright 2004--2006 Institute for Astronomy, University of Hawaii
    1010 */
    … …  
    3030// image.  The bias image is subtracted in place from the input image.
    3131bool pmBiasSubtractFrame(pmReadout *in, // Input readout
    32                         pmReadout *sub, // Readout to be subtracted from input
    33                         float scale   // Scale to apply before subtracting
     32                        pmReadout *sub, // Readout to be subtracted from input
     33                        float scale   // Scale to apply before subtracting
    3434    );
    3535
    36 bool pmBiasSubtractScan_Threaded (psThreadJob *job);
     36bool pmBiasSubtractScan_Threaded(const psThreadJob *job);
    3737bool pmBiasSubtractScan (pmReadout *in, pmReadout *sub, float scan, int xOffset, int yOffset, int rowStart, int rowStop);
    3838

  • trunk/psModules/src/detrend/pmDark.c

    r18893 r18960  
    102102bool pmDarkCombinePrepare(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept)
    103103{
    104     psArray *values = psArrayAlloc(inputs->n); 
     104    psArray *values = psArrayAlloc(inputs->n);
    105105    psVector *roMask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for bad readouts
    106106    psVector *norm = normConcept ? psVectorAlloc(inputs->n, PS_TYPE_F32) : NULL; // Normalizations for each
    … …  
    117117        if (!norm) continue;
    118118
    119         pmReadout *readout = inputs->data[i]; // Readout of interest
    120         float normValue;            // Normalisation value
    121         if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) {
    122             psWarning("Unable to find value of %s for readout %d", normConcept, i);
    123             roMask->data.U8[i] = 0xff;
    124             norm->data.F32[i] = NAN;
    125             numBadInputs++;
    126             continue;
    127         }
    128         if (normValue == 0.0) {
    129             psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i);
    130             roMask->data.U8[i] = 0xff;
    131             norm->data.F32[i] = NAN;
    132             numBadInputs++;
    133             continue;
    134         }
    135         norm->data.F32[i] = 1.0 / normValue;
     119        pmReadout *readout = inputs->data[i]; // Readout of interest
     120        float normValue;            // Normalisation value
     121        if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) {
     122            psWarning("Unable to find value of %s for readout %d", normConcept, i);
     123            roMask->data.U8[i] = 0xff;
     124            norm->data.F32[i] = NAN;
     125            numBadInputs++;
     126            continue;
     127        }
     128        if (normValue == 0.0) {
     129            psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i);
     130            roMask->data.U8[i] = 0xff;
     131            norm->data.F32[i] = NAN;
     132            numBadInputs++;
     133            continue;
     134        }
     135        norm->data.F32[i] = 1.0 / normValue;
    136136    }
    137137
    … …  
    144144            psFree(roMask);
    145145            psFree(orders);
    146             psFree(norm);
     146            psFree(norm);
    147147            return false;
    148148        }
    … …  
    215215        }
    216216
    217         pmReadoutStackDefineOutput(readout, col0, row0, numCols, numRows, false, false, 0);
    218         psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", col0, row0);
     217        pmReadoutStackDefineOutput(readout, col0, row0, numCols, numRows, false, false, 0);
     218        psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", col0, row0);
    219219    }
    220220
    … …  
    269269    bool mdok = false;
    270270
    271     // retrieve the required parameter vectors 
    272     psArray *values      = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");  psAssert (values, "values not supplied");
    273     psVector *roMask     = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK"); psAssert (roMask, "roMask not supplied");
     271    // retrieve the required parameter vectors
     272    psArray *values      = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");  psAssert (values, "values not supplied");
     273    psVector *roMask     = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK"); psAssert (roMask, "roMask not supplied");
    274274    psPolynomialMD *poly = psMetadataLookupPtr(&mdok, output->analysis, "DARK.POLY");    psAssert (poly, "orders not supplied");
    275275
    276276    // retrieve the norm vector, if supplied
    277     psVector *norm       = psMetadataLookupPtr(&mdok, output->analysis, "DARK.NORM");   
     277    psVector *norm       = psMetadataLookupPtr(&mdok, output->analysis, "DARK.NORM");
    278278
    279279    // retrieve the 'counts' and 'sigma' images
    … …  
    367367}
    368368
    369 bool pmDarkApplyScan_Threaded (psThreadJob *job) {
     369bool pmDarkApplyScan_Threaded(const psThreadJob *job)
     370{
     371    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     372
    370373    pmReadout *readout   = job->args->data[0];
    371374    pmCell *dark         = job->args->data[1];
    … …  
    478481    bool threaded = true;
    479482    int scanRows = pmDetrendGetScanRows();
    480     if (scanRows == 0) { 
    481         threaded = false;
    482         scanRows = readout->image->numRows;
     483    if (scanRows == 0) {
     484        threaded = false;
     485        scanRows = readout->image->numRows;
    483486    }
    484487
    … …  
    486489      int rowStop = PS_MIN (rowStart + scanRows, readout->image->numRows);
    487490
    488 # define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \
    489       psScalar *scalar = psScalarAlloc(VALUE, TYPE); \
    490       psArrayAdd(ARRAY, 1, scalar); \
    491       psFree (scalar); }
    492      
    493491      if (threaded) {
    494           // allocate a job, construct the arguments for this job
    495           psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_DARK");
    496           psArrayAdd (job->args, 1, readout);
    497           psArrayAdd (job->args, 1, dark);
    498           psArrayAdd (job->args, 1, poly);
    499           psArrayAdd (job->args, 1, values);
    500           PS_ARRAY_ADD_SCALAR (job->args, bad, PS_TYPE_MASK);
    501           PS_ARRAY_ADD_SCALAR (job->args, doNorm, PS_TYPE_U8);
    502           PS_ARRAY_ADD_SCALAR (job->args, norm, PS_TYPE_F32);
    503           PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
    504           PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
    505 
    506           // ppImageDetrendReadout(config, options, view)
    507           if (!psThreadJobAddPending (job)) {
    508               return false;
    509           }
     492          // allocate a job, construct the arguments for this job
     493          psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_DARK");
     494          psArrayAdd (job->args, 1, readout);
     495          psArrayAdd (job->args, 1, dark);
     496          psArrayAdd (job->args, 1, poly);
     497          psArrayAdd (job->args, 1, values);
     498          PS_ARRAY_ADD_SCALAR (job->args, bad, PS_TYPE_MASK);
     499          PS_ARRAY_ADD_SCALAR (job->args, doNorm, PS_TYPE_U8);
     500          PS_ARRAY_ADD_SCALAR (job->args, norm, PS_TYPE_F32);
     501          PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
     502          PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
     503
     504          // ppImageDetrendReadout(config, options, view)
     505          if (!psThreadJobAddPending (job)) {
     506              psFree(job);
     507              return false;
     508          }
     509          psFree(job);
    510510      } else {
    511           pmDarkApplyScan (readout, dark, poly, values, bad, doNorm, norm, rowStart, rowStop);
     511          pmDarkApplyScan (readout, dark, poly, values, bad, doNorm, norm, rowStart, rowStop);
    512512      }
    513513    }
    514514
    515515    if (threaded) {
    516         // wait here for the threaded jobs to finish
    517         if (!psThreadPoolWait ()) {
    518             psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    519             return false;
    520         }
    521         fprintf (stderr, "success for threaded jobs\n");
    522 
    523         // free the done jobs
    524         psThreadJob *job = NULL;
    525         while ((job = psThreadJobGetDone()) != NULL) {
    526             psFree (job);
    527         }
     516        // wait here for the threaded jobs to finish
     517        if (!psThreadPoolWait(false)) {
     518            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     519            return false;
     520        }
     521        fprintf (stderr, "success for threaded jobs\n");
     522
     523        // free the done jobs
     524        psThreadJob *job = NULL;
     525        while ((job = psThreadJobGetDone()) != NULL) {
     526            psFree (job);
     527        }
    528528    }
    529529

  • trunk/psModules/src/detrend/pmDark.h

    r18893 r18960  
    2727// Combine darks -- preparation step
    2828bool pmDarkCombinePrepare(pmCell *output,      // Output cell; readouts will be attached
    29                           const psArray *inputs, // Input readouts for combination
    30                           psArray *ordinates,  // Ordinates for fitting
    31                           const char *normConcept // Concept name to use to divide input pixel values
     29                          const psArray *inputs, // Input readouts for combination
     30                          psArray *ordinates,  // Ordinates for fitting
     31                          const char *normConcept // Concept name to use to divide input pixel values
    3232    );
    3333
    … …  
    4040    );
    4141
    42 bool pmDarkApplyScan_Threaded (psThreadJob *job);
     42bool pmDarkApplyScan_Threaded(const psThreadJob *job);
    4343bool pmDarkApplyScan (pmReadout *readout, const pmCell *dark, psPolynomialMD *poly, psVector *values, psMaskType bad, bool doNorm, float norm, int rowStart, int rowStop);
    4444

  • trunk/psModules/src/detrend/pmDetrendThreads.c

    r18893 r18960  
    1818#include "pmDetrendThreads.h"
    1919
    20 static int scanRows = 0;
     20static int scanRows = 0;                // Number of rows to work on at once
    2121
    22 int pmDetrendGetScanRows () {
     22int pmDetrendGetScanRows(void)
     23{
    2324    return scanRows;
    24 }   
     25}
    2526
    26 bool pmDetrendSetThreadTasks (int newScanRows) {
     27bool pmDetrendSetThreadTasks (int newScanRows)
     28{
     29    psAssert(scanRows == 0, "programming error: program called pmDetrendSetThreadTasks twice");
    2730
    28     psAssert (scanRows == 0, "programming error: program called pmDetrendSetThreadTasks twice");
     31    PS_ASSERT_INT_POSITIVE(newScanRows, false);
    2932    scanRows = newScanRows;
    3033
    31     psThreadTask *task = NULL;
     34    {
     35        psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_BIAS", 7);
     36        task->function = &pmBiasSubtractScan_Threaded;
     37        psThreadTaskAdd(task);
     38        psFree(task);
     39    }
    3240
    33     task = psThreadTaskAlloc ("PSMODULES_DETREND_BIAS", 7);
    34     task->function = &pmBiasSubtractScan_Threaded;
    35     psThreadTaskAdd (task);
     41    {
     42        psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_DARK", 9);
     43        task->function = &pmDarkApplyScan_Threaded;
     44        psThreadTaskAdd(task);
     45        psFree(task);
     46    }
    3647
    37     task = psThreadTaskAlloc ("PSMODULES_DETREND_DARK", 9);
    38     task->function = &pmDarkApplyScan_Threaded;
    39     psThreadTaskAdd (task);
     48    {
     49        psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_SHUTTER", 7);
     50        task->function = &pmShutterCorrectionApplyScan_Threaded;
     51        psThreadTaskAdd(task);
     52        psFree(task);
     53    }
    4054
    41     task = psThreadTaskAlloc ("PSMODULES_DETREND_SHUTTER", 7);
    42     task->function = &pmShutterCorrectionApplyScan_Threaded;
    43     psThreadTaskAdd (task);
    44 
    45     task = psThreadTaskAlloc ("PSMODULES_DETREND_FLAT", 9);
    46     task->function = &pmFlatFieldScan_Threaded;
    47     psThreadTaskAdd (task);
     55    {
     56        psThreadTask *task = psThreadTaskAlloc("PSMODULES_DETREND_FLAT", 9);
     57        task->function = &pmFlatFieldScan_Threaded;
     58        psThreadTaskAdd(task);
     59        psFree(task);
     60    }
    4861
    4962    return true;

  • trunk/psModules/src/detrend/pmFlatField.c

    r18893 r18960  
    1313#include "pmDetrendThreads.h"
    1414
    15 bool pmFlatFieldScan_Threaded (psThreadJob *job) {
     15bool pmFlatFieldScan_Threaded(const psThreadJob *job)
     16{
     17    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     18
    1619    psImage *inImage   = job->args->data[0];
    1720    psImage *inMask    = job->args->data[1];
    … …  
    2932
    3033// Macro for all PS types
    31 #define FLAT_DIVISION_CASE(TYPE, SPECIAL)       \
    32     case PS_TYPE_##TYPE:                        \
     34#define FLAT_DIVISION_CASE(TYPE, SPECIAL)       \
     35    case PS_TYPE_##TYPE:                        \
    3336    for (int j = rowStart; j < rowStop; j++) { \
    3437        for (int i = 0; i < inImage->numCols; i++) { \
    … …  
    129132    bool threaded = true;
    130133    int scanRows = pmDetrendGetScanRows();
    131     if (scanRows == 0) { 
    132         threaded = false;
    133         scanRows = inImage->numRows;
     134    if (scanRows == 0) {
     135        threaded = false;
     136        scanRows = inImage->numRows;
    134137    }
    135138
    … …  
    137140      int rowStop = PS_MIN (rowStart + scanRows, inImage->numRows);
    138141
    139 # define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \
    140       psScalar *scalar = psScalarAlloc(VALUE, TYPE); \
    141       psArrayAdd(ARRAY, 1, scalar); \
    142       psFree (scalar); }
    143      
    144142      if (threaded) {
    145           // allocate a job, construct the arguments for this job
    146           psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_FLAT");
    147           psArrayAdd (job->args, 1, inImage);
    148           psArrayAdd (job->args, 1, inMask);
    149           psArrayAdd (job->args, 1, flatImage);
    150           psArrayAdd (job->args, 1, flatMask);
    151           PS_ARRAY_ADD_SCALAR (job->args, badFlat, PS_TYPE_U8);
    152           PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32);
    153           PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32);
    154           PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
    155           PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
     143          // allocate a job, construct the arguments for this job
     144          psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_FLAT");
     145          psArrayAdd (job->args, 1, inImage);
     146          psArrayAdd (job->args, 1, inMask);
     147          psArrayAdd (job->args, 1, flatImage);
     148          psArrayAdd (job->args, 1, flatMask);
     149          PS_ARRAY_ADD_SCALAR (job->args, badFlat, PS_TYPE_U8);
     150          PS_ARRAY_ADD_SCALAR (job->args, xOffset, PS_TYPE_S32);
     151          PS_ARRAY_ADD_SCALAR (job->args, yOffset, PS_TYPE_S32);
     152          PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
     153          PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
    156154
    157           // ppImageDetrendReadout(config, options, view)
    158           if (!psThreadJobAddPending (job)) {
    159               return false;
    160           }
     155          // ppImageDetrendReadout(config, options, view)
     156          if (!psThreadJobAddPending (job)) {
     157              psFree(job);
     158              return false;
     159          }
     160          psFree(job);
    161161      } else {
    162           pmFlatFieldScan (inImage, inMask, flatImage, flatMask, badFlat, xOffset, yOffset, rowStart, rowStop);
     162          pmFlatFieldScan (inImage, inMask, flatImage, flatMask, badFlat, xOffset, yOffset, rowStart, rowStop);
    163163      }
    164164    }
    165165
    166166    if (threaded) {
    167         // wait here for the threaded jobs to finish
    168         if (!psThreadPoolWait ()) {
    169             psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    170             return false;
    171         }
    172         fprintf (stderr, "success for threaded jobs\n");
     167        // wait here for the threaded jobs to finish
     168        if (!psThreadPoolWait(false)) {
     169            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     170            return false;
     171        }
     172        fprintf (stderr, "success for threaded jobs\n");
    173173
    174         // free done jobs
    175         psThreadJob *job = NULL;
    176         while ((job = psThreadJobGetDone()) != NULL) {
    177             psFree (job);
    178         }
     174        // free done jobs
     175        psThreadJob *job = NULL;
     176        while ((job = psThreadJobGetDone()) != NULL) {
     177            psFree (job);
     178        }
    179179    }
    180180

  • trunk/psModules/src/detrend/pmFlatField.h

    r18893 r18960  
    55 * @author Paul Price, IfA
    66 *
    7  * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2008-08-05 01:24:47 $
     7 * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2008-08-08 18:09:07 $
    99 * Copyright 2004-2006 Institute for Astronomy, University of Hawaii
    1010 */
    … …  
    2828                );
    2929
    30 bool pmFlatFieldScan_Threaded (psThreadJob *job);
     30bool pmFlatFieldScan_Threaded(const psThreadJob *job);
    3131bool pmFlatFieldScan (psImage *inImage, psImage *inMask, psImage *flatImage, psImage *flatMask, psMaskType badFlag, int xOffset, int yOffset, int rowStart, int rowStop);
    3232

  • trunk/psModules/src/detrend/pmShutterCorrection.c

    r18893 r18960  
    122122
    123123    // Iterate only
    124     for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++); 
     124    for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++);
    125125
    126126    if (index == N - 1) {
    … …  
    344344    psVector *resid = psVectorAlloc (exptime->n, PS_TYPE_F32);
    345345    for (int i = 0; i < exptime->n; i++) {
    346         float fitCounts = corr->scale * (exptime->data.F32[i] + corr->offset) / (exptime->data.F32[i] + corr->offref);
    347         resid->data.F32[i] = counts->data.F32[i] - fitCounts;
    348     }
    349    
     346        float fitCounts = corr->scale * (exptime->data.F32[i] + corr->offset) / (exptime->data.F32[i] + corr->offref);
     347        resid->data.F32[i] = counts->data.F32[i] - fitCounts;
     348    }
     349
    350350    psStats *rawStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    351351    psStats *resStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    … …  
    353353    psVectorStats (resStats, resid, NULL, NULL, 0);
    354354
    355     // XXX temporary hard-wired minimum stdev improvement factor 
     355    // XXX temporary hard-wired minimum stdev improvement factor
    356356    psTrace("psModules.detrend", 3, "raw scatter %f vs res scatter %f\n", rawStats->sampleStdev, resStats->sampleStdev);
    357357    if (rawStats->sampleStdev / resStats->sampleStdev < 1.5) corr->valid = false;
    … …  
    433433            numRows = image->numRows;
    434434            numCols = image->numCols;
    435             // define the reference region : a box of size 'size' at the center
     435            // define the reference region : a box of size 'size' at the center
    436436            refRegion = psRegionForSquare(0.5 * numCols, 0.5 * numRows, size);
    437437            // Set up the sample regions : boxes of size 'size' at the 4 image corners
    … …  
    656656
    657657
    658 bool pmShutterCorrectionApplyScan_Threaded (psThreadJob *job) {
     658bool pmShutterCorrectionApplyScan_Threaded(const psThreadJob *job)
     659{
     660    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     661
    659662    psImage *image        = job->args->data[0];
    660663    psImage *shutterImage = job->args->data[1];
    … …  
    672675{
    673676    for (int y = rowStart; y < rowStop; y++) {
    674         for (int x = 0; x < image->numCols; x++) {
    675             if (mask && !isfinite(shutterImage->data.F32[y][x])) {
    676                 mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
    677                 image->data.F32[y][x] = NAN;
    678                 continue;
    679             }
    680             image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]);
    681         }
     677        for (int x = 0; x < image->numCols; x++) {
     678            if (mask && !isfinite(shutterImage->data.F32[y][x])) {
     679                mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
     680                image->data.F32[y][x] = NAN;
     681                continue;
     682            }
     683            image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]);
     684        }
    682685    }
    683686    return true;
    … …  
    732735    bool threaded = true;
    733736    int scanRows = pmDetrendGetScanRows();
    734     if (scanRows == 0) { 
    735         threaded = false;
    736         scanRows = image->numRows;
     737    if (scanRows == 0) {
     738        threaded = false;
     739        scanRows = image->numRows;
    737740    }
    738741
    … …  
    756759        psFree (line);
    757760    } else {
    758         for (int rowStart = 0; rowStart < image->numRows; rowStart += scanRows) {
    759             int rowStop = PS_MIN (rowStart + scanRows, image->numRows);
    760 
    761 # define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \
    762       psScalar *scalar = psScalarAlloc(VALUE, TYPE); \
    763       psArrayAdd(ARRAY, 1, scalar); \
    764       psFree (scalar); }
    765      
    766             if (threaded) {
    767                 // allocate a job, construct the arguments for this job
    768                 psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_SHUTTER");
    769                 psArrayAdd (job->args, 1, image);
    770                 psArrayAdd (job->args, 1, shutterImage);
    771                 psArrayAdd (job->args, 1, mask);
    772                 PS_ARRAY_ADD_SCALAR (job->args, exptime, PS_TYPE_F32);
    773                 PS_ARRAY_ADD_SCALAR (job->args, blank, PS_TYPE_MASK);
    774                 PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
    775                 PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
    776 
    777                 // ppImageDetrendReadout(config, options, view)
    778                 if (!psThreadJobAddPending (job)) {
    779                     return false;
    780                 }
    781             } else {
    782                 pmShutterCorrectionApplyScan (image, shutterImage, mask, exptime, blank, rowStart, rowStop);
    783             }
    784         }
    785         if (threaded) {
    786             // wait here for the threaded jobs to finish
    787             if (!psThreadPoolWait ()) {
    788                 psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    789                 return false;
    790             }
    791             fprintf (stderr, "success for threaded jobs\n");
    792 
    793             // free the done jobs
    794             psThreadJob *job = NULL;
    795             while ((job = psThreadJobGetDone()) != NULL) {
    796                 psFree (job);
    797             }
    798         }
     761        for (int rowStart = 0; rowStart < image->numRows; rowStart += scanRows) {
     762            int rowStop = PS_MIN (rowStart + scanRows, image->numRows);
     763
     764            if (threaded) {
     765                // allocate a job, construct the arguments for this job
     766                psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_SHUTTER");
     767                psArrayAdd (job->args, 1, image);
     768                psArrayAdd (job->args, 1, shutterImage);
     769                psArrayAdd (job->args, 1, mask);
     770                PS_ARRAY_ADD_SCALAR (job->args, exptime, PS_TYPE_F32);
     771                PS_ARRAY_ADD_SCALAR (job->args, blank, PS_TYPE_MASK);
     772                PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
     773                PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
     774
     775                // ppImageDetrendReadout(config, options, view)
     776                if (!psThreadJobAddPending (job)) {
     777                    psFree(job);
     778                    return false;
     779                }
     780                psFree(job);
     781            } else {
     782                pmShutterCorrectionApplyScan (image, shutterImage, mask, exptime, blank, rowStart, rowStop);
     783            }
     784        }
     785        if (threaded) {
     786            // wait here for the threaded jobs to finish
     787            if (!psThreadPoolWait(false)) {
     788                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     789                return false;
     790            }
     791            fprintf (stderr, "success for threaded jobs\n");
     792
     793            // free the done jobs
     794            psThreadJob *job = NULL;
     795            while ((job = psThreadJobGetDone()) != NULL) {
     796                psFree (job);
     797            }
     798        }
    799799    }
    800800    psFree(shutterImage);
    … …  
    952952        stdev->data.F32[stdev->n] = psStatsGetValue(stats, stdevStat) * refValue;
    953953
    954         psTrace("psModules.detrend", 5, "input shutter image sample value %d: %f +/- %f  ->  %f +/- %f\n", j,
    955                 psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat),
    956                 mean->data.F32[mean->n], stdev->data.F32[stdev->n]);
     954        psTrace("psModules.detrend", 5, "input shutter image sample value %d: %f +/- %f  ->  %f +/- %f\n", j,
     955                psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat),
     956                mean->data.F32[mean->n], stdev->data.F32[stdev->n]);
    957957
    958958        data->mean->data[j] = psVectorExtend(mean, IMAGES_BUFFER, 1);
    … …  
    978978    // reshuffle exptimes to new sequence (this is only a local value)
    979979    for (int j = 0; j < newtimes->n; j++) {
    980         newtimes->data.F32[j] = data->exptimes->data.F32[index->data.S32[j]];
     980        newtimes->data.F32[j] = data->exptimes->data.F32[index->data.S32[j]];
    981981    }
    982982
    … …  
    984984    int numGood = 0;                    // Number of good measurements
    985985    for (int i = 0; i < MEASURE_SAMPLES; i++) {
    986         psVector *newcounts = psVectorAlloc (data->exptimes->n, PS_TYPE_F32);
    987         psVector *newerrors = psVectorAlloc (data->exptimes->n, PS_TYPE_F32);
    988         psVector *counts = data->mean->data[i];
    989         psVector *errors = data->stdev->data[i];
    990 
    991         for (int j = 0; j < newcounts->n; j++) {
    992             newcounts->data.F32[j] = counts->data.F32[index->data.S32[j]];
    993             newerrors->data.F32[j] = errors->data.F32[index->data.S32[j]];
    994         }
    995 
    996         // use the sorted exptime, counts, and errors for the measurements
     986        psVector *newcounts = psVectorAlloc (data->exptimes->n, PS_TYPE_F32);
     987        psVector *newerrors = psVectorAlloc (data->exptimes->n, PS_TYPE_F32);
     988        psVector *counts = data->mean->data[i];
     989        psVector *errors = data->stdev->data[i];
     990
     991        for (int j = 0; j < newcounts->n; j++) {
     992            newcounts->data.F32[j] = counts->data.F32[index->data.S32[j]];
     993            newerrors->data.F32[j] = errors->data.F32[index->data.S32[j]];
     994        }
     995
     996        // use the sorted exptime, counts, and errors for the measurements
    997997        pmShutterCorrection *guess = pmShutterCorrectionGuess(newtimes, newcounts); // Guess at correction
    998         psTrace("psModules.detrend", 5, "Shutter correction guess: scale: %f, offset: %f, offref: %f\n", guess->scale, guess->offset, guess->offref);
     998        psTrace("psModules.detrend", 5, "Shutter correction guess: scale: %f, offset: %f, offref: %f\n", guess->scale, guess->offset, guess->offref);
    999999
    10001000        pmShutterCorrection *corr = pmShutterCorrectionFullFit(newtimes, newcounts, newerrors, guess); // The actual correction
    1001         psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref);
     1001        psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref);
    10021002
    10031003        if (corr && isfinite(corr->offref) && corr->valid) {
    … …  
    10091009        psFree(corr);
    10101010        psFree(guess);
    1011         psFree (newcounts);
    1012         psFree (newerrors);
     1011        psFree (newcounts);
     1012        psFree (newerrors);
    10131013    }
    10141014    psFree (newtimes);
    … …  
    10401040    pmReadoutStackDefineOutput(shutter, col0, row0, numCols, numRows, false, false, maskVal);
    10411041    if (pattern) {
    1042         pmReadoutStackDefineOutput(pattern, col0, row0, numCols, numRows, false, false, maskVal);
     1042        pmReadoutStackDefineOutput(pattern, col0, row0, numCols, numRows, false, false, maskVal);
    10431043    }
    10441044
    … …  
    10721072
    10731073    pattern->data_exists = true;
    1074     if (pattern->parent) { 
    1075         pattern->parent->data_exists = true;
    1076         if (pattern->parent->parent) {
    1077             pattern->parent->parent->data_exists = true;
    1078         }
     1074    if (pattern->parent) {
     1075        pattern->parent->data_exists = true;
     1076        if (pattern->parent->parent) {
     1077            pattern->parent->parent->data_exists = true;
     1078        }
    10791079    }
    10801080
    … …  
    11371137                    errors->data.F32[r] = sqrtf(readout->weight->data.F32[yIn][xIn]) * ref;
    11381138                } else {
    1139                     // XXX guess that the input data is Poisson distributed; if we go negative, force high
     1139                    // XXX guess that the input data is Poisson distributed; if we go negative, force high
    11401140                    errors->data.F32[r] = sqrtf(fabs(image->data.F32[yIn][xIn])) * ref;
    11411141                }

  • trunk/psModules/src/detrend/pmShutterCorrection.h

    r18893 r18960  
    55 * @author Paul Price, IfA
    66 *
    7  * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2008-08-05 01:24:47 $
     7 * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2008-08-08 18:09:07 $
    99 * Copyright 2006 Institute for Astronomy, University of Hawaii
    1010 */
    … …  
    6060/// Shutter correction parameters, applicable for a single pixel
    6161typedef struct {
    62     double scale;                       ///< The normalisation for an exposure, A(k) or f'(x,y) 
     62    double scale;                       ///< The normalisation for an exposure, A(k) or f'(x,y)
    6363    double offset;                      ///< The time offset, dTk
    6464    double offref;                      ///< The reference time offset, dTo
    6565    int num;                            ///< Number of points used
    6666    float stdev;                        ///< Standard deviation
    67     bool valid;                         // is the fitted shutter correction valid (produce a significant improvement?)
     67    bool valid;                         // is the fitted shutter correction valid (produce a significant improvement?)
    6868} pmShutterCorrection;
    6969
    … …  
    124124    );
    125125
    126 bool pmShutterCorrectionApplyScan_Threaded (psThreadJob *job);
     126bool pmShutterCorrectionApplyScan_Threaded(const psThreadJob *job);
    127127bool pmShutterCorrectionApplyScan(psImage *image, psImage *shutterImage, psImage *mask, float exptime, psMaskType blank, int rowStart, int rowStop);
    128128
Note: See TracChangeset for help on using the changeset viewer.