IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2764


Ignore:
Timestamp:
Dec 20, 2004, 6:15:28 PM (22 years ago)
Author:
Paul Price
Message:

Running faster now by only transforming rejection masks where we're interested. To do: only transform and combine in the second iteration for areas where we're interested

Location:
trunk/stac/src
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/Makefile

    r2751 r2764  
    11SHELL = /bin/sh
    22CC = gcc
    3 CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include # -DCRFLUX -DTESTING
     3CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DCRFLUX -DTESTING
    44PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
    55LDFLAGS += $(PSLIB)
    66
    77OBJECTS = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
    8         stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o
     8        stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o stacAreaOfInterest.o
    99
    1010TARGET = stac

  • trunk/stac/src/stac.c

    r2751 r2764  
    7272#endif
    7373
     74    // Get regions of interest in the source frame
     75    psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest
     76    for (int i = 0; i < inputs->n; i++) {
     77        fprintf(stderr, "Finding area of interest for image %d\n", i);
     78        regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
     79                                              ((psImage*)(inputs->data[i]))->numCols,
     80                                              ((psImage*)(inputs->data[i]))->numRows);
     81#ifdef TESTING
     82        char regionFile[MAXCHAR];       // Filename of region image
     83        sprintf(regionFile,"%s.region",config->inputs->data[i]);
     84        psImageWriteSection(regions->data[i], 0, 0, 0, NULL, 0, regionFile);
     85#endif
     86    }
     87
    7488    // Transform rejected pixels to source frame
    75     psArray *rejectedSource = stacRejection(inputs, transformed, combined, maps, inverseMaps, rejected, config);
     89    psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config);
    7690
    7791    // Redo transformation with the mask
    … …  
    93107    // Free everything I've used
    94108    stacConfigFree(config);
     109    psFree(regions);
    95110    psFree(inputs);
    96111    psFree(maps);

  • trunk/stac/src/stac.h

    r2670 r2764  
    137137// Transform the rejection masks back to the source frame
    138138psArray *stacRejection(psArray *inputs, // Input images
    139                        psArray *transformed, // Transformed images
    140                        psImage *combined, // Combined image
     139                       psArray *rejected, // Rejected images
     140                       psArray *regions, // Regions of interest
    141141                       psArray *maps,   // Maps from input to transformed image
    142142                       psArray *inverseMaps, // Maps from transformed to input image
    143                        psArray *rejected, // Rejected images
    144143                       stacConfig *config // Configuration
    145144    );
    146145
     146/************************************************************************************************************/
     147
     148// stacAreaOfInterest.c
     149
     150// Returns the change in x' and y' for a change in x and y of 1
     151bool stacPlaneTransformDeriv(psPlane *out, // Output derivative, assumed already allocated
     152                             psPlaneTransform *transform, // The transform for which to obtain the derivative
     153                             psPlane *in // The position at which to obtain the derivative
     154    );
     155
     156// Return a mask image designating the region of interest in the source frame
     157// Basically, uses derivatives of the transformation to work out what pixels in the source might be of interest
     158psImage *stacAreaOfInterest(psImage *mask, // Mask that is to be used to determine the area of interest
     159                            psPlaneTransform *map, // Map from the mask to the area of interest
     160                            int nxOut, int nyOut // Size of the area of interest
     161    );
     162
     163/************************************************************************************************************/
    147164
    148165#endif

  • trunk/stac/src/stacCombine.c

    r2670 r2764  
    106106        }
    107107        for (int i = 0; i < nImages; i++) {
    108             (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Create rejection mask
     108            (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8); // Create rejection mask
    109109            // Zero out the mask
    110110            for (int r = 0; r < numRows; r++) {
    111111                for (int c = 0; c < numCols; c++) {
    112                     ((psImage*)((*rejected)->data[i]))->data.F32[r][c] = 0.0;
     112                    ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0;
    113113                }
    114114            }
    … …  
    154154           
    155155            // Rejection iterations
    156             for (int rejNum = 0; rejNum < nReject; rejNum++) {
    157                 float max = 0.0;
    158                 int maxIndex = -1;
     156            bool keepGoing = true;      // Keep going with rejection?
     157            for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
     158                float max = 0.0;        // Maximum deviation
     159                int maxIndex = -1;      // Index of the maximum deviation
    159160                for (int i = 0; i < nImages; i++) {
    160161                    if (!mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
    … …  
    164165                }
    165166                // Reject the pixel with the maximum deviation
    166                 if (max > reject) {   
     167                if (max > reject) {
    167168                    mask->data.U8[maxIndex] = 1;
    168                     ((psImage*)((*rejected)->data[maxIndex]))->data.F32[y][x] += 1.0;
     169                    ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1.0;
    169170                    // Re-do combination following rejection
    170171                    average = stacCombineMean(pixels, deltas, mask);
    171                 }
    172             } // Rejection iterations
    173            
    174             combined->data.F32[y][x] = stacCombineMean(pixels, deltas, mask);
     172                } else {
     173                    keepGoing = false;
     174                }
     175            }
     176           
     177            combined->data.F32[y][x] = average;
    175178        }
    176179    } // Iterating over output pixels

  • trunk/stac/src/stacRead.c

    r2500 r2764  
    8181
    8282        // Mask out the cross-terms
    83         map->x->mask[1][1] = 0;
    84         map->y->coeff[1][1] = 0;
     83        map->x->mask[1][1] = 1;
     84        map->y->mask[1][1] = 1;
    8585        // Set the cross-terms to zero, just in case...
    8686        map->x->coeff[1][1] = 0.0;

  • trunk/stac/src/stacRejection.c

    r2751 r2764  
    4242}
    4343
    44 #if 0
    45 float stacGradient(psImage *image,      // Input for which to measure the gradient
    46                    int x, int y         // Coordinates at which to measure the gradient
    47                    )
    48 {
    49     float sum = 0.0;                    // The sum of surrounding pixels
    50     float maxDiff = 0.0;
    51     int num = 0;
    52     // Get limits
    53     int xMin = MAX(x - 1, 0);
    54     int xMax = MIN(x + 1, image->numCols - 1);
    55     int yMin = MAX(y - 1, 0);
    56     int yMax = MIN(y + 1, image->numRows - 1);
    57     for (int j = yMin; j <= yMax; j++) {
    58         for (int i = xMin; i <= xMax; i++) {
    59             if (image->data.F32[j][i] - image->data.F32[y][x] < maxDiff) {
    60                 maxDiff = image->data.F32[j][i] - image->data.F32[y][x];
    61             }
    62             sum += image->data.F32[j][i];
    63             num++;
    64         }
    65     }
    66     sum -= image->data.F32[y][x];
    67     sum /= (float)(num-1);
    68     sum -= image->data.F32[y][x];
    69 //    return sum / image->data.F32[y][x];
    70     return maxDiff / image->data.F32[y][x];
    71 }   
    72 #endif
    73 
    74 
    7544psArray *stacRejection(psArray *inputs, // Input images
    76                        psArray *transformed, // Transformed images
    77                        psImage *combined, // Combined image
     45                       psArray *rejected, // Rejected images
     46                       psArray *regions, // Regions of interest
    7847                       psArray *maps,   // Maps from input to transformed image
    7948                       psArray *inverseMaps, // Maps from transformed to input image
    80                        psArray *rejected, // Rejected images
    8149                       stacConfig *config // Configuration
    8250    )
    … …  
    8553
    8654    // Check inputs
    87     if (inputs->n != transformed->n) {
    88         psError("stac.rejection", "Number of input (%d) and transformed (%d) images does not match.\n",
    89                 inputs->n, transformed->n);
     55    if (inputs->n != rejected->n) {
     56        psError("stac.rejection", "Number of input images (%d) and rejection masks (%d) does not match.\n",
     57                inputs->n, rejected->n);
     58        return NULL;
     59    }
     60    if (inputs->n != regions->n) {
     61        psError("stac.rejection", "Number of input images (%d) and region masks (%d) does not match.\n",
     62                inputs->n, regions->n);
     63        return NULL;
     64    }
     65    if (inputs->n != maps->n) {
     66        psError("stac.rejection", "Number of input images (%d) and transformation maps (%d) does "
     67                "not match.\n", inputs->n, maps->n);
     68        return NULL;
     69    }
     70    if (inputs->n != inverseMaps->n) {
     71        psError("stac.rejection", "Number of input images (%d) and inverse transformation maps (%d) does "
     72                "not match.\n", inputs->n, inverseMaps->n);
    9073        return NULL;
    9174    }
    9275    for (int i = 0; i < nImages; i++) {
    93         if ((((psImage*)(transformed->data[i]))->numRows != combined->numRows) ||
    94             (((psImage*)(transformed->data[i]))->numCols != combined->numCols)) {
     76        if ((((psImage*)(inputs->data[i]))->numRows != ((psImage*)(regions->data[i]))->numRows) ||
     77            (((psImage*)(inputs->data[i]))->numCols != ((psImage*)(regions->data[i]))->numCols)) {
    9578            psError("stac.rejection",
    96                     "Sizes of transformed (%dx%d) and combined (%dx%d) images do not match.\n",
    97                     ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
    98                     combined->numCols, combined->numRows);
     79                    "Sizes of input image (%dx%d) and region mask (%dx%d) for input %d do not match.\n",
     80                    ((psImage*)(inputs->data[i]))->numCols, ((psImage*)(inputs->data[i]))->numRows,
     81                    ((psImage*)(regions->data[i]))->numCols, ((psImage*)(regions->data[i]))->numRows, i);
    9982            return NULL;
    10083        }
    101         if ((((psImage*)(transformed->data[i]))->numRows != ((psImage*)(rejected->data[i]))->numRows) ||
    102             (((psImage*)(transformed->data[i]))->numCols != ((psImage*)(rejected->data[i]))->numCols)) {
    103             psError("stac.rejection",
    104                     "Sizes of transformed (%dx%d) and rejected (%dx%d) images do not match.\n",
    105                     ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
    106                     ((psImage*)(rejected->data[i]))->numCols, ((psImage*)(rejected->data[i]))->numRows);
    107             return NULL;
    108         }
    109     }
    110 
    111     // Size of output images
    112     int nxOutput = combined->numCols;
    113     int nyOutput = combined->numRows;
     84    }
    11485
    11586    // Stuff for the transformations
    … …  
    132103        psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
    133104        psPlaneTransform *map = maps->data[i]; // The map from input to output
     105        psImage *region = regions->data[i]; // The region of interest for this image
    134106
    135107        psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
    … …  
    153125        for (int y = 0; y < nyInput; y++) {
    154126            for (int x = 0; x < nxInput; x++) {
    155                 inCoords->x = (double)x;
    156                 inCoords->y = (double)y;
    157                 (void)psPlaneTransformApply(outCoords, map, inCoords);
    158                 float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
    159                                                                NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     127
     128                // Only transform pixels of interest
     129                if (region->data.U8[y][x]) {
     130
     131                    inCoords->x = (double)x;
     132                    inCoords->y = (double)y;
     133                    (void)psPlaneTransformApply(outCoords, map, inCoords);
     134                    float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
     135                                                                   NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
    160136#ifdef TESTING
    161                 rejmap->data.F32[y][x] = maskVal;
    162 #endif
    163 
    164 #ifdef TESTING
    165                 // Check gradient
    166                 float meanGrads = 0.0;
    167                 int numGrads = 0;
    168                 for (int j = 0; j < nImages; j++) {
    169                     if (i != j) {
    170                         // Get coordinates for that image
    171                         (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
    172                         int xPix = (int)(inCoords->x + 0.5);
    173                         int yPix = (int)(inCoords->y + 0.5);
    174                         if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
    175                             (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
    176                             // Calculate the gradient
    177                             meanGrads += stacGradient(inputs->data[j], xPix, yPix);
    178                             numGrads++;
     137                    rejmap->data.F32[y][x] = maskVal;
     138
     139                    // Calculate gradient
     140                    float meanGrads = 0.0;
     141                    int numGrads = 0;
     142                    for (int j = 0; j < nImages; j++) {
     143                        if (i != j) {
     144                            // Get coordinates for that image
     145                            (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
     146                            int xPix = (int)(inCoords->x + 0.5);
     147                            int yPix = (int)(inCoords->y + 0.5);
     148                            if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
     149                                (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
     150                                // Calculate the gradient
     151                                meanGrads += stacGradient(inputs->data[j], xPix, yPix);
     152                                numGrads++;
     153                            }
    179154                        }
    180155                    }
    181                 }
    182                 if (numGrads > 0) {
    183                     meanGrads /= (float)numGrads;
    184                 } else {
    185                     meanGrads = 0;
    186                 }
    187 #else
    188                 float meanGrads = 10000.0;
    189 #endif
    190 
    191                
    192 #ifdef TESTING
    193                 grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
    194 #endif
    195                 if (maskVal > config->frac) {
     156                    if (numGrads > 0) {
     157                        meanGrads /= (float)numGrads;
     158                    } else {
     159                        meanGrads = 0;
     160                    }
     161                    grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
     162#endif
    196163
    197164                    if ((maskVal > config->frac) &&
    … …  
    200167                        nBad++;
    201168#ifdef CRFLUX
    202                         fprintf(crs, "%d %d --> %f %f %f\n", x, y, ((psImage*)(inputs->data[i]))->data.F32[y][x],
    203                                 maskVal, stacGradient(inputs->data[i], x, y));
    204 #endif
     169                        fprintf(crs, "%d %d --> %f %f %f\n", x, y,
     170                                ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal,
     171                                stacGradient(inputs->data[i], x, y));
     172#endif
     173                    } else {
     174                        mask->data.U8[y][x] = 0;
    205175                    }
    206                 } else {
    207                     mask->data.U8[y][x] = 0;
    208                 }
     176
     177                } // Only touching pixels of interest
    209178
    210179            }
Note: See TracChangeset for help on using the changeset viewer.