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/stacCombine.c

    r5743 r5745  
    88#define ABS(x) ((x) >= 0 ? (x) : -(x))
    99
    10 float stacCombineMean(psVector *values, // Values for which to take the mean
    11                       psVector *errors, // Errors in the values
    12                       psVector *masks   // Masks for the values, 1 = don't use, 0 = use
     10float stacCombineMean(psVector *values, // Values for which to take the mean
     11                      psVector *errors, // Errors in the values
     12                      psVector *masks   // Masks for the values, 1 = don't use, 0 = use
    1313    )
    1414{
    … …  
    2626    int num = values->n;
    2727    for (int i = 0; i < num; i++) {
    28         if (! masks->data.U8[i]) {
    29             // "error" here is the variance
    30             sum += values->data.F32[i] / errors->data.F32[i];
    31             weights += 1.0 / errors->data.F32[i];
    32         }
     28        if (! masks->data.U8[i]) {
     29            // "error" here is the variance
     30            sum += values->data.F32[i] / errors->data.F32[i];
     31            weights += 1.0 / errors->data.F32[i];
     32        }
    3333    }
    3434    if (weights > 0.0) {
    35         return (float)(sum/weights);
     35        return (float)(sum/weights);
    3636    } else {
    37         return NAN;
     37        return NAN;
    3838    }
    3939#endif
    … …  
    4343
    4444float stacCombineMedian(psVector *values, // Values for which to take the median
    45                         psVector *errors, // Errors in the values
    46                         psVector *masks // Masks for the values, 0 = don't use, 1 = use
     45                        psVector *errors, // Errors in the values
     46                        psVector *masks // Masks for the values, 0 = don't use, 1 = use
    4747    )
    4848{
    … …  
    5656
    5757
    58 bool stacCombine(psImage **combined,    // The combined image for output
    59                  psArray **rejected,    // Array of rejection masks
    60                  psArray *images,       // Array of transformed images
    61                  psArray *errors,       // Array of transformed error images
    62                  int nReject,           // Number of rejection iterations
    63                  psImage *region,       // Region to combine
    64                  psVector *saturated,   // Saturation limits for each image
    65                  psVector *bad,         // Bad pixel limits for each image
    66                  float reject           // Rejection (k-sigma)
     58bool stacCombine(psImage **combined,    // The combined image for output
     59                 psArray **rejected,    // Array of rejection masks
     60                 psArray *images,       // Array of transformed images
     61                 psArray *errors,       // Array of transformed error images
     62                 int nReject,           // Number of rejection iterations
     63                 psImage *region,       // Region to combine
     64                 psVector *saturated,   // Saturation limits for each image
     65                 psVector *bad,         // Bad pixel limits for each image
     66                 float reject           // Rejection (k-sigma)
    6767    )
    6868{
    … …  
    7676    assert(bad->n == images->n);
    7777
    78     int nImages = images->n;            // Number of images
    79     int numRows = ((psImage*)images->data[0])->numRows; // Image size
     78    int nImages = images->n;            // Number of images
     79    int numRows = ((psImage*)images->data[0])->numRows; // Image size
    8080    int numCols = ((psImage*)images->data[0])->numCols; // Image size
    81    
     81
    8282    // Check dimensions for consistency
    8383    for (int i = 0; i < nImages; i++) {
    84         psImage *image = (psImage *)images->data[i]; // The image
    85         psImage *error = (psImage *)errors->data[i]; // The error image
    86 
    87         assert(image->numCols == numCols && image->numRows == numRows);
    88         assert(error->numCols == numCols && error->numRows == numRows);
     84        psImage *image = (psImage *)images->data[i]; // The image
     85        psImage *error = (psImage *)errors->data[i]; // The error image
     86
     87        assert(image->numCols == numCols && image->numRows == numRows);
     88        assert(error->numCols == numCols && error->numRows == numRows);
    8989    }
    9090
    9191    // Check combined image
    92     assert(!*combined || ((*combined)->numRows == numRows) && ((*combined)->numCols == numCols));
     92    assert(!*combined || ((*combined)->numRows == numRows && (*combined)->numCols == numCols));
    9393    if (*combined == NULL) {
    94         *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
     94        *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
    9595    }
    9696
    9797    // Check area of interest
    98     assert(!region || (region->numRows == numRows) && (region->numCols == numCols));
     98    assert(!region || (region->numRows == numRows && region->numCols == numCols));
    9999
    100100    psTrace("stac.combine", 1, "Combining images....\n");
    … …  
    106106    // Set up rejection masks
    107107    if (nReject > 0) {
    108         if (*rejected == NULL) {
    109             // Allocate the rejection masks, if required
    110             *rejected = psArrayAlloc(nImages);
    111         } else {
    112             assert((*rejected)->n != nImages);
    113         }
    114 
    115         // Create and initialise rejection masks
    116         for (int i = 0; i < nImages; i++) {
    117             (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
    118             for (int r = 0; r < numRows; r++) {
    119                 for (int c = 0; c < numCols; c++) {
    120                     ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0;
    121                 }
    122             }
    123         }
     108        if (*rejected == NULL) {
     109            // Allocate the rejection masks, if required
     110            *rejected = psArrayAlloc(nImages);
     111        } else {
     112            assert((*rejected)->n != nImages);
     113        }
     114
     115        // Create and initialise rejection masks
     116        for (int i = 0; i < nImages; i++) {
     117            (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
     118            for (int r = 0; r < numRows; r++) {
     119                for (int c = 0; c < numCols; c++) {
     120                    ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0;
     121                }
     122            }
     123        }
    124124    }
    125125
    … …  
    127127    // chi^2 image
    128128    psImage *chi2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    129     static int iteration = 0;           // Number of times function has been called
    130 #endif
    131    
     129    static int iteration = 0;           // Number of times function has been called
     130#endif
     131
    132132    for (int y = 0; y < numRows; y++) {
    133         for (int x = 0; x < numCols; x++) {
    134 
    135             // Only combine those pixels requested
    136             if (!region || (region && region->data.U8[y][x])) {
    137            
    138                 // Export pixels into the vector and get stats
    139                 for (int i = 0; i < nImages; i++) {
    140                     float pixel = ((psImage*)images->data[i])->data.F32[y][x];
    141                     float delta = ((psImage*)errors->data[i])->data.F32[y][x]; // This is the variance!
    142                     pixels->data.F32[i] = pixel;
    143                     deltas->data.F32[i] = delta;
    144                     if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
    145                         (! isfinite(pixel)) || (! isfinite(delta))) {
    146                         mask->data.U8[i] = (psU8)1; // Don't use!
    147                     } else {
    148                         mask->data.U8[i] = (psU8)0; // Use.
    149                     }
    150                 }
    151                
    152                 float average = stacCombineMean(pixels, deltas, mask); // Combined value
    153 
    154                 // We set the value BEFORE the rejection iteration because not all pixels that we reject
    155                 // here will be rejected in the final cut.
    156                 (*combined)->data.F32[y][x] = average;
     133        for (int x = 0; x < numCols; x++) {
     134
     135            // Only combine those pixels requested
     136            if (!region || (region && region->data.U8[y][x])) {
     137
     138                // Export pixels into the vector and get stats
     139                for (int i = 0; i < nImages; i++) {
     140                    float pixel = ((psImage*)images->data[i])->data.F32[y][x];
     141                    float delta = ((psImage*)errors->data[i])->data.F32[y][x]; // This is the variance!
     142                    pixels->data.F32[i] = pixel;
     143                    deltas->data.F32[i] = delta;
     144                    if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
     145                        (! isfinite(pixel)) || (! isfinite(delta))) {
     146                        mask->data.U8[i] = (psU8)1; // Don't use!
     147                    } else {
     148                        mask->data.U8[i] = (psU8)0; // Use.
     149                    }
     150                }
     151
     152                float average = stacCombineMean(pixels, deltas, mask); // Combined value
     153
     154                // We set the value BEFORE the rejection iteration because not all pixels that we reject
     155                // here will be rejected in the final cut.
     156                (*combined)->data.F32[y][x] = average;
    157157
    158158#ifdef TESTING
    159                 // Calculate chi^2
    160                 chi2->data.F32[y][x] = 0.0;
    161                 int numGoodPix = 0;
    162                 for (int i = 0; i < nImages; i++) {
    163                     if (! mask->data.U8[i]) {
    164                         chi2->data.F32[y][x] += SQUARE(pixels->data.F32[i] - average) / deltas->data.F32[i];
    165                         numGoodPix++;
    166                     }
    167                 }
    168                 chi2->data.F32[y][x] /= (float)numGoodPix;
    169 #endif
    170                
    171                 // Rejection iterations
    172                 bool keepGoing = true;  // Keep going with rejection?
    173                 for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
    174                     float max = 0.0;    // Maximum deviation
    175                     int maxIndex = -1;  // Index of the maximum deviation
    176                     for (int i = 0; i < nImages; i++) {
    177                         if (!mask->data.U8[i] &&
    178                             ((pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]) > max)) {
    179                             max = (pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]);
    180                             maxIndex = i;
    181                         }
    182                     }
    183                     // Reject the pixel with the maximum deviation
    184                     if (max > reject) {
    185                         mask->data.U8[maxIndex] = 1;
    186                         ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1;
    187                         // Re-do combination following rejection
    188                         average = stacCombineMean(pixels, deltas, mask);
    189                     } else {
    190                         keepGoing = false;
    191                     }
    192                 }
    193            
    194             } // Pixels of interest
    195 
    196         }
     159                // Calculate chi^2
     160                chi2->data.F32[y][x] = 0.0;
     161                int numGoodPix = 0;
     162                for (int i = 0; i < nImages; i++) {
     163                    if (! mask->data.U8[i]) {
     164                        chi2->data.F32[y][x] += SQUARE(pixels->data.F32[i] - average) / deltas->data.F32[i];
     165                        numGoodPix++;
     166                    }
     167                }
     168                chi2->data.F32[y][x] /= (float)numGoodPix;
     169#endif
     170
     171                // Rejection iterations
     172                bool keepGoing = true;  // Keep going with rejection?
     173                for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
     174                    float max = 0.0;    // Maximum deviation
     175                    int maxIndex = -1;  // Index of the maximum deviation
     176                    for (int i = 0; i < nImages; i++) {
     177                        if (!mask->data.U8[i] &&
     178                            ((pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]) > max)) {
     179                            max = (pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]);
     180                            maxIndex = i;
     181                        }
     182                    }
     183                    // Reject the pixel with the maximum deviation
     184                    if (max > reject) {
     185                        mask->data.U8[maxIndex] = 1;
     186                        ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1;
     187                        // Re-do combination following rejection
     188                        average = stacCombineMean(pixels, deltas, mask);
     189                    } else {
     190                        keepGoing = false;
     191                    }
     192                }
     193
     194            } // Pixels of interest
     195
     196        }
    197197    } // Iterating over output pixels
    198198
    … …  
    200200    // Write chi^2 image
    201201    iteration++;
    202     char chiName[MAXCHAR];              // Filename of chi^2 image
     202    char chiName[MAXCHAR];              // Filename of chi^2 image
    203203    sprintf(chiName,"chi2_%d.fits",iteration);
    204204    psFits *chiFile = psFitsAlloc(chiName);
    205205    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0)) {
    206         psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
     206        psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
    207207    }
    208208    psTrace("stac.combine", 1, "Chi^2 image written to %s\n", chiName);
Note: See TracChangeset for help on using the changeset viewer.