IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2670


Ignore:
Timestamp:
Dec 8, 2004, 3:12:31 PM (22 years ago)
Author:
Paul Price
Message:

More robust gradient measurement

Location:
trunk/stac/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/Makefile

    r2661 r2670  
    11SHELL = /bin/sh
    22CC = gcc
    3 CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DTESTING -DCRFLUX
     3CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DCRFLUX
    44PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
    55LDFLAGS += $(PSLIB)
    … …  
    2121                $(CC) $(CFLAGS) -o $@ calcGrad.o $(LDFLAGS) $(OPTFLAGS)
    2222
     23median:         median.o
     24                $(CC) $(CFLAGS) -o $@ median.o $(LDFLAGS) $(OPTFLAGS)
     25
     26
    2327clean:
    2428                -$(RM) *.o gmon.* profile.txt
    … …  
    3741                -$(RM) test_[0-3].fits.shifterr.*
    3842                -$(RM) leaks.dat
    39                 ./stac -v -k 2.5 -f 0.3 -G 0.6 testout.fits test_0.fits test_1.fits test_2.fits test_3.fits
     43                ./stac -v -k 3 -f 0.7 -G 0.77 testout.fits test_0.fits test_1.fits test_2.fits test_3.fits
    4044
    4145# Run profiling.

  • trunk/stac/src/stac.c

    r2500 r2670  
    5858
    5959    // Transform rejected pixels to source frame
    60     psArray *rejectedSource = stacRejection(inputs, transformed, combined, maps, rejected, config);
     60    psArray *rejectedSource = stacRejection(inputs, transformed, combined, maps, inverseMaps, rejected, config);
    6161
    6262    // Redo transformation with the mask

  • trunk/stac/src/stac.h

    r2661 r2670  
    140140                       psImage *combined, // Combined image
    141141                       psArray *maps,   // Maps from input to transformed image
     142                       psArray *inverseMaps, // Maps from transformed to input image
    142143                       psArray *rejected, // Rejected images
    143144                       stacConfig *config // Configuration

  • trunk/stac/src/stacCombine.c

    r2661 r2670  
    1414#if 0
    1515    // Would like to use psVectorStats, but it doesn't have errors built in yet
    16     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     16    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    1717    (void)psVectorStats(stats, values, masks, 1);
    1818    float mean = stats->sampleMean;
    … …  
    2525    int num = values->n;
    2626    for (int i = 0; i < num; i++) {
    27         if (masks->data.U8[i]) {
     27        if (! masks->data.U8[i]) {
    2828            sum += values->data.F32[i] / SQUARE(errors->data.F32[i]);
    2929            weights += 1.0 / SQUARE(errors->data.F32[i]);
    … …  
    4646{
    4747    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    48     (void)psVectorStats(stats, values, masks, 0);
     48    (void)psVectorStats(stats, values, masks, 1);
    4949    float median = stats->sampleMedian;
    5050    psFree(stats);
    … …  
    132132                deltas->data.F32[i] = delta;
    133133                if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
    134                     mask->data.U8[i] = (psU8)0; // Don't use!
     134                    mask->data.U8[i] = (psU8)1; // Don't use!
    135135                } else {
    136                     mask->data.U8[i] = (psU8)1; // Use.
     136                    mask->data.U8[i] = (psU8)0; // Use.
    137137                }
    138138            }
    139139           
    140             float average = stacCombineMedian(pixels, deltas, mask); // Combined value
     140            float average = stacCombineMean(pixels, deltas, mask); // Combined value
    141141
    142142#ifdef TESTING
    … …  
    145145            int numGoodPix = 0;
    146146            for (int i = 0; i < nImages; i++) {
    147                 if (mask->data.U8[i]) {
     147                if (! mask->data.U8[i]) {
    148148                    chi2->data.F32[y][x] += SQUARE((pixels->data.F32[i] - average) / deltas->data.F32[i]);
    149149                    numGoodPix++;
    … …  
    158158                int maxIndex = -1;
    159159                for (int i = 0; i < nImages; i++) {
    160                     if (mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
     160                    if (!mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
    161161                        max = (pixels->data.F32[i] - average) / deltas->data.F32[i];
    162162                        maxIndex = i;
    … …  
    165165                // Reject the pixel with the maximum deviation
    166166                if (max > reject) {   
    167                     mask->data.U8[maxIndex] = 0;
     167                    mask->data.U8[maxIndex] = 1;
    168168                    ((psImage*)((*rejected)->data[maxIndex]))->data.F32[y][x] += 1.0;
    169169                    // Re-do combination following rejection

  • trunk/stac/src/stacRejection.c

    r2661 r2670  
    1111    )
    1212{
    13     float sum = 0.0;                    // The sum of surrounding pixels
    1413    int num = 0;
    15     psVector *pixels = psVectorAlloc(9, PS_TYPE_F32); // Array of pixels
    16     psVector *mask = psVectorAlloc(9, PS_TYPE_U8); // Corresponding mask
     14    psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels
     15    psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask
    1716
    1817    // Get limits
    … …  
    2322    for (int j = yMin; j <= yMax; j++) {
    2423        for (int i = xMin; i <= xMax; i++) {
    25             sum += image->data.F32[j][i];
    26             pixels->data.F32[num] = image->data.F32[j][i];
    27             mask->data.U8[num] = 1;
    28             num++;
    29         }
    30     }
    31 #if 0
    32     sum -= image->data.F32[y][x];
    33     sum /= (float)(num-1);
    34 //    sum -= image->data.F32[y][x];
    35     return sum / image->data.F32[y][x];
    36 #endif
    37 
    38     // Fill out the array
    39     for (int i = num; i < 9; i++) {
    40         pixels->data.F32[i] = 0;
    41         mask->data.U8[i] = 0;
    42     }
     24            if ((i != x) && (j != y)) {
     25                pixels->data.F32[num] = image->data.F32[j][i];
     26                mask->data.U8[num] = 0;
     27                num++;
     28            }
     29        }
     30    }
     31    pixels->n = num;
     32    mask->n = num;
     33
     34    // Get the median
    4335    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    44     (void)psVectorStats(stats, pixels, mask, 0);
     36    (void)psVectorStats(stats, pixels, mask, 1);
    4537    float median = stats->sampleMedian;
    4638    psFree(stats);
    … …  
    8577                       psImage *combined, // Combined image
    8678                       psArray *maps,   // Maps from input to transformed image
     79                       psArray *inverseMaps, // Maps from transformed to input image
    8780                       psArray *rejected, // Rejected images
    8881                       stacConfig *config // Configuration
    … …  
    167160#ifdef TESTING
    168161                rejmap->data.F32[y][x] = maskVal;
    169                 grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y);
    170 #endif
    171                 // Check threshold, then check gradient
    172                 if ((maskVal > config->frac) && (stacGradient(inputs->data[i], x, y) < config->grad)) {
    173                     mask->data.U8[y][x] = 1;
    174                     nBad++;
     162#endif
     163
     164                // Check gradient
     165                float meanGrads = 0.0;
     166                int numGrads = 0;
     167                for (int j = 0; j < nImages; j++) {
     168                    if (i != j) {
     169                        // Get coordinates for that image
     170                        (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
     171                        int xPix = (int)(inCoords->x + 0.5);
     172                        int yPix = (int)(inCoords->y + 0.5);
     173                        if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
     174                            (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
     175                            // Calculate the gradient
     176                            meanGrads += stacGradient(inputs->data[j], xPix, yPix);
     177                            numGrads++;
     178                        }
     179                    }
     180                }
     181                if (numGrads > 0) {
     182                    meanGrads /= (float)numGrads;
     183                } else {
     184                    meanGrads = 0;
     185                }
     186               
     187#ifdef TESTING
     188                grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
     189#endif
     190                if (maskVal > config->frac) {
     191
     192                    if ((maskVal > config->frac) &&
     193                        (stacGradient(inputs->data[i], x, y)) < config->grad * meanGrads) {
     194                        mask->data.U8[y][x] = 1;
     195                        nBad++;
    175196#ifdef CRFLUX
    176                     fprintf(crs, "%d %d --> %f %f %f\n", x, y, ((psImage*)(inputs->data[i]))->data.F32[y][x],
    177                             maskVal, stacGradient(inputs->data[i], x, y));
     197                        fprintf(crs, "%d %d --> %f %f %f\n", x, y, ((psImage*)(inputs->data[i]))->data.F32[y][x],
     198                                maskVal, stacGradient(inputs->data[i], x, y));
     199                    }
    178200#endif
    179201                } else {
    … …  
    189211#endif
    190212
    191         psTrace("stac.rejection", 1, "%d pixels masked in image %d.\n", nBad, i);
     213        psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
    192214        // Clip the image, and convert to suitable mask format
    193215
Note: See TracChangeset for help on using the changeset viewer.