IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 7, 2005, 4:04:22 PM (21 years ago)
Author:
Paul Price
Message:

Using autoconf; some bug fixes to the polynomial order

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/stacScales.c

    r5743 r5745  
    66#define SQUARE(x) ((x)*(x))
    77
    8 #define SAMPLE 10                       // Subsample rate for images
     8#define SAMPLE 10                       // Subsample rate for images
    99
    1010float stacBackground(const psImage *image, // Image for which to get the background
    11                      int sample         // Sample in increments of this value
     11                     int sample         // Sample in increments of this value
    1212    )
    1313{
    … …  
    1515
    1616#if 1
    17     int size = image->numCols * image->numRows; // Number of pixels in image
    18     int numSamples = size / sample;     // Number of samples in image
     17    int size = image->numCols * image->numRows; // Number of pixels in image
     18    int numSamples = size / sample;     // Number of samples in image
    1919    psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample
    20     psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
    21 
    22     int offset = 0;                     // Offset from start of the row
    23     int index = 0;                      // Sample number
     20    psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
     21
     22    int offset = 0;                     // Offset from start of the row
     23    int index = 0;                      // Sample number
    2424    for (int row = 0; row < image->numRows; row++) {
    25         // I'd cast this as a "for", but this makes it a bit easier to understand.
    26         int col = offset;
    27         while (col < image->numCols) {
    28             values->data.F32[index] = image->data.F32[row][col];
    29             if (isnan(values->data.F32[index])) {
    30                 mask->data.U8[index] = 1;
    31             } else {
    32                 mask->data.U8[index] = 0;
    33             }
    34             col += sample;
    35             index++;
    36         }
    37         offset = col - image->numCols;
     25        // I'd cast this as a "for", but this makes it a bit easier to understand.
     26        int col = offset;
     27        while (col < image->numCols) {
     28            values->data.F32[index] = image->data.F32[row][col];
     29            if (isnan(values->data.F32[index])) {
     30                mask->data.U8[index] = 1;
     31            } else {
     32                mask->data.U8[index] = 0;
     33            }
     34            col += sample;
     35            index++;
     36        }
     37        offset = col - image->numCols;
    3838    }
    3939
    … …  
    5858
    5959
    60 bool stacScales(psVector **scalesPtr,   // Scales to return
    61                 psVector **offsetsPtr,  // Offsets to return
    62                 const psArray *images,  // Images on which to measure the scales and offsets
    63                 const char *starFile,   // File containing coordinates to photometer
    64                 const char *starMapFile, // Map for coodinates to the common output frame
    65                 float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
    66                 float aper              // Aperture to use for photometry (radius)
     60bool stacScales(psVector **scalesPtr,   // Scales to return
     61                psVector **offsetsPtr,  // Offsets to return
     62                const psArray *images,  // Images on which to measure the scales and offsets
     63                const char *starFile,   // File containing coordinates to photometer
     64                const char *starMapFile, // Map for coodinates to the common output frame
     65                float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
     66                float aper              // Aperture to use for photometry (radius)
    6767    )
    6868{
    … …  
    7171    assert(images);
    7272    for (int i = 0; i < images->n; i++) {
    73         psImage *image = images->data[i];
    74         assert(image->type.type == PS_TYPE_F32);
     73        psImage *image = images->data[i];
     74        assert(image->type.type == PS_TYPE_F32);
    7575    }
    7676
    … …  
    7878    psVector *offsets = NULL; // Offsets between images (ADU)
    7979    if (*scalesPtr) {
    80         scales = *scalesPtr;
    81         assert(scales);
    82         assert(scales->n == images->n);
    83         assert(scales->type.type == PS_TYPE_F32);
     80        scales = *scalesPtr;
     81        assert(scales);
     82        assert(scales->n == images->n);
     83        assert(scales->type.type == PS_TYPE_F32);
    8484    } else {
    85         scales = psVectorAlloc(images->n, PS_TYPE_F32);
    86         *scalesPtr = scales;
     85        scales = psVectorAlloc(images->n, PS_TYPE_F32);
     86        *scalesPtr = scales;
    8787    }
    8888    if (*offsetsPtr) {
    89         offsets = *offsetsPtr;
    90         assert(offsets);
    91         assert(offsets->n == images->n);
    92         assert(offsets->type.type == PS_TYPE_F32);
     89        offsets = *offsetsPtr;
     90        assert(offsets);
     91        assert(offsets->n == images->n);
     92        assert(offsets->type.type == PS_TYPE_F32);
    9393    } else {
    94         offsets = psVectorAlloc(images->n, PS_TYPE_F32);
    95         *offsetsPtr = offsets;
     94        offsets = psVectorAlloc(images->n, PS_TYPE_F32);
     95        *offsetsPtr = offsets;
    9696    }
    9797
    … …  
    100100    double startTime = getTime();
    101101    for (int i = 0; i < images->n; i++) {
    102         offsets->data.F32[i] = stacBackground(images->data[i], SAMPLE);
    103         psTrace("stac.scales", 5, "Background in image %d is %f\n", i, offsets->data.F32[i]);
    104         double time = getTime();
    105         psTrace("stac.scales", 10, "Took %f sec\n", time - startTime);
    106         startTime = time;
     102        offsets->data.F32[i] = stacBackground(images->data[i], SAMPLE);
     103        psTrace("stac.scales", 5, "Background in image %d is %f\n", i, offsets->data.F32[i]);
     104        double time = getTime();
     105        psTrace("stac.scales", 10, "Took %f sec\n", time - startTime);
     106        startTime = time;
    107107    }
    108108
    109109    // Now the scales
    110110    if (starFile == NULL || starMapFile == NULL) {
    111         psLogMsg("stac.scales", PS_LOG_INFO,
    112                 "No coordinates available to set scales --- assuming all are identical.\n");
    113         for (int i = 0; i < images->n; i++) {
    114             scales->data.F32[i] = 1.0;
    115             psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
    116         }
     111        psLogMsg("stac.scales", PS_LOG_INFO,
     112                "No coordinates available to set scales --- assuming all are identical.\n");
     113        for (int i = 0; i < images->n; i++) {
     114            scales->data.F32[i] = 1.0;
     115            psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
     116        }
    117117    } else {
    118         // Read star coordinates and map
    119         psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates
    120         psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates
    121 
    122         // Transform the stellar positions to match the transformed reference frame
    123         psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
    124         // Fix up difference between map and output frame
    125         starMap->x->coeff[0][0] -= xMapDiff;
    126         starMap->y->coeff[0][0] -= yMapDiff;
    127         for (int i = 0; i < starCoords->n; i++) {
    128             starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]);
    129         }
    130         psFree(starMap);
    131 
    132         psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
    133         psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
    134 
    135         // Set scales relative to the first image
    136         scales->data.F32[0] = 1.0;
    137 
    138         // Iterate over images
    139         for (int i = 0; i < scales->n; i++) {
    140             psImage *image = images->data[i]; // The image we're working with
    141 
    142             float background = offsets->data.F32[i]; // Background in image
    143 
    144             // Do photometry on transformed image i
    145             psTrace("stac.scales", 3, "Doing photometry on image %d....\n", i);
    146             psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
    147             psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
    148             for (int j = 0; j < starCoordsTransformed->n; j++) {
    149                 psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
    150                
    151                 if (coords->x < aper || coords->y < aper ||
    152                     coords->x + aper > image->numCols - 1 ||
    153                     coords->y + aper > image->numRows) {
    154                     mask->data.U8[j] = 1;
    155                 } else {
    156                     // Sum flux within the aperture
    157                     float sum = 0.0;
    158                     int numPix = 0;
    159                     float aper2 = SQUARE(aper);
    160                     for (int y = (int)floorf(coords->y - aper);
    161                          y <= (int)ceilf(coords->y + aper); y++) {
    162                         for (int x = (int)floorf(coords->x - aper);
    163                              x <= (int)ceilf(coords->x + aper); x++) {
    164                             if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <=
    165                                 aper2) {
    166                                 sum += image->data.F32[y][x];
    167                                 numPix++;
    168                             }
    169                         }
    170                     }
    171                     // Subtract background, renormalise to account for circular aperture
    172                     if (numPix > 0 && sum > 0) {
    173                         sum -= offsets->data.F32[i] * (float)numPix;
    174                         photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix;
    175                         if (photometry->data.F32[j] > 0 && finite(photometry->data.F32[j])) {
    176                             mask->data.U8[j] = 1;
    177                             psTrace("stac.scales", 8, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
    178                         } else {
    179                             mask->data.U8[j] = 0;
    180                         }
    181                     } else {
    182                         mask->data.U8[j] = 0;
    183                     }
    184                 }
    185             }
    186             stars->data[i] = photometry;
    187             masks->data[i] = mask;
    188         }
    189         psFree(starCoords);
    190         psFree(starCoordsTransformed);
    191 
    192         // Get the scales
    193         psVector *ref = stars->data[0]; // The reference photometry
    194         psVector *refMask = masks->data[0]; // The reference mask
     118        // Read star coordinates and map
     119        psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates
     120        psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates
     121
     122        // Transform the stellar positions to match the transformed reference frame
     123        psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
     124        // Fix up difference between map and output frame
     125        starMap->x->coeff[0][0] -= xMapDiff;
     126        starMap->y->coeff[0][0] -= yMapDiff;
     127        for (int i = 0; i < starCoords->n; i++) {
     128            starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]);
     129        }
     130        psFree(starMap);
     131
     132        psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
     133        psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
     134
     135        // Set scales relative to the first image
     136        scales->data.F32[0] = 1.0;
     137
     138        // Iterate over images
     139        for (int i = 0; i < scales->n; i++) {
     140            psImage *image = images->data[i]; // The image we're working with
     141
     142            // Do photometry on transformed image i
     143            psTrace("stac.scales", 3, "Doing photometry on image %d....\n", i);
     144            psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
     145            psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
     146            for (int j = 0; j < starCoordsTransformed->n; j++) {
     147                psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
     148
     149                if (coords->x < aper || coords->y < aper ||
     150                    coords->x + aper > image->numCols - 1 ||
     151                    coords->y + aper > image->numRows) {
     152                    mask->data.U8[j] = 1;
     153                } else {
     154                    // Sum flux within the aperture
     155                    float sum = 0.0;
     156                    int numPix = 0;
     157                    float aper2 = SQUARE(aper);
     158                    for (int y = (int)floorf(coords->y - aper);
     159                         y <= (int)ceilf(coords->y + aper); y++) {
     160                        for (int x = (int)floorf(coords->x - aper);
     161                             x <= (int)ceilf(coords->x + aper); x++) {
     162                            if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <=
     163                                aper2) {
     164                                sum += image->data.F32[y][x];
     165                                numPix++;
     166                            }
     167                        }
     168                    }
     169                    // Subtract background, renormalise to account for circular aperture
     170                    if (numPix > 0 && sum > 0) {
     171                        sum -= offsets->data.F32[i] * (float)numPix;
     172                        photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix;
     173                        if (photometry->data.F32[j] > 0 && finite(photometry->data.F32[j])) {
     174                            mask->data.U8[j] = 1;
     175                            psTrace("stac.scales", 8, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
     176                        } else {
     177                            mask->data.U8[j] = 0;
     178                        }
     179                    } else {
     180                        mask->data.U8[j] = 0;
     181                    }
     182                }
     183            }
     184            stars->data[i] = photometry;
     185            masks->data[i] = mask;
     186        }
     187        psFree(starCoords);
     188        psFree(starCoordsTransformed);
     189
     190        // Get the scales
     191        psVector *ref = stars->data[0]; // The reference photometry
     192        psVector *refMask = masks->data[0]; // The reference mask
    195193#if 0
    196         psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics
    197         stats->clipSigma = 2.5;
    198         stats->clipIter = 3;
     194        psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics
     195        stats->clipSigma = 2.5;
     196        stats->clipIter = 3;
    199197#else
    200         psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
    201 #endif
    202         for (int i = 1; i < scales->n; i++) {
    203             psVector *input = stars->data[i];   // The comparison photometry
    204             psVector *inputMask = masks->data[i]; // The comparison mask
    205 
    206             psVector *compare = (psVector*)psBinaryOp(NULL, input, "/", ref);
    207             psVector *compareMask = (psVector*)psBinaryOp(NULL, inputMask, "*", refMask);
    208             (void)psBinaryOp(compareMask, psScalarAlloc(1, PS_TYPE_U8), "-", compareMask);
     198        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
     199#endif
     200        for (int i = 1; i < scales->n; i++) {
     201            psVector *input = stars->data[i];   // The comparison photometry
     202            psVector *inputMask = masks->data[i]; // The comparison mask
     203
     204            psVector *compare = (psVector*)psBinaryOp(NULL, input, "/", ref);
     205            psVector *compareMask = (psVector*)psBinaryOp(NULL, inputMask, "*", refMask);
     206            (void)psBinaryOp(compareMask, psScalarAlloc(1, PS_TYPE_U8), "-", compareMask);
    209207
    210208#if 0
    211             psTrace("stac.scales", 9, "Getting scale for image %d...\n", i);
    212             for (int j = 0; j < compare->n; j++) {
    213                 if (compareMask->data.U8[j] & 1) {
    214                     psTrace("stac.scales", 9, "Bad star %d: %f %f --> %f %d\n", j, input->data.F32[j],
    215                             ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
    216                 } else {
    217                     psTrace("stac.scales", 9, "Good star %d: %f %f --> %f %d\n", j, input->data.F32[j],
    218                             ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
    219                 }
    220             }
    221 #endif
    222 
    223             psVectorStats(stats, compare, NULL, compareMask, 1);
     209            psTrace("stac.scales", 9, "Getting scale for image %d...\n", i);
     210            for (int j = 0; j < compare->n; j++) {
     211                if (compareMask->data.U8[j] & 1) {
     212                    psTrace("stac.scales", 9, "Bad star %d: %f %f --> %f %d\n", j, input->data.F32[j],
     213                            ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
     214                } else {
     215                    psTrace("stac.scales", 9, "Good star %d: %f %f --> %f %d\n", j, input->data.F32[j],
     216                            ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
     217                }
     218            }
     219#endif
     220
     221            psVectorStats(stats, compare, NULL, compareMask, 1);
    224222
    225223#if 0
    226             scales->data.F32[i] = stats->clippedMean;
     224            scales->data.F32[i] = stats->clippedMean;
    227225#else
    228             scales->data.F32[i] = stats->sampleMedian;
    229 #endif
    230             psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
    231             psFree(compare);
    232             psFree(compareMask);
    233         }
    234         psFree(stats);
    235 
    236         psFree(stars);
    237         psFree(masks);
     226            scales->data.F32[i] = stats->sampleMedian;
     227#endif
     228            psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
     229            psFree(compare);
     230            psFree(compareMask);
     231        }
     232        psFree(stats);
     233
     234        psFree(stars);
     235        psFree(masks);
    238236    }
    239237
    … …  
    242240
    243241
    244 bool stacRescale(psArray *images,       // Images to rescale
    245                  psArray *errImages,    // Variance images to rescale
    246                  const psImage *mask,   // Mask for pixels to scale
    247                 const psVector *scales,// Scales for images
    248                 const psVector *offsets // Offsets for images
     242bool stacRescale(psArray *images,       // Images to rescale
     243                 psArray *errImages,    // Variance images to rescale
     244                 const psImage *mask,   // Mask for pixels to scale
     245                const psVector *scales,// Scales for images
     246                const psVector *offsets // Offsets for images
    249247    )
    250248{
    … …  
    258256    assert(offsets->type.type == PS_TYPE_F32);
    259257    for (int i = 0; i < images->n; i++) {
    260         psImage *image = images->data[i]; // Image of interest
    261         assert(image->type.type == PS_TYPE_F32);
    262         assert(scales->data.F32[i] != 0);
    263         if (mask) {
    264             assert(mask->type.type == PS_TYPE_U8);
    265             assert(image->numCols == mask->numCols && image->numRows == mask->numRows);
    266         }
    267         if (errImages) {
    268             psImage *errImage = errImages->data[i];
    269             assert(errImage->type.type == PS_TYPE_F32);
    270             assert(errImage->numCols == image->numCols && errImage->numRows == image->numRows);
    271         }
    272     }
    273 
    274     for (int i = 0; i < images->n; i++) {
    275         psImage *image = images->data[i]; // Image to rescale
    276         psImage *errImage = NULL;       // Variance image to rescale
    277         if (errImages) {
    278             errImage = errImages->data[i];
    279         }
    280         float scale = scales->data.F32[i]; // Scale to use
    281         float offset = offsets->data.F32[i]; // Offset to use
    282         for (int y = 0; y < image->numRows; y++) {
    283             for (int x = 0; x < image->numCols; x++) {
    284                 if (!mask || mask->data.F32[y][x]) {
    285                     image->data.F32[y][x] = (image->data.F32[y][x] - offset) / scale;
    286                     if (errImage) {
    287                         errImage->data.F32[y][x] /= SQUARE(scale);
    288                     }
    289                 }
    290             }
    291         }
     258        psImage *image = images->data[i]; // Image of interest
     259        assert(image->type.type == PS_TYPE_F32);
     260        assert(scales->data.F32[i] != 0);
     261        if (mask) {
     262            assert(mask->type.type == PS_TYPE_U8);
     263            assert(image->numCols == mask->numCols && image->numRows == mask->numRows);
     264        }
     265        if (errImages) {
     266            psImage *errImage = errImages->data[i];
     267            assert(errImage->type.type == PS_TYPE_F32);
     268            assert(errImage->numCols == image->numCols && errImage->numRows == image->numRows);
     269        }
     270    }
     271
     272    for (int i = 0; i < images->n; i++) {
     273        psImage *image = images->data[i]; // Image to rescale
     274        psImage *errImage = NULL;       // Variance image to rescale
     275        if (errImages) {
     276            errImage = errImages->data[i];
     277        }
     278        float scale = scales->data.F32[i]; // Scale to use
     279        float offset = offsets->data.F32[i]; // Offset to use
     280        for (int y = 0; y < image->numRows; y++) {
     281            for (int x = 0; x < image->numCols; x++) {
     282                if (!mask || mask->data.F32[y][x]) {
     283                    image->data.F32[y][x] = (image->data.F32[y][x] - offset) / scale;
     284                    if (errImage) {
     285                        errImage->data.F32[y][x] /= SQUARE(scale);
     286                    }
     287                }
     288            }
     289        }
    292290    }
    293291
Note: See TracChangeset for help on using the changeset viewer.