IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 25, 2009, 2:00:56 PM (17 years ago)
Author:
eugene
Message:

merging changes from head

Location:
branches/eam_branches/20090522
Files:
26 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090522

  • branches/eam_branches/20090522/psModules

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/eam_branches/20090522/psModules/src/camera/pmFPABin.c

    r24003 r24557  
    6565            }
    6666
    67             float imageValue, maskValue;// Values to set
     67            // Values to set
     68            float imageValue;
     69            psImageMaskType maskValue;
    6870            if (numPix > 0) {
    6971                imageValue = sum / numPix;

  • branches/eam_branches/20090522/psModules/src/camera/pmFPACopy.c

    r23761 r24557  
    7070    psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
    7171    *target = psImageAlloc(binning->nXruff, binning->nYruff, source->type.type);
     72    psImageInit (*target, 0.0);
    7273    return;
    7374}
    … …  
    221222    psTrace("psModules.camera", 3, "xFlip: %d; yFlip: %d\n", xFlip, yFlip);
    222223
     224    // Blow away extant readouts
     225    for (int i = 0; i < target->readouts->n; i++) {
     226        psFree(target->readouts->data[i]);
     227        target->readouts->data[i] = NULL;
     228    }
     229    target->readouts->n = 0;
     230
    223231    // Perform deep copy of the images.  I would prefer *not* to do a deep copy, in the interests of speed (we
    224232    // still need to do another deep copy into the HDU for when we write out), but this is the only way I can
    … …  
    242250        readoutCopyComponent(&targetReadout->variance, sourceReadout->variance, binning, xFlip, yFlip,
    243251                             pixels);
     252        // Copy covariance matrix: doesn't care about flips, etc.
     253        if (sourceReadout->covariance) {
     254            if (targetReadout->covariance) {
     255                psFree(targetReadout->covariance);
     256            }
     257            targetReadout->covariance = psKernelCopy(sourceReadout->covariance);
     258#if 0
     259            if (binning) {
     260                // XXX This isn't strictly correct, but we don't have a function that bins covariance matrices
     261                // with unequal binning factors.
     262                psKernel *covar = psImageCovarianceBin(PS_MAX(binning->nXbin, binning->nYbin),
     263                                                       targetReadout->covariance);
     264                psFree(targetReadout->covariance);
     265                targetReadout->covariance = covar;
     266            }
     267#endif
     268        }
    244269
    245270        // Copy bias

  • branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.c

    r23989 r24557  
    199199}
    200200
    201 bool pmReadoutSetVariance(pmReadout *readout, bool poisson)
     201bool pmReadoutSetVariance(pmReadout *readout, const psImage *noiseMap, bool poisson)
    202202{
    203203    PS_ASSERT_PTR_NON_NULL(readout, false);
     204    // check that the noiseMap (if it exists) matches the readout variance size)
    204205
    205206    pmCell *cell = readout->parent;     // The parent cell
    … …  
    228229
    229230        // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1.
     231        // XXX this calculation is wrong: limit is 1 e-, but this is in DN
    230232        readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs");
    231233        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max",
    … …  
    239241    }
    240242
    241     readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+",
    242                                            psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
     243    // apply a supplied readnoise map (NOTE: in DN, not electrons):
     244    if (noiseMap) {
     245        psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
     246        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
     247        psFree (rdVar);
     248    } else {
     249        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
     250    }
    243251
    244252    return true;
    … …  
    247255// this function creates the variance pixels, or uses the existing variance pixels.  it will set
    248256// the noise pixel values only if the variance image is not supplied
    249 bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson)
     257bool pmReadoutGenerateVariance(pmReadout *readout, const psImage *noiseMap, bool poisson)
    250258{
    251259    PS_ASSERT_PTR_NON_NULL(readout, false);
    … …  
    291299    readout->variance = variance;
    292300
    293     return pmReadoutSetVariance(readout, poisson);
    294 }
    295 
    296 bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
     301    return pmReadoutSetVariance(readout, noiseMap, poisson);
     302}
     303
     304bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson)
    297305{
    298306    PS_ASSERT_PTR_NON_NULL(readout, false);
    … …  
    301309
    302310    success &= pmReadoutGenerateMask(readout, satMask, badMask);
    303     success &= pmReadoutGenerateVariance(readout, poisson);
     311    success &= pmReadoutGenerateVariance(readout, noiseMap, poisson);
    304312
    305313    return success;
    306314}
    307315
    308 bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
     316bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, const psImage *noiseMap, bool poisson)
    309317{
    310318    PS_ASSERT_PTR_NON_NULL(cell, false);
    … …  
    314322    for (int i = 0; i < readouts->n; i++) {
    315323        pmReadout *readout = readouts->data[i]; // The readout
    316         success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask);
     324        success &= pmReadoutGenerateMaskVariance(readout, satMask, badMask, noiseMap, poisson);
    317325    }
    318326

  • branches/eam_branches/20090522/psModules/src/camera/pmFPAMaskWeight.h

    r21363 r24557  
    5454/// can't be generated.
    5555bool pmReadoutSetVariance(pmReadout *readout, ///< Readout for which to set variance
     56                          const psImage *noiseMap, ///< 2D image of the read noise in DN
    5657                          bool poisson    ///< Include poisson variance (in addition to read noise)?
    5758    );
    … …  
    7273/// with HDU entry).  This is intended for most operations.
    7374bool pmReadoutGenerateVariance(pmReadout *readout, ///< Readout for which to generate variance
     75                          const psImage *noiseMap, ///< 2D image of the read noise in DN
    7476                               bool poisson    ///< Include poisson variance (in addition to read noise)?
    7577    );
    … …  
    8183                                   psImageMaskType sat, ///< Mask value to give saturated pixels
    8284                                   psImageMaskType bad, ///< Mask value to give bad (low) pixels
     85                                   const psImage *noiseMap, ///< 2D image of the read noise in DN
    8386                                   bool poisson ///< Include poisson variance (in addition to read noise)?
    8487    );
    … …  
    9093                                psImageMaskType sat, ///< Mask value to give saturated pixels
    9194                                psImageMaskType bad, ///< Mask value to give bad (low) pixels
     95                                const psImage *noiseMap, ///< 2D image of the read noise in DN
    9296                                bool poisson ///< Include poisson variance (in addition to read noise)?
    9397    );

  • branches/eam_branches/20090522/psModules/src/camera/pmHDU.c

    r21363 r24557  
    9393    }
    9494
     95    psTrace("psModules.camera", 5, "Reading the header...\n");
     96
     97    // The header may already exist (e.g., from doing concept writing at the PHU level) so we need to be
     98    // careful.  We read into a separate container and copy that over the top of anything that's already read.
     99    psMetadata *header = psFitsReadHeader(NULL, fits);
     100    if (!header) {
     101        psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname);
     102        return false;
     103    }
     104
    95105    if (!hdu->header) {
    96         psTrace("psModules.camera", 5, "Reading the header...\n");
    97         hdu->header = psFitsReadHeader(hdu->header, fits);
    98         if (! hdu->header) {
    99             psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname);
    100             return false;
    101         }
    102     }
     106        hdu->header = header;
     107        return true;
     108    }
     109
     110    psMetadataIterator *iter = psMetadataIteratorAlloc(header, PS_LIST_HEAD, NULL); // Iterator
     111    psMetadataItem *item;           // Item from iteration
     112    while ((item = psMetadataGetAndIncrement(iter))) {
     113        const char *name = item->name; // Name of item
     114        if (psMetadataLookup(hdu->header, name)) {
     115            // It exists; clobber
     116            psMetadataRemoveKey(hdu->header, name);
     117        }
     118        psMetadataAddItem(hdu->header, item, PS_LIST_TAIL, 0);
     119    }
     120    psFree(iter);
     121    psFree(header);
    103122
    104123    return true;

  • branches/eam_branches/20090522/psModules/src/concepts/pmConceptsStandard.c

    r23819 r24557  
    3636// Format type for time
    3737typedef enum {
    38   TIME_FORMAT_YYYYMMDD,                 // Date stored in YYYY-MM-DD order (ISO-standard)
    39   TIME_FORMAT_DDMMYYYY,                 // Date stored in DD-MM-YYYY order
    40   TIME_FORMAT_MMDDYYYY,                 // Date stored in MM-DD-YYYY order
    41   TIME_FORMAT_JD,                       // Date stored as JD
    42   TIME_FORMAT_MJD,                      // Date stored as MJD
     38  TIME_FORMAT_YYYYMMDD,                 // Date stored in YYYY-MM-DD order (ISO-standard)
     39  TIME_FORMAT_DDMMYYYY,                 // Date stored in DD-MM-YYYY order
     40  TIME_FORMAT_MMDDYYYY,                 // Date stored in MM-DD-YYYY order
     41  TIME_FORMAT_JD,                       // Date stored as JD
     42  TIME_FORMAT_MJD,                      // Date stored as MJD
    4343} conceptTimeFormatType;
    4444
    … …  
    391391}
    392392
     393
    393394// FPA.RA and FPA.DEC
    394395psMetadataItem *p_pmConceptFormat_FPA_Coords(const psMetadataItem *concept,
    … …  
    410411    // How to interpret the coordinates
    411412    bool mdok = true;                   // Status of MD lookup
    412     psMetadata *formats = psMetadataLookupMetadata(&mdok,
    413                                                    cameraFormat,
    414                                                    "FORMATS");
     413    psMetadata *formats = psMetadataLookupMetadata(&mdok, cameraFormat, "FORMATS");
     414    bool sexagesimal = false;           // Write sexagesimal format?
    415415    if (mdok && formats) {
    416         psString format = psMetadataLookupStr(&mdok,formats, concept->name);
     416        psString format = psMetadataLookupStr(&mdok, formats, concept->name);
    417417        if (mdok && strlen(format) > 0) {
    418418            if (strcasecmp(format, "HOURS") == 0) {
    … …  
    428428            coords /= defaultCoordScaling(concept);
    429429        }
     430
     431        psString ra = psMetadataLookupStr(&mdok, formats, "FPA.RA"); // Format for RA
     432        psString dec = psMetadataLookupStr(&mdok, formats, "FPA.DEC"); // Format for Dec
     433        if (ra && strcasecmp(ra, "HOURS") == 0 && dec && strcasecmp(dec, "DEGREES") == 0) {
     434            sexagesimal = true;
     435        }
    430436    } else {
    431437        coords /= defaultCoordScaling(concept);
    432438    }
    433439
    434     // We choose to write sexagesimal format
    435     int big, medium;                    // Degrees and minutes
    436     float small;                        // Seconds
    437     bool negative = (coords < 0);       // Are we working below zero?
    438     coords = fabs(coords);
    439     big = (int)abs(coords);
    440     coords -= big;
    441     medium = 60.0 * coords;
    442     coords -= medium / 60.0;
    443     small = 3600.0 * coords;
    444     small = (float)((int)(small * 1000.0)) / 1000.0;
    445     if (negative) {
    446         big *= -1;
    447     }
    448     psString coordString = NULL;        // String with the coordinates in sexagesimal format
    449     psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small);
    450     psMetadataItem *coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString);
    451     psFree(coordString);
     440    psMetadataItem *coordItem = NULL;   // Item with coordinates, to return
     441    if (sexagesimal) {
     442        int big, medium;                    // Degrees and minutes
     443        float small;                        // Seconds
     444        bool negative = (coords < 0);       // Are we working below zero?
     445        coords = fabs(coords);
     446        big = (int)abs(coords);
     447        coords -= big;
     448        medium = 60.0 * coords;
     449        coords -= medium / 60.0;
     450        small = 3600.0 * coords;
     451        small = (float)((int)(small * 1000.0)) / 1000.0;
     452        if (negative) {
     453            big *= -1;
     454        }
     455        psString coordString = NULL;        // String with the coordinates in sexagesimal format
     456        psStringAppend(&coordString, "%d:%02d:%06.3f", big, medium, small);
     457        coordItem = psMetadataItemAllocStr(concept->name, concept->comment, coordString);
     458        psFree(coordString);
     459    } else {
     460        coordItem = psMetadataItemAllocF64(concept->name, concept->comment, coords);
     461    }
    452462
    453463    return coordItem;

  • branches/eam_branches/20090522/psModules/src/config/pmConfig.c

    r23748 r24557  
    243243    unsigned int numBadLines = 0;
    244244    struct stat filestat;
     245
     246    pmErrorRegister();
    245247
    246248    psTrace("psModules.config", 3, "Loading %s configuration from file %s\n",
    … …  
    14961498                }
    14971499            }
     1500            globfree(&globList);
    14981501        }
    14991502        psFree (words);

  • branches/eam_branches/20090522/psModules/src/config/pmErrorCodes.dat

    r23688 r24557  
    1212SKY                     Problem in sky
    1313STAMPS                  Unable to select stamps for PSF-matching
     14SMALL_AREA              Small area for PSF-matching
    1415# these errors correspond to standard exit conditions
    1516ARGUMENTS               Incorrect arguments

  • branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.c

    r23600 r24557  
    105105        DETREND_STRING_CASE(FRINGE);
    106106        DETREND_STRING_CASE(ASTROM);
     107        DETREND_STRING_CASE(NOISEMAP);
    107108    default:
    108109        return NULL;

  • branches/eam_branches/20090522/psModules/src/detrend/pmDetrendDB.h

    r23487 r24557  
    3535    PM_DETREND_TYPE_BACKGROUND,
    3636    PM_DETREND_TYPE_ASTROM,
     37    PM_DETREND_TYPE_NOISEMAP,
    3738} pmDetrendType;
    3839

  • branches/eam_branches/20090522/psModules/src/detrend/pmDetrendThreads.c

    r21457 r24557  
    6060    }
    6161
     62    // NOISEMAP : for now, not applied in the threaded loop
     63
    6264    return true;
    6365}

  • branches/eam_branches/20090522/psModules/src/imcombine/pmSubtraction.c

    r23989 r24557  
    196196{
    197197    psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
    198     int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
    199198    for (int y = -footprint; y <= footprint; y++) {
    200         // Convolution with a delta function is just the value specified by the offset
    201         memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes);
     199        for (int x = -footprint; x <= footprint; x++) {
     200            convolved->kernel[y][x] = image->kernel[y - v][x - u];
     201        }
    202202    }
    203203    return convolved;
    204204}
    205205
    206 // Take a subset of an image
    207 static inline psImage *convolveSubsetAlloc(psImage *image, // Image to be convolved
    208                                            psRegion region, // Region of interest (with border)
    209                                            bool threaded // Are we running threaded?
    210                                            )
    211 {
    212     if (!image) {
    213         return NULL;
    214     }
    215     // XXX if (threaded) {
    216     // XXX     psMutexLock(image);
    217     // XXX }
    218     psImage *subset = psImageSubset(image, region); // Subset image, to return
    219     // XXX if (threaded) {
    220     // XXX     psMutexUnlock(image);
    221     // XXX }
    222     return subset;
    223 }
    224 
    225 // Free the subset of an image
    226 static inline void convolveSubsetFree(psImage *parent, // Parent image
    227                                       psImage *child, // Child (subset) image
    228                                       bool threaded // Are we running threaded?
    229                                       )
    230 {
    231     if (!child || !parent) {
    232         return;
    233     }
    234     // XXX if (threaded) {
    235     // XXX     psMutexLock(parent);
    236     // XXX }
    237     psFree(child);
    238     // XXX if (threaded) {
    239     // XXX     psMutexUnlock(parent);
    240     // XXX }
    241     return;
    242 }
    243206
    244207// Convolve an image using FFT
    … …  
    256219                                  region.y0 - size, region.y1 + size); // Add a border
    257220
    258     bool threaded = pmSubtractionThreaded(); // Are we running threaded?
    259 
    260     psImage *subImage = convolveSubsetAlloc(image, border, threaded); // Subimage to convolve
    261     psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Subimage mask
     221    psImage *subImage = image ? psImageSubset(image, border) : NULL; // Subimage to convolve
     222    psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask
    262223
    263224    psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution
    264225
    265     convolveSubsetFree(image, subImage, threaded);
    266     convolveSubsetFree(mask, subMask, threaded);
     226    psFree(subImage);
     227    psFree(subMask);
    267228
    268229    // Now, we have to stick it in where it belongs
    … …  
    300261                                  region.y0 - size, region.y1 + size); // Add a border
    301262
    302     bool threaded = pmSubtractionThreaded(); // Are we running threaded?
    303 
    304     psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map
    305     psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image
    306     psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask
     263    psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map
     264    psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image
     265    psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask
    307266
    308267    // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
    … …  
    310269    psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
    311270
    312     convolveSubsetFree(variance, subVariance, threaded);
    313     convolveSubsetFree(sys, subSys, threaded);
    314     convolveSubsetFree(mask, subMask, threaded);
     271    psFree(subVariance);
     272    psFree(subSys);
     273    psFree(subMask);
    315274
    316275    // Now, we have to stick it in where it belongs
    … …  
    420379        if (box > 0) {
    421380            int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
    422             bool threaded = pmSubtractionThreaded(); // Are we running threaded?
    423 
    424381            psRegion region = psRegionSet(colMin - box, colMax + box,
    425382                                          rowMin - box, rowMax + box); // Region to convolve
    426383
    427             psImage *image = convolveSubsetAlloc(subMask, region, threaded); // Mask to convolve
     384            psImage *image = subMask ? psImageSubset(subMask, region) : NULL; // Mask to convolve
    428385
    429386            psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad,
    430387                                                     -box, box, -box, box); // Convolved subtraction mask
    431388
    432             convolveSubsetFree(subMask, image, threaded);
     389            psFree(image);
    433390
    434391            psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns");
    … …  
    653610                  float value = 0.0;    // Value of convolved pixel
    654611                  int uMin = x - size, uMax = x + size; // Range for u
    655                   psF32 *xKernelData = xKernel->data.F32; // Kernel values
     612                  psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values
    656613                  psF32 *imageData = &image->kernel[y][uMin]; // Image values
    657                   for (int u = uMin; u <= uMax; u++, xKernelData++, imageData++) {
     614                  for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
    658615                      value += *xKernelData * *imageData;
    659616                  }
    … …  
    668625                  float value = 0.0;    // Value of convolved pixel
    669626                  int vMin = y - size, vMax = y + size; // Range for v
    670                   psF32 *yKernelData = yKernel->data.F32; // Kernel values
     627                  psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values
    671628                  psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
    672                   for (int v = vMin; v <= vMax; v++, yKernelData++, imageData++) {
     629                  for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
    673630                      value += *yKernelData * *imageData;
    674631                  }

  • branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionEquation.c

    r24093 r24557  
    1515#include "pmSubtractionVisual.h"
    1616
    17 //#define TESTING
     17//#define TESTING                         // TESTING output for debugging; may not work with threads!
     18
     19#define USE_VARIANCE                    // Include variance in equation?
    1820
    1921//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    3133    for (int y = - footprint; y <= footprint; y++) {
    3234        for (int x = - footprint; x <= footprint; x++) {
    33             sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x];
     35            double value = image1->kernel[y][x] * image2->kernel[y][x];
     36#ifdef USE_VARIANCE
     37            value /= variance->kernel[y][x];
     38#endif
     39            sum += value;
    3440        }
    3541    }
    … …  
    195201            for (int y = - footprint; y <= footprint; y++) {
    196202                for (int x = - footprint; x <= footprint; x++) {
    197                     sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x];
     203                    double value = conv->kernel[y][x];
     204#ifdef USE_VARIANCE
     205                    value /= variance->kernel[y][x];
     206#endif
     207                    sumC += value;
    198208                }
    199209            }
    … …  
    218228        for (int y = - footprint; y <= footprint; y++) {
    219229            for (int x = - footprint; x <= footprint; x++) {
    220                 double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x];
     230                double invNoise2 = 1.0;
     231#ifdef USE_VARIANCE
     232                invNoise2 /= variance->kernel[y][x];
     233#endif
    221234                double value = input->kernel[y][x] * invNoise2;
    222235                sumI += value;
    … …  
    277290        for (int y = - footprint; y <= footprint; y++) {
    278291            for (int x = - footprint; x <= footprint; x++) {
    279                 sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x];
     292                double value = target->kernel[y][x] * conv->kernel[y][x];
     293#ifdef USE_VARIANCE
     294                value /= variance->kernel[y][x];
     295#endif
     296                sumTC += value;
    280297            }
    281298        }
    … …  
    297314        for (int y = - footprint; y <= footprint; y++) {
    298315            for (int x = - footprint; x <= footprint; x++) {
    299                 float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x];
     316                double value = target->kernel[y][x];
     317#ifdef USE_VARIANCE
     318                value /= variance->kernel[y][x];
     319#endif
    300320                sumIT += value * input->kernel[y][x];
    301321                sumT += value;
    … …  
    366386        for (int y = - footprint; y <= footprint; y++) {
    367387            for (int x = - footprint; x <= footprint; x++) {
    368                 sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x];
     388                double value = conv->kernel[y][x];
     389#ifdef USE_VARIANCE
     390                value /= variance->kernel[y][x];
     391#endif
     392                sumC += value;
    369393            }
    370394        }

  • branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionKernels.c

    r20421 r24557  
    145145
    146146                // Calculate moments
    147                 double sum = 0.0;       // Sum of kernel component, for normalisation
    148147                double moment = 0.0;    // Moment, for penalty
    149148                for (int v = -size, y = 0; v <= size; v++, y++) {
    150149                    for (int u = -size, x = 0; u <= size; u++, x++) {
    151150                        double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
    152                         sum += value;
    153151                        moment += value * (PS_SQR(u) + PS_SQR(v));
    154152                    }
    … …  
    163161                        }
    164162                    }
    165                     sum = 1.0 / sum;
     163                    sum = 1.0 / sqrt(sum);
    166164                    psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    167165                    psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    168                     moment *= sum;
     166                    moment *= PS_SQR(sum);
    169167                }
    170168

  • branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMask.c

    r23851 r24557  
    66#include <pslib.h>
    77
     8#include "pmErrorCodes.h"
    89#include "pmSubtraction.h"
    910#include "pmSubtractionKernels.h"
    … …  
    7879        }
    7980        if (numBad > badFrac * numCols * numRows) {
    80             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     81            psError(PM_ERR_SMALL_AREA, true,
    8182                    "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n",
    8283                    numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac);

  • branches/eam_branches/20090522/psModules/src/imcombine/pmSubtractionMatch.c

    r24181 r24557  
    246246                                         badFrac, mode); // Subtraction mask
    247247    if (!subMask) {
    248         psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
     248        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    249249        psFree(kernels);
    250250        psFree(regions);
    … …  
    337337        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false);
    338338    }
    339     PS_ASSERT_INT_POSITIVE(iter, false);
     339    PS_ASSERT_INT_NONNEGATIVE(iter, false);
    340340    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
    341341
    … …  
    379379                                badFrac, subMode);
    380380    if (!subMask) {
    381         psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
     381        psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
    382382        goto MATCH_ERROR;
    383383    }

  • branches/eam_branches/20090522/psModules/src/objects/Makefile.am

    r23990 r24557  
    3838        pmSourceIO_PS1_CAL_0.c \
    3939        pmSourceIO_CMF_PS1_V1.c \
     40        pmSourceIO_CMF_PS1_V2.c \
    4041        pmSourceIO_MatchedRefs.c \
    4142        pmSourcePlots.c \

  • branches/eam_branches/20090522/psModules/src/objects/pmFootprintArraysMerge.c

    r20937 r24557  
    2828 */
    2929psArray *pmFootprintArraysMerge(const psArray *footprints1, // one set of footprints
    30                                 const psArray *footprints2, // the other set
    31                                 const int includePeaks) { // which peaks to set? 0x1 => footprints1, 0x2 => 2
    32     assert (footprints1->n == 0 || pmFootprintTest(footprints1->data[0]));
    33     assert (footprints2->n == 0 || pmFootprintTest(footprints2->data[0]));
     30                                const psArray *footprints2, // the other set
     31                                const int includePeaks // which peaks to set? 0x1 => footprints1, 0x2 => 2
     32    )
     33{
     34    if (!footprints1 && !footprints2) {
     35        // No footprints in merged array
     36        return psArrayAllocEmpty(0);
     37    }
    3438
    35     if (footprints1->n == 0 || footprints2->n == 0) {           // nothing to do but put copies on merged
    36         const psArray *old = (footprints1->n == 0) ? footprints2 : footprints1;
     39    assert(!footprints1 || footprints1->n == 0 || pmFootprintTest(footprints1->data[0]));
     40    assert(!footprints2 || footprints2->n == 0 || pmFootprintTest(footprints2->data[0]));
    3741
    38         psArray *merged = psArrayAllocEmpty(old->n);
    39         for (int i = 0; i < old->n; i++) {
    40             psArrayAdd(merged, 1, old->data[i]);
    41         }
    42        
    43         return merged;
     42    if (!footprints1 || footprints1->n == 0 || !footprints2 || footprints2->n == 0) {
     43        // nothing to do but put copies on merged
     44        const psArray *old = (!footprints1 || footprints1->n == 0) ? footprints2 : footprints1;
     45
     46        psArray *merged = psArrayAllocEmpty(old->n);
     47        for (int i = 0; i < old->n; i++) {
     48            psArrayAdd(merged, 1, old->data[i]);
     49        }
     50
     51        return merged;
    4452    }
    4553    /*
    … …  
    4856     */
    4957    {
    50         pmFootprint *fp1 = footprints1->data[0];
    51         pmFootprint *fp2 = footprints2->data[0];
    52         if (fp1->region.x0 != fp2->region.x0 ||
    53             fp1->region.x1 != fp2->region.x1 ||
    54             fp1->region.y0 != fp2->region.y0 ||
    55             fp1->region.y1 != fp2->region.y1) {
    56             psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    57                     "The two pmFootprint arrays correspnond to different-sized regions");
    58             return NULL;
    59         }
     58        pmFootprint *fp1 = footprints1->data[0];
     59        pmFootprint *fp2 = footprints2->data[0];
     60        if (fp1->region.x0 != fp2->region.x0 ||
     61            fp1->region.x1 != fp2->region.x1 ||
     62            fp1->region.y0 != fp2->region.y0 ||
     63            fp1->region.y1 != fp2->region.y1) {
     64            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     65                    "The two pmFootprint arrays correspnond to different-sized regions");
     66            return NULL;
     67        }
    6068    }
    6169    /*
    … …  
    7280     * Now assign the peaks appropriately.  We could do this more efficiently
    7381     * using idImage (which we just freed), but this is easy and probably fast enough
    74      */ 
     82     */
    7583    if (includePeaks & 0x1) {
    76         psArray *peaks = pmFootprintArrayToPeaks(footprints1);
    77         pmFootprintsAssignPeaks(merged, peaks);
    78         psFree(peaks);
     84        psArray *peaks = pmFootprintArrayToPeaks(footprints1);
     85        pmFootprintsAssignPeaks(merged, peaks);
     86        psFree(peaks);
    7987    }
    8088
    8189    if (includePeaks & 0x2) {
    82         psArray *peaks = pmFootprintArrayToPeaks(footprints2);
    83         pmFootprintsAssignPeaks(merged, peaks);
    84         psFree(peaks);
     90        psArray *peaks = pmFootprintArrayToPeaks(footprints2);
     91        pmFootprintsAssignPeaks(merged, peaks);
     92        psFree(peaks);
    8593    }
    86    
     94
    8795    return merged;
    8896}

  • branches/eam_branches/20090522/psModules/src/objects/pmMoments.c

    r23487 r24557  
    5151    tmp->Sky = 0.0;
    5252    tmp->nPixels = 0;
     53    tmp->SN = 0;
    5354
    5455    psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);

  • branches/eam_branches/20090522/psModules/src/objects/pmSource.c

    r24529 r24557  
    480480
    481481        if (!psVectorStats (stats, tmpSx, NULL, NULL, 0)) {
    482             psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats");
     482            psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats");
    483483            return (emptyClump);
    484         }
     484        }
    485485        psfClump.X  = stats->clippedMean;
    486486        psfClump.dX = stats->clippedStdev;
    487487
    488488        if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) {
    489             psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats");
     489            psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats");
    490490            return (emptyClump);
    491         }
     491        }
    492492        psfClump.Y  = stats->clippedMean;
    493493        psfClump.dY = stats->clippedStdev;
    … …  
    548548        pmSource *source = (pmSource *) sources->data[i];
    549549
    550         // psf clumps are found for image subregions:
    551         // skip sources not in this region
     550        // psf clumps are found for image subregions:
     551        // skip sources not in this region
    552552        if (source->peak->x <  region->x0) continue;
    553553        if (source->peak->x >= region->x1) continue;
    554554        if (source->peak->y <  region->y0) continue;
    555         if (source->peak->y >= region->y1) continue;
    556 
    557         // should be set by pmSourceAlloc
     555        if (source->peak->y >= region->y1) continue;
     556
     557        // should be set by pmSourceAlloc
    558558        psAssert (source->type == PM_SOURCE_TYPE_UNKNOWN, "source type was not init-ed?");
    559559
    … …  
    561561        if (!source->moments) {
    562562            source->type = PM_SOURCE_TYPE_STAR;
    563             psAssert (source->mode & noMoments, "why is this source missing moments?");
     563            psAssert (source->mode & noMoments, "why is this source missing moments?");
    564564            Nstar++;
    565565            continue;
    … …  
    596596        }
    597597
    598         // likely defect (too small to be stellar) (push out to 3 sigma)
    599         // low S/N objects which are small are probably stellar
    600         // only set candidate defects if
    601         // XXX these limits are quite arbitrary
    602         if ((sigX < 0.05) || (sigY < 0.05)) {
    603             source->type = PM_SOURCE_TYPE_DEFECT;
    604             source->mode |= PM_SOURCE_MODE_DEFECT;
    605             Ncr ++;
    606             continue;
    607         }
    608 
    609         // likely unsaturated extended source (too large to be stellar)
    610         if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    611             source->type = PM_SOURCE_TYPE_EXTENDED;
    612             Next ++;
    613             continue;
    614         }
    615 
    616         // the rest are probable stellar objects
    617         starsn_moments->data.F32[starsn_moments->n] = source->moments->SN;
    618         starsn_moments->n ++;
    619         starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN;
    620         starsn_peaks->n ++;
    621         Nstar ++;
    622 
    623         // PSF star (within 1.5 sigma of clump center, S/N > limit)
    624         psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    625         if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
    626             source->type = PM_SOURCE_TYPE_STAR;
    627             source->mode |= PM_SOURCE_MODE_PSFSTAR;
    628             Npsf ++;
    629             continue;
     598        // The following determinations require the use of moments
     599        if (!(source->mode & noMoments)) {
     600            // likely defect (too small to be stellar) (push out to 3 sigma)
     601            // low S/N objects which are small are probably stellar
     602            // XXX these limits are quite arbitrary
     603            if (sigX < 0.05 || sigY < 0.05) {
     604                source->type = PM_SOURCE_TYPE_DEFECT;
     605                source->mode |= PM_SOURCE_MODE_DEFECT;
     606                Ncr ++;
     607                continue;
     608            }
     609
     610            // likely unsaturated extended source (too large to be stellar)
     611            if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) {
     612                source->type = PM_SOURCE_TYPE_EXTENDED;
     613                Next ++;
     614                continue;
     615            }
     616
     617            // the rest are probable stellar objects
     618            starsn_moments->data.F32[starsn_moments->n] = source->moments->SN;
     619            starsn_moments->n ++;
     620            starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN;
     621            starsn_peaks->n ++;
     622            Nstar ++;
     623
     624            // PSF star (within 1.5 sigma of clump center, S/N > limit)
     625            psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
     626            if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
     627                source->type = PM_SOURCE_TYPE_STAR;
     628                source->mode |= PM_SOURCE_MODE_PSFSTAR;
     629                Npsf ++;
     630                continue;
     631            }
    630632        }
    631633
    … …  
    642644
    643645        if (!psVectorStats (stats, starsn_moments, NULL, NULL, 0)) {
    644             psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
    645             psFree (stats);
    646             psFree (starsn_peaks);
    647             return false;
     646            psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
     647            psFree (stats);
     648            psFree (starsn_peaks);
     649            return false;
    648650        }
    649651        psLogMsg ("pmObjects", 3, "SN range (moments): %f - %f\n", stats->min, stats->max);
    … …  
    656658        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    657659        if (!psVectorStats (stats, starsn_peaks, NULL, NULL, 0)) {
    658             psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
    659             psFree (stats);
    660             psFree (starsn_peaks);
    661             return false;
     660            psError(PS_ERR_UNKNOWN, false, "failed to measure SN / moments stats");
     661            psFree (stats);
     662            psFree (starsn_peaks);
     663            return false;
    662664        }
    663665        psLogMsg ("psModules.objects", 3, "SN range (peaks)  : %f - %f (%ld)\n",

  • branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.c

    r23990 r24557  
    496496                status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
    497497            }
     498            if (!strcmp (exttype, "PS1_V2")) {
     499                status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
     500            }
    498501            if (xsrcname) {
    499502              if (!strcmp (exttype, "PS1_DEV_1")) {
    … …  
    506509                  status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
    507510              }
     511              if (!strcmp (exttype, "PS1_V2")) {
     512                  status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
     513              }
    508514            }
    509515            if (xfitname) {
    … …  
    517523                  status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
    518524              }
     525              if (!strcmp (exttype, "PS1_V2")) {
     526                  status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
     527              }
    519528            }
    520529            if (!status) {
    … …  
    944953                sources = pmSourcesRead_CMF_PS1_V1 (file->fits, hdu->header);
    945954            }
     955            if (!strcmp (exttype, "PS1_V2")) {
     956                sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
     957            }
    946958        }
    947959

  • branches/eam_branches/20090522/psModules/src/objects/pmSourceIO.h

    r23990 r24557  
    3939bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname);
    4040
     41bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
     42bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
     43bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname);
     44
    4145bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    4246
    … …  
    4852psArray *pmSourcesRead_PS1_CAL_0 (psFits *fits, psMetadata *header);
    4953psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
     54psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
    5055
    5156bool pmSourcesWritePSFs (psArray *sources, char *filename);

  • branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_MatchedRefs.c

    r24025 r24557  
    4949    pmCell *cell = NULL;
    5050    pmReadout *readout = NULL;
    51     pmFPAview *view = pmFPAviewAlloc (0);
    5251
    5352    // first, check if there are any matches to be written
    … …  
    5756    psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSASTRO");
    5857    if (!status) {
    59         psError(PS_ERR_UNKNOWN, true, "missing recipe PSASTRO in config data");
    60         return false;
     58        psError(PS_ERR_UNKNOWN, true, "missing recipe PSASTRO in config data");
     59        return false;
    6160    }
    6261
    … …  
    6564
    6665    psArray *table = psArrayAllocEmpty (0x1000);
     66    pmFPAview *view = pmFPAviewAlloc (0);
    6767
    6868    // this loop selects the matched stars for all chips
    … …  
    7070        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    7171        if (!chip->process || !chip->file_exists) continue;
    72        
    73         char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
    7472
    75         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
     73        char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
     74
     75        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
    7676            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    7777            if (!cell->process || !cell->file_exists) continue;
    7878
    79             // process each of the readouts
    80             // XXX there can only be one readout per chip, right?
    81             while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
    82                 if (! readout->data_exists) continue;
     79            // process each of the readouts
     80            // XXX there can only be one readout per chip, right?
     81            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
     82                if (! readout->data_exists) continue;
    8383
    84                 // select the raw objects for this readout
    85                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    86                 if (rawstars == NULL) continue;
     84                // select the raw objects for this readout
     85                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     86                if (rawstars == NULL) continue;
    8787
    88                 // select the raw objects for this readout
    89                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    90                 if (refstars == NULL) continue;
    91                 psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
     88                // select the raw objects for this readout
     89                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     90                if (refstars == NULL) continue;
     91                psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
    9292
    93                 psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
    94                 if (matches == NULL) continue;
     93                psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
     94                if (matches == NULL) continue;
    9595
    96                 for (int i = 0; i < matches->n; i++) {
    97                     pmAstromMatch *match = matches->data[i];
     96                for (int i = 0; i < matches->n; i++) {
     97                    pmAstromMatch *match = matches->data[i];
    9898
    99                     pmAstromObj *raw = rawstars->data[match->raw];
    100                     pmAstromObj *ref = refstars->data[match->ref];
     99                    pmAstromObj *raw = rawstars->data[match->raw];
     100                    pmAstromObj *ref = refstars->data[match->ref];
    101101
    102                     psMetadata *row = psMetadataAlloc ();
    103                     psMetadataAdd (row, PS_LIST_TAIL, "RA",       PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r);
    104                     psMetadataAdd (row, PS_LIST_TAIL, "DEC",      PS_DATA_F64, "declination (deg, J2000)",     PM_DEG_RAD*ref->sky->d);
    105                     psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP",   PS_DATA_F32, "x coord on chip",              raw->chip->x);
    106                     psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP",   PS_DATA_F32, "y coord on chip",              raw->chip->y);   
    107                     psMetadataAdd (row, PS_LIST_TAIL, "X_FPA",    PS_DATA_F32, "x coord on focal plane",       raw->FP->x);
    108                     psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA",    PS_DATA_F32, "y coord on focal plane",       raw->FP->y);
    109                     psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude",       raw->Mag);
    110                     psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF",  PS_DATA_F32, "reference star magnitude",     ref->Mag);
    111                     psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF",PS_DATA_F32, "reference star color",         ref->Color);
    112                     psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID",  PS_DATA_STRING, "chip identifier",           chipName);
    113                     // XXX need to add the reference color, but this needs getstar / dvo.photcodes for the reference to be refined.
     102                    psMetadata *row = psMetadataAlloc ();
     103                    psMetadataAdd (row, PS_LIST_TAIL, "RA",       PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r);
     104                    psMetadataAdd (row, PS_LIST_TAIL, "DEC",      PS_DATA_F64, "declination (deg, J2000)",     PM_DEG_RAD*ref->sky->d);
     105                    psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP",   PS_DATA_F32, "x coord on chip",              raw->chip->x);
     106                    psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP",   PS_DATA_F32, "y coord on chip",              raw->chip->y);
     107                    psMetadataAdd (row, PS_LIST_TAIL, "X_FPA",    PS_DATA_F32, "x coord on focal plane",       raw->FP->x);
     108                    psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA",    PS_DATA_F32, "y coord on focal plane",       raw->FP->y);
     109                    psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude",       raw->Mag);
     110                    psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF",  PS_DATA_F32, "reference star magnitude",     ref->Mag);
     111                    psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF",PS_DATA_F32, "reference star color",         ref->Color);
     112                    psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID",  PS_DATA_STRING, "chip identifier",           chipName);
     113                    // XXX need to add the reference color, but this needs getstar / dvo.photcodes for the reference to be refined.
    114114
    115                     psArrayAdd (table, 100, row);
    116                     psFree (row);
    117                 }
    118             }
    119         }
     115                    psArrayAdd (table, 100, row);
     116                    psFree (row);
     117                }
     118            }
     119        }
    120120    }
    121121    psFree (view);
    122122
    123123    if (table->n == 0) {
    124         psFree(table);
    125         return true;
     124        psFree(table);
     125        return true;
    126126    }
    127127

  • branches/eam_branches/20090522/psModules/src/objects/pmSourceMatch.c

    r23989 r24557  
    4545                         psVector **x, psVector **y, // Coordinate vectors to return
    4646                         psVector **mag, psVector **magErr, // Magnitude and error vectors to return
     47                         psVector **indices, // Indices for sources
    4748                         const psArray *sources // Input sources
    4849    )
    … …  
    5859    *mag = psVectorRecycle(*mag, numSources, PS_TYPE_F32);
    5960    *magErr = psVectorRecycle(*magErr, numSources, PS_TYPE_F32);
     61    *indices = psVectorRecycle(*indices, numSources, PS_TYPE_S32);
    6062    float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources
    6163    long num = 0;                       // Number of valid sources
    6264    for (long i = 0; i < numSources; i++) {
    6365        pmSource *source = sources->data[i]; // Source of interest
    64         if (!source) continue;
    65         if (source->mode & SOURCE_MASK) continue;
    66         if (!isfinite(source->psfMag)) continue;
    67         if (!isfinite(source->errMag)) continue;
    68         if (source->psfMag > SOURCE_FAINTEST) continue;
     66        if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) ||
     67            !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) {
     68            continue;
     69        }
    6970
    7071        float xSrc, ySrc;               // Coordinates of source
    … …  
    7879        (*y)->data.F32[num] = ySrc;
    7980        (*mag)->data.F32[num] = source->psfMag;
    80         (*magErr)->data.F32[num] = source->errMag ;
     81        (*magErr)->data.F32[num] = source->errMag;
     82        (*indices)->data.S32[num] = i;
    8183        num++;
    8284    }
    … …  
    8587    (*mag)->n = num;
    8688    (*magErr)->n = num;
     89    (*indices)->n = num;
    8790
    8891    if (*bounds) {
    … …  
    175178        psVector *xImage = NULL, *yImage = NULL; // Coordinates of sources
    176179        psVector *magImage = NULL, *magErrImage = NULL; // Magnitude and mag
    177 
    178         int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage,
     180        psVector *indices = NULL;     // Indices for sources
     181
     182        int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage, &indices,
    179183                                      sources); // Number of sources
    180184
    … …  
    187191            for (int j = 0; j < numSources; j++) {
    188192                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    189                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     193                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     194                                 i, indices->data.S32[j]);
    190195                matches->data[j] = match;
    191196            }
    … …  
    206211            for (int j = 0, k = numMaster; j < numSources; j++, k++) {
    207212                pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    208                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     213                pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     214                                 i, indices->data.S32[j]);
    209215                matches->data[k] = match;
    210216            }
    … …  
    231237                    // Record the match
    232238                    pmSourceMatch *match = matches->data[index]; // Match data
    233                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     239                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     240                                     i, indices->data.S32[j]);
    234241                    numMatch++;
    235242                } else {
    236243                    // Add to the master list
    237244                    pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
    238                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j);
     245                    pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
     246                                     i, indices->data.S32[j]);
    239247                    xMaster->data.F32[numMaster] = xImage->data.F32[j];
    240248                    yMaster->data.F32[numMaster] = yImage->data.F32[j];
Note: See TracChangeset for help on using the changeset viewer.