IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 18, 2006, 12:20:45 PM (20 years ago)
Author:
Paul Price
Message:

Updating to use psLib rel11. Main problem was that vector and array lengths ('n' element) are no longer set to equal the number of allocated values ('nalloc' element) when allocated.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/stacRejection.c

    r5743 r6887  
    88#define MIN(x,y) ((x) < (y) ? (x) : (y))
    99
    10 float stacGradient(psImage *image,      // Input for which to measure the gradient
    11                    int x, int y         // Coordinates at which to measure the gradient
     10float stacGradient(psImage *image,      // Input for which to measure the gradient
     11                   int x, int y         // Coordinates at which to measure the gradient
    1212    )
    1313{
    … …  
    2222    int yMax = MIN(y + 1, image->numRows - 1);
    2323    for (int j = yMin; j <= yMax; j++) {
    24         for (int i = xMin; i <= xMax; i++) {
    25             if ((i != x) && (j != y)) {
    26                 pixels->data.F32[num] = image->data.F32[j][i];
    27                 mask->data.U8[num] = 0;
    28                 num++;
    29             }
    30         }
     24        for (int i = xMin; i <= xMax; i++) {
     25            if ((i != x) && (j != y)) {
     26                pixels->data.F32[num] = image->data.F32[j][i];
     27                mask->data.U8[num] = 0;
     28                num++;
     29            }
     30        }
    3131    }
    3232    pixels->n = num;
    … …  
    4343}
    4444
    45 psArray *stacRejection(psArray *inputs, // Input images
    46                        psArray *rejected, // Rejected images
    47                        psArray *regions, // Regions of interest
    48                        psArray *maps,   // Maps from input to transformed image
    49                        psArray *inverseMaps, // Maps from transformed to input image
    50                        float frac,      // Fraction of pixel rejected before looking more carefully
    51                        float grad,      // Gradient limit for rejection
    52                        psArray *names   // Names of original images (only for writing out when TESTING)
     45psArray *stacRejection(psArray *inputs, // Input images
     46                       psArray *rejected, // Rejected images
     47                       psArray *regions, // Regions of interest
     48                       psArray *maps,   // Maps from input to transformed image
     49                       psArray *inverseMaps, // Maps from transformed to input image
     50                       float frac,      // Fraction of pixel rejected before looking more carefully
     51                       float grad,      // Gradient limit for rejection
     52                       psArray *names   // Names of original images (only for writing out when TESTING)
    5353    )
    5454{
    55     int nImages = inputs->n;            // Number of input images
     55    int nImages = inputs->n;            // Number of input images
    5656
    5757    // Check inputs
    … …  
    6363
    6464    for (int i = 0; i < nImages; i++) {
    65         psImage *input = inputs->data[i];
    66         psImage *region = regions->data[i];
    67         assert(input->numRows == region->numRows && input->numCols == region->numCols);
     65        psImage *input = inputs->data[i];
     66        psImage *region = regions->data[i];
     67        assert(input->numRows == region->numRows && input->numCols == region->numCols);
    6868    }
    6969
    … …  
    7777    psVector *grads = psVectorAlloc(nImages, PS_TYPE_F32); // Gradient for each image
    7878    psVector *gradsMask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask for gradient vector
     79    grads->n = nImages;
     80    gradsMask->n = nImages;
    7981
    8082    // Transform rejection masks back to source
    8183    psArray *inputRej = psArrayAlloc(nImages);
     84    inputRej->n = nImages;
    8285    for (int i = 0; i < nImages; i++) {
    83         // Size of input image
    84         int nxInput = ((psImage*)(inputs->data[i]))->numCols;
    85         int nyInput = ((psImage*)(inputs->data[i]))->numRows;
    86         psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_U8); // The pixel mask for the input
    87 #ifdef TESTING
    88         psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
    89         psImage *gradient = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);        // The gradient image
    90 #endif
    91         psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
    92         psPlaneTransform *map = maps->data[i]; // The map from input to output
    93         psImage *region = regions->data[i]; // The region of interest for this image
    94 
    95         psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
    96         int nBad = 0;                   // Number of bad pixels
     86        // Size of input image
     87        int nxInput = ((psImage*)(inputs->data[i]))->numCols;
     88        int nyInput = ((psImage*)(inputs->data[i]))->numRows;
     89        psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_U8); // The pixel mask for the input
     90#ifdef TESTING
     91        psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
     92        psImage *gradient = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);        // The gradient image
     93#endif
     94        psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
     95        psPlaneTransform *map = maps->data[i]; // The map from input to output
     96        psImage *region = regions->data[i]; // The region of interest for this image
     97
     98        psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
     99        int nBad = 0;                   // Number of bad pixels
    97100
    98101#ifdef CRFLUX
    99         // Set up CR output
    100         FILE *crs = NULL;               // File for outputting details of rejected pixels
    101         char crfile[MAXCHAR];   // Filename
    102         sprintf(crfile,"%s.cr",names->data[i]);
    103         if ((crs = fopen(crfile, "w")) == NULL) {
    104             fprintf(stderr, "Unable to open file for detailed output, %s\n");
    105             return NULL;
    106         }
    107 #endif
    108 
    109         // Transform the mask
    110         // Optimisation option is to only transform the pixels that have been rejected in the output,
    111         // calculate derivatives of the map, and use that as a buffer around the transformed position
    112         // in the input image.
    113         psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    114         for (int y = 0; y < nyInput; y++) {
    115             for (int x = 0; x < nxInput; x++) {
    116 
    117                 // Only transform pixels of interest
    118                 if (region->data.U8[y][x]) {
    119 
    120                     inCoords->x = (double)x + 0.5;
    121                     inCoords->y = (double)y + 0.5;
    122                     (void)psPlaneTransformApply(outCoords, map, inCoords);
    123                     float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
    124                                                                    NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
    125 #ifdef TESTING
    126                     rejmap->data.F32[y][x] = maskVal;
    127 #endif
    128 
    129                     if (maskVal > frac) {
    130                         // Calculate mean gradient on other images
    131                         float meanGrads = 0.0;
    132                         int numGrads = 0;
    133                         for (int j = 0; j < nImages; j++) {
    134                             if (i != j) {
    135                                 // Get coordinates for that image
    136                                 (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
    137                                 int xPix = (int)(inCoords->x + 0.5);
    138                                 int yPix = (int)(inCoords->y + 0.5);
    139                                 if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
    140                                     (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
    141                                     // Calculate the gradient
    142                                     grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
    143                                     gradsMask->data.U8[j] = 0;
    144                                     numGrads++;
    145                                 } else {
    146                                     gradsMask->data.U8[j] = 1; // Mask this one
    147                                 }
    148                             } else {
    149                                 gradsMask->data.U8[j] = 1; // Mask this one
    150                             }
    151                         }
    152                         if (numGrads > 0) {
    153                             (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
    154                             meanGrads = median->sampleMedian;
    155                         } else {
    156                             meanGrads = 0.0;
    157                         }
    158 
    159 #ifdef TESTING
    160                         //gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
    161                         gradient->data.F32[y][x] = meanGrads;
    162 #endif
    163 
    164                         //if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
    165                         if (meanGrads > grad) {
    166                             mask->data.U8[y][x] = 1;
    167                             nBad++;
     102        // Set up CR output
     103        FILE *crs = NULL;               // File for outputting details of rejected pixels
     104        char crfile[MAXCHAR];   // Filename
     105        sprintf(crfile,"%s.cr",names->data[i]);
     106        if ((crs = fopen(crfile, "w")) == NULL) {
     107            fprintf(stderr, "Unable to open file for detailed output, %s\n");
     108            return NULL;
     109        }
     110#endif
     111
     112        // Transform the mask
     113        // Optimisation option is to only transform the pixels that have been rejected in the output,
     114        // calculate derivatives of the map, and use that as a buffer around the transformed position
     115        // in the input image.
     116        psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     117        for (int y = 0; y < nyInput; y++) {
     118            for (int x = 0; x < nxInput; x++) {
     119
     120                // Only transform pixels of interest
     121                if (region->data.U8[y][x]) {
     122
     123                    inCoords->x = (double)x + 0.5;
     124                    inCoords->y = (double)y + 0.5;
     125                    (void)psPlaneTransformApply(outCoords, map, inCoords);
     126                    float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
     127                                                                   NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     128#ifdef TESTING
     129                    rejmap->data.F32[y][x] = maskVal;
     130#endif
     131
     132                    if (maskVal > frac) {
     133                        // Calculate mean gradient on other images
     134                        float meanGrads = 0.0;
     135                        int numGrads = 0;
     136                        for (int j = 0; j < nImages; j++) {
     137                            if (i != j) {
     138                                // Get coordinates for that image
     139                                (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
     140                                int xPix = (int)(inCoords->x + 0.5);
     141                                int yPix = (int)(inCoords->y + 0.5);
     142                                if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
     143                                    (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
     144                                    // Calculate the gradient
     145                                    grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
     146                                    gradsMask->data.U8[j] = 0;
     147                                    numGrads++;
     148                                } else {
     149                                    gradsMask->data.U8[j] = 1; // Mask this one
     150                                }
     151                            } else {
     152                                gradsMask->data.U8[j] = 1; // Mask this one
     153                            }
     154                        }
     155                        if (numGrads > 0) {
     156                            (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
     157                            meanGrads = median->sampleMedian;
     158                        } else {
     159                            meanGrads = 0.0;
     160                        }
     161
     162#ifdef TESTING
     163                        //gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
     164                        gradient->data.F32[y][x] = meanGrads;
     165#endif
     166
     167                        //if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
     168                        if (meanGrads > grad) {
     169                            mask->data.U8[y][x] = 1;
     170                            nBad++;
    168171
    169172#ifdef CRFLUX
    170                             fprintf(crs, "%d %d --> %f %f %f\n", x, y,
    171                                     ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal,
    172                                     stacGradient(inputs->data[i], x, y));
    173 #endif
    174 
    175                         } else {
    176                             mask->data.U8[y][x] = 0;
    177                         } // Gradient threshold
    178                     } else {
    179                         mask->data.U8[y][x] = 0;
    180                     } // Rejection threshold
    181 
    182                 } else {
    183                     mask->data.U8[y][x] = 0;
    184                 } // Only touching pixels of interest
    185 
    186             }
    187         } // Iterating over pixels
    188         psFree(median);
     173                            fprintf(crs, "%d %d --> %f %f %f\n", x, y,
     174                                    ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal,
     175                                    stacGradient(inputs->data[i], x, y));
     176#endif
     177
     178                        } else {
     179                            mask->data.U8[y][x] = 0;
     180                        } // Gradient threshold
     181                    } else {
     182                        mask->data.U8[y][x] = 0;
     183                    } // Rejection threshold
     184
     185                } else {
     186                    mask->data.U8[y][x] = 0;
     187                } // Only touching pixels of interest
     188
     189            }
     190        } // Iterating over pixels
     191        psFree(median);
    189192
    190193#ifdef CRFLUX
    191         // Close file used for detailed output
    192         fclose(crs);
    193 #endif
    194 
    195         psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
    196         // Clip the image, and convert to suitable mask format
    197 
    198 #ifdef TESTING
    199         // Write error image out to check
    200         char maskName[MAXCHAR];         // Filename of mask image
    201         char rejmapName[MAXCHAR];       // Filename of rejection image
    202         char gradName[MAXCHAR];         // Filename of gradient image
    203         sprintf(maskName, "%s.mask", names->data[i]);
    204         sprintf(rejmapName, "%s.rejmap", names->data[i]);
    205         sprintf(gradName, "%s.grad", names->data[i]);
    206 
    207         psFits *maskFile = psFitsAlloc(maskName);
    208         psFits *rejmapFile = psFitsAlloc(rejmapName);
    209         psFits *gradFile = psFitsAlloc(gradName);
    210         if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
    211             psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
    212         }
    213         psTrace("stac", 1, "Mask image written to %s\n", maskName);
    214         if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0)) {
    215             psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
    216         }
    217         psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
    218         if (!psFitsWriteImage(gradFile, NULL, gradient, 0)) {
    219             psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
    220         }
    221         psTrace("stac", 1, "Gradient image written to %s\n", gradName);
    222         psFree(maskFile);
    223         psFree(rejmapFile);
    224         psFree(gradFile);
    225         psFree(rejmap);
    226         psFree(gradient);
    227 #endif
    228 
    229         // Stuff into the array
    230         inputRej->data[i] = mask;
     194        // Close file used for detailed output
     195        fclose(crs);
     196#endif
     197
     198        psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
     199        // Clip the image, and convert to suitable mask format
     200
     201#ifdef TESTING
     202        // Write error image out to check
     203        char maskName[MAXCHAR];         // Filename of mask image
     204        char rejmapName[MAXCHAR];       // Filename of rejection image
     205        char gradName[MAXCHAR];         // Filename of gradient image
     206        sprintf(maskName, "%s.mask", (char*)names->data[i]);
     207        sprintf(rejmapName, "%s.rejmap", (char*)names->data[i]);
     208        sprintf(gradName, "%s.grad", (char*)names->data[i]);
     209
     210        psFits *maskFile = psFitsOpen(maskName, "w");
     211        psFits *rejmapFile = psFitsOpen(rejmapName, "w");
     212        psFits *gradFile = psFitsOpen(gradName, "w");
     213        if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
     214            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
     215        }
     216        psTrace("stac", 1, "Mask image written to %s\n", maskName);
     217        if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0, NULL)) {
     218            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
     219        }
     220        psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
     221        if (!psFitsWriteImage(gradFile, NULL, gradient, 0, NULL)) {
     222            psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
     223        }
     224        psTrace("stac", 1, "Gradient image written to %s\n", gradName);
     225        psFitsClose(maskFile);
     226        psFitsClose(rejmapFile);
     227        psFitsClose(gradFile);
     228        psFree(rejmap);
     229        psFree(gradient);
     230#endif
     231
     232        // Stuff into the array
     233        inputRej->data[i] = mask;
    231234
    232235    }
Note: See TracChangeset for help on using the changeset viewer.