IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2661


Ignore:
Timestamp:
Dec 7, 2004, 3:38:05 PM (22 years ago)
Author:
Paul Price
Message:

Working, but not getting faint CRs

Location:
trunk/stac/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/Makefile

    r2500 r2661  
    11SHELL = /bin/sh
    22CC = gcc
    3 CFLAGS = -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include # -DTESTING
    4 PSLIB = -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
    5 LDFLAGS = $(PSLIB)
     3CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DTESTING -DCRFLUX
     4PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
     5LDFLAGS += $(PSLIB)
    66
    77OBJECTS = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
    … …  
    1818                $(CC) $(CFLAGS) -o $@ $(OBJECTS) $(LDFLAGS) $(OPTFLAGS)
    1919
     20calcGrad:       calcGrad.o
     21                $(CC) $(CFLAGS) -o $@ calcGrad.o $(LDFLAGS) $(OPTFLAGS)
     22
    2023clean:
    2124                -$(RM) *.o gmon.* profile.txt
    … …  
    2730                -$(RM) testout.fits
    2831                -$(RM) testout.fits.pre
     32                -$(RM) chi2_*.fits
    2933                -$(RM) test_[0-3].fits.err
    3034                -$(RM) test_[0-3].fits.mask
    31                 -$(RM) test_[0-3].fits.shift.?
     35                -$(RM) test_[0-3].fits.shift.*
    3236                -$(RM) test_[0-3].fits.shiftrej
    3337                -$(RM) test_[0-3].fits.shifterr.*
    3438                -$(RM) leaks.dat
    35                 ./stac -v testout.fits test_0.fits test_1.fits test_2.fits test_3.fits
     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
    3640
    3741# Run profiling.
    … …  
    4953optimise:       clean
    5054                -$(RM) $(TARGET)
    51                 export OPTFLAGS=-fprofile-arcs ; $(MAKE) test
     55                export OPTFLAGS=-fprofile-generate ; $(MAKE) test
    5256                -$(RM) $(TARGET)
    5357                $(MAKE) clean
    54                 export OPTFLAGS=-fbranch-probabilities ; $(MAKE) $(TARGET)
     58                export OPTFLAGS=-fbranch-use ; $(MAKE) $(TARGET)
    5559                -$(RM) *.da
    5660

  • trunk/stac/src/stac.h

    r2500 r2661  
    2727    float frac;                         // Fraction of input pixel that must be masked before the pixel is
    2828                                        // considered bad
     29    float grad;                         // Multiplier of the gradient
    2930} stacConfig;
    3031
    … …  
    129130// stacRejection.c
    130131
     132// Return the relative gradient for a given pixel
     133float stacGradient(psImage *image,      // Input for which to measure the gradient
     134                   int x, int y         // Coordinates at which to measure the gradient
     135    );
     136
    131137// Transform the rejection masks back to the source frame
    132138psArray *stacRejection(psArray *inputs, // Input images

  • trunk/stac/src/stacCombine.c

    r2500 r2661  
    1616    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    1717    (void)psVectorStats(stats, values, masks, 1);
     18    float mean = stats->sampleMean;
    1819    psFree(stats);
    19     return stats->sampleMean;
     20    return mean;
    2021#else
    2122    // Instead, do it ourselves
    … …  
    4546{
    4647    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    47     (void)psVectorStats(stats, values, masks, 1);
     48    (void)psVectorStats(stats, values, masks, 0);
     49    float median = stats->sampleMedian;
    4850    psFree(stats);
    49     return stats->sampleMedian;
     51    return median;
    5052}
    5153
    … …  
    113115        }
    114116    }
    115    
     117
     118#ifdef TESTING
     119    // chi^2 image
     120    psImage *chi2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     121    static int iteration = 0;           // Number of times function has been called
     122#endif
    116123   
    117124    for (int y = 0; y < numRows; y++) {
    … …  
    131138            }
    132139           
    133             float average = stacCombineMean(pixels, deltas, mask); // Combined value
     140            float average = stacCombineMedian(pixels, deltas, mask); // Combined value
     141
     142#ifdef TESTING
     143            // Calculate chi^2
     144            chi2->data.F32[y][x] = 0.0;
     145            int numGoodPix = 0;
     146            for (int i = 0; i < nImages; i++) {
     147                if (mask->data.U8[i]) {
     148                    chi2->data.F32[y][x] += SQUARE((pixels->data.F32[i] - average) / deltas->data.F32[i]);
     149                    numGoodPix++;
     150                }
     151            }
     152            chi2->data.F32[y][x] /= (float)numGoodPix;
     153#endif
    134154           
    135155            // Rejection iterations
    136156            for (int rejNum = 0; rejNum < nReject; rejNum++) {
    137157                float max = 0.0;
    138                 int maxIndex = 0;
     158                int maxIndex = -1;
    139159                for (int i = 0; i < nImages; i++) {
    140                     if (mask->data.U8[i] && (ABS(pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
    141                         max = ABS(pixels->data.F32[i] - average) / deltas->data.F32[i];
     160                    if (mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
     161                        max = (pixels->data.F32[i] - average) / deltas->data.F32[i];
    142162                        maxIndex = i;
    143163                    }
    … …  
    152172            } // Rejection iterations
    153173           
    154             combined->data.F32[y][x] = average;
     174            combined->data.F32[y][x] = stacCombineMean(pixels, deltas, mask);
    155175        }
    156176    } // Iterating over output pixels
    … …  
    160180    if (nReject > 0) {
    161181        for (int i = 0; i < nImages; i++) {
    162             char rejfile[MAXCHAR];              // Filename of rejection image
     182            char rejfile[MAXCHAR];      // Filename of rejection image
    163183            sprintf(rejfile,"%s.shiftrej",config->inputs->data[i]);
    164184            psImageWriteSection((psImage*)((*rejected)->data[i]),0,0,0,NULL,0,rejfile);
    165185        }
    166186    }
     187    // Write chi^2 image
     188    iteration++;
     189    char chifile[MAXCHAR];              // Filename of chi^2 image
     190    sprintf(chifile,"chi2_%d.fits",iteration);
     191    psImageWriteSection(chi2,0,0,0,NULL,0,chifile);
     192    psFree(chi2);
    167193#endif
    168194

  • trunk/stac/src/stacConfig.c

    r2500 r2661  
    99{
    1010    fprintf (stderr, "STAC: Simultaneous Telescope Array Combination\n"
    11              "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-k FRAC] OUT IN1 IN2...\n"
     11             "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
    1212             "where\n"
    1313             "\t-h           Help (this info)\n"
    … …  
    1919             "\t-b BAD       Bad level (0)\n"
    2020             "\t-k REJ       Rejection level (k-sigma; 3.5)\n"
    21              "\t-f FRAC      Fraction of pixel to be marked before considered bad (0.8)\n"
     21             "\t-f FRAC      Fraction of pixel to be marked before considered bad (0.5)\n"
     22             "\t-d GRAD      Gradient threshold for pixel to be marked (0.0)\n"
    2223             "\tOUT          Output image\n"
    2324             "\tIN1, IN2...  Input images, which have associated .map files.\n",
    … …  
    3738    config->inputs = NULL;
    3839    config->output = NULL;
    39     config->outnx = 512;
    40     config->outny = 512;
     40    config->outnx = 1024;
     41    config->outny = 1024;
    4142    config->saturated = 65536.0;
    4243    config->bad = 0.0;
    43     config->reject = 2.75;
     44    config->reject = 3.5;
    4445    config->frac = 0.5;
     46    config->grad = 0.4;
    4547
    4648    return config;
    … …  
    7274
    7375    /* Parse command-line arguments using getopt */
    74     while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:f:")) != -1) {
     76    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:f:G:")) != -1) {
    7577        switch (opt) {
    7678          case 'h':
    … …  
    121123            }
    122124            break;
     125          case 'G':
     126            if (sscanf(optarg, "%f", &config->grad) != 1) {
     127                help(programName);
     128            }
     129            break;
    123130          default:
    124131            help(programName);

  • trunk/stac/src/stacRejection.c

    r2500 r2661  
    22#include "pslib.h"
    33#include "stac.h"
     4
     5
     6#define MAX(x,y) ((x) > (y) ? (x) : (y))
     7#define MIN(x,y) ((x) < (y) ? (x) : (y))
     8
     9float stacGradient(psImage *image,      // Input for which to measure the gradient
     10                   int x, int y         // Coordinates at which to measure the gradient
     11    )
     12{
     13    float sum = 0.0;                    // The sum of surrounding pixels
     14    int num = 0;
     15    psVector *pixels = psVectorAlloc(9, PS_TYPE_F32); // Array of pixels
     16    psVector *mask = psVectorAlloc(9, PS_TYPE_U8); // Corresponding mask
     17
     18    // Get limits
     19    int xMin = MAX(x - 1, 0);
     20    int xMax = MIN(x + 1, image->numCols - 1);
     21    int yMin = MAX(y - 1, 0);
     22    int yMax = MIN(y + 1, image->numRows - 1);
     23    for (int j = yMin; j <= yMax; j++) {
     24        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    }
     43    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     44    (void)psVectorStats(stats, pixels, mask, 0);
     45    float median = stats->sampleMedian;
     46    psFree(stats);
     47    psFree(pixels);
     48    psFree(mask);
     49    return median / image->data.F32[y][x];
     50}
     51
     52#if 0
     53float stacGradient(psImage *image,      // Input for which to measure the gradient
     54                   int x, int y         // Coordinates at which to measure the gradient
     55                   )
     56{
     57    float sum = 0.0;                    // The sum of surrounding pixels
     58    float maxDiff = 0.0;
     59    int num = 0;
     60    // Get limits
     61    int xMin = MAX(x - 1, 0);
     62    int xMax = MIN(x + 1, image->numCols - 1);
     63    int yMin = MAX(y - 1, 0);
     64    int yMax = MIN(y + 1, image->numRows - 1);
     65    for (int j = yMin; j <= yMax; j++) {
     66        for (int i = xMin; i <= xMax; i++) {
     67            if (image->data.F32[j][i] - image->data.F32[y][x] < maxDiff) {
     68                maxDiff = image->data.F32[j][i] - image->data.F32[y][x];
     69            }
     70            sum += image->data.F32[j][i];
     71            num++;
     72        }
     73    }
     74    sum -= image->data.F32[y][x];
     75    sum /= (float)(num-1);
     76    sum -= image->data.F32[y][x];
     77//    return sum / image->data.F32[y][x];
     78    return maxDiff / image->data.F32[y][x];
     79}   
     80#endif
    481
    582
    … …  
    55132        int nxInput = ((psImage*)(inputs->data[i]))->numCols;
    56133        int nyInput = ((psImage*)(inputs->data[i]))->numRows;
    57         psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The mask in the source frame
     134        psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_U8); // The pixel mask for the input
     135#ifdef TESTING
     136        psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
     137        psImage *grad = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);    // The gradient image
     138#endif
    58139        psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
    59140        psPlaneTransform *map = maps->data[i]; // The map from input to output
    60141
    61142        psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
     143        int nBad = 0;                   // Number of bad pixels
     144
     145#ifdef CRFLUX
     146        // Set up CR output
     147        FILE *crs = NULL;               // File for outputting details of rejected pixels
     148        char crfile[MAXCHAR];   // Filename
     149        sprintf(crfile,"%s.cr",config->inputs->data[i]);
     150        if ((crs = fopen(crfile, "w")) == NULL) {
     151            psError("stac.rejection","Unable to open file for detailed output, %s\n");
     152            return NULL;
     153        }
     154#endif
    62155
    63156        // Transform the mask
    … …  
    67160        for (int y = 0; y < nyInput; y++) {
    68161            for (int x = 0; x < nxInput; x++) {
    69                 inCoords->x = (double)x + 0.5;
    70                 inCoords->y = (double)y + 0.5;
     162                inCoords->x = (double)x;
     163                inCoords->y = (double)y;
    71164                (void)psPlaneTransformApply(outCoords, map, inCoords);
    72                 mask->data.F32[y][x] = psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
     165                float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
    73166                                                               NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     167#ifdef TESTING
     168                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++;
     175#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));
     178#endif
     179                } else {
     180                    mask->data.U8[y][x] = 0;
     181                }
    74182
    75183            }
    76         }
     184        } // Iterating over pixels
     185
     186#ifdef CRFLUX
     187        // Close file used for detailed output
     188        fclose(crs);
     189#endif
     190
     191        psTrace("stac.rejection", 1, "%d pixels masked in image %d.\n", nBad, i);
     192        // Clip the image, and convert to suitable mask format
    77193
    78194#ifdef TESTING
    79195        // Write error image out to check
    80         char maskfile[MAXCHAR]; // Filename of mask image
     196        char maskfile[MAXCHAR];         // Filename of mask image
     197        char rejmapfile[MAXCHAR];       // Filename of rejection image
     198        char gradfile[MAXCHAR];         // Filename of gradient image
    81199        sprintf(maskfile,"%s.mask",config->inputs->data[i]);
     200        sprintf(rejmapfile,"%s.rejmap",config->inputs->data[i]);
     201        sprintf(gradfile,"%s.grad",config->inputs->data[i]);
    82202        psImageWriteSection(mask,0,0,0,NULL,0,maskfile);
    83 #endif
    84 
    85         // Clip the image, and convert to suitable mask format
    86         (void)psImageClip(mask, config->frac, 0.0, config->frac, 1.0);
    87         psImage *maskU8 = psImageCopy(NULL, mask, PS_TYPE_U8);
    88         psFree(mask);
    89    
     203        psImageWriteSection(rejmap,0,0,0,NULL,0,rejmapfile);
     204        psImageWriteSection(grad,0,0,0,NULL,0,gradfile);
     205        psFree(rejmap);
     206        psFree(grad);
     207#endif
     208
    90209        // Stuff into the array
    91         inputRej->data[i] = maskU8;
    92 
    93     }
     210        inputRej->data[i] = mask;
     211
     212    }
     213
    94214
    95215    psFree(inCoords);

  • trunk/stac/src/stacTransform.c

    r2500 r2661  
    2020                                                   unsigned int maskVal,
    2121                                                   psF64 unexposedValue)
    22 { 
     22{
    2323    double floorX = floor((psF64)(x) - 0.5);
    2424    double floorY = floor((psF64)(y) - 0.5);
    … …  
    206206                (void)psPlaneTransformApply(detector, map, sky);
    207207
    208                 outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x, detector->y,
    209                                                                           mask, 1, 0.0,
     208                outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x + 0.5,
     209                                                                          detector->y + 0.5, mask, 1, 0.0,
    210210                                                                          PS_INTERPOLATE_BILINEAR);
    211                 outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, detector->x,
    212                                                                                         detector->y, mask, 1,
    213                                                                                         0.0);
     211                outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
     212                                                                                        detector->x + 0.5,
     213                                                                                        detector->y + 0.5,
     214                                                                                        mask, 1, 0.0);
    214215                outError->data.F32[y][x] = sqrtf(outError->data.F32[y][x]);
    215216            }
Note: See TracChangeset for help on using the changeset viewer.