IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 13, 2010, 2:34:04 PM (16 years ago)
Author:
Paul Price
Message:

It's bothered me for a while that the noise model that we have been using is only kind of spatially variable. The variance map is intended to track the change in the noise properties as a function of position, and while it does keep track of different noise levels due to, e.g., different detectors, vignetting, etc., it has not tracked the different noise that comes from convolution kernels changing over an image. This had been an intentional oversight --- the prescription we were using for calculating the variance and covariance removed the spatial variation (by prescribing that the variance level under a convolution remain unchanged) and we were sweeping it under the carpet so we could calculate the covariance pseudo-matrix and use that for setting the absolute level of the noise.

I've now figured out how to have our cake and eat it too. By changing where the scalings are applied, the variance map can track the changing noise introduced by a varying convolution kernel, while the covariance pseudo-matrix tracks an absolute change. This fix is purely algorithmic --- it does not affect any of the variances/covariances stored on disk, nor is it affected by them. I hope that this means that any Magic recertification required is either minimal or completely unnecessary.

Again, all that has changed is where the scalings are applied.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/archive/noise_model/simulate.c

    r28173 r28667  
    66#define DET_SIZE 2000
    77#define SKY_SIZE 1000
    8 #define SIZE 1024
    98#define OUTPUT_ROOT "test"
    10 #define SCALE 0.654321
     9#define SCALE 0.7654321
    1110#define ROT M_PI_2
    12 #define INTERPOLATION PS_INTERPOLATE_LANCZOS4
     11#define INTERPOLATION PS_INTERPOLATE_LANCZOS3
    1312#define OFFSET 16
    1413#define SMOOTH_SIGMA 6.54321
    … …  
    1615#define DUAL_KERNEL "sub.subkernel"
    1716#define WARP_NUM 10000
     17#define CONV_NUM 1000
    1818
    1919static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 };
    … …  
    2121static const char *rootNames[] = { "o5298g0209o.fake.warp", "o5298g0210o.fake.warp", "o5298g0211o.fake.warp",
    2222                                   "o5298g0212o.fake.warp", "o5298g0213o.fake.warp", "o5298g0214o.fake.warp",
     23
    2324                                   "o5298g0215o.fake.warp", "o5298g0216o.fake.warp" };
     25#if 1
    2426static const char *kernels[] = { "test.5.kernel", "test.11.kernel", "test.17.kernel", "test.23.kernel",
    2527                                 "test.29.kernel", "test.35.kernel", "test.41.kernel", "test.47.kernel" };
    26 
     28#else
     29static const char *kernels[] = { "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel",
     30                                 "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel" };
     31#endif
    2732
    2833void writeImage(const psImage *image, const char *suffix)
    … …  
    7681float meanVar(const psImage *var, const psImage *mask, const psKernel *covar)
    7782{
    78     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    7983    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    80     psImageBackground(stats, NULL, var, mask, 0xFF, rng);
     84    psImageStats(stats, var, mask, 0xFF);
    8185    float variance = stats->sampleMedian * psImageCovarianceFactor(covar);
    8286    psFree(stats);
    83     psFree(rng);
    8487    return variance;
    8588}
    … …  
    9396        for (int x = 0; x < image->numCols; x++) {
    9497            sn->data.F32[y][x] = image->data.F32[y][x] / sqrtf(var->data.F32[y][x] * varFactor);
     98            if (!isfinite(sn->data.F32[y][x])) {
     99                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0xFF;
     100            }
    95101        }
    96102    }
    … …  
    100106    float offset = stats->sampleMean;
    101107    psFree(stats);
     108
     109    static int number = 0;
     110
     111    psString snName = NULL;
     112    psStringAppend(&snName, "sn.%d.fits", number);
     113    writeImage(sn, snName);
     114    fprintf(stderr, "Writing S/N image %d\n", number);
    102115
    103116    psHistogram *hist = psHistogramAlloc(-5, +5, 101);
    … …  
    121134    psFree(sn);
    122135
    123     static int number = 0;
    124136    psString name = NULL;
    125137    fprintf(stderr, "Writing histogram %d\n", number);
    126     psStringAppend(&name, "hist_%d.dat", number++);
     138    psStringAppend(&name, "hist_%d.dat", number);
    127139    FILE *file = fopen(name, "w");
    128140    psFree(name);
    … …  
    136148    psFree(hist);
    137149
     150    number++;
     151
    138152    return noise;
    139153}
    … …  
    144158{
    145159    psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
    146     double sum = 0.0, sum2 = 0.0;
    147160    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    148161        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    149             sum += kernel->kernel[y][x];
    150             sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
    151         }
    152     }
    153     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    154         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    155             kernel2->kernel[y][x] /= sum2;
     162            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
    156163        }
    157164    }
    … …  
    176183
    177184    psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
     185    double sum2 = 0.0;                                               // Sum of kernel squared
     186    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
     187        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
     188            sum2 += PS_SQR(kernel->kernel[y][x]);
     189        }
     190    }
     191    float factor = 1.0 / sum2;
    178192    psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar);
    179193    psFree(kernel);
    180     psImageCovarianceTransfer(smoothVariance, smoothCovar);
     194
     195    // Apply square root of significance image scaling factor to image
     196    psBinaryOp(smoothImage, smoothImage, "*", psScalarAlloc(sqrtf(factor), PS_TYPE_F32));
     197
    181198#else
    182199    psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
    183200    psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
    184     double sum = 0.0, sum2 = 0.0;
    185201    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    186202        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    187             sum += kernel->kernel[y][x];
    188             sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
    189         }
    190     }
    191     fprintf(stderr, "Kernel sum: %f %f\n", sum, sum2);
    192     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    193         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    194             kernel2->kernel[y][x] /= sum2;
     203            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
    195204        }
    196205    }
    … …  
    332341}
    333342
     343float covarianceSum(const psKernel *covar)
     344{
     345    if (!covar) {
     346        return 1.0;
     347    }
     348    double sum = 0.0;
     349    for (int y = covar->yMin; y <= covar->yMax; y++) {
     350        for (int x = covar->xMin; x <= covar->xMax; x++) {
     351            sum += covar->kernel[y][x];
     352        }
     353    }
     354    return sum;
     355}
     356
     357
    334358int main(int argc, char *argv[])
    335359{
    … …  
    383407        }
    384408
    385         fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n",
     409        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
    386410                i,
    387411                signoise(inImage, inMask, inVariance, inCovar),
    388412                psImageCovarianceFactor(inCovar),
    389                 meanVar(inVariance, inMask, inCovar));
     413                meanVar(inVariance, inMask, inCovar),
     414                covarianceSum(inCovar));
    390415
    391416        phot(inImage, inMask, inVariance, inCovar);
    … …  
    417442        }
    418443
    419         fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n",
     444        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
    420445                i,
    421446                signoise(warpImage, warpMask, warpVariance, warpCovar),
    422447                psImageCovarianceFactor(warpCovar),
    423                 meanVar(warpVariance, warpMask, warpCovar));
     448                meanVar(warpVariance, warpMask, warpCovar),
     449                covarianceSum(warpCovar));
    424450
    425451        phot(warpImage, warpMask, warpVariance, warpCovar);
    … …  
    434460        pmReadoutReadSubtractionKernels(ro, fits);
    435461        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
     462        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels);
     463        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels);
     464#if 1
     465//        kernels->solution1->data.F64[normIndex] += 1.0;
     466        kernels->solution1->data.F64[bgIndex] = 0.0;
     467#endif
     468        fprintf(stderr, "Norm: %f BG: %f\n", kernels->solution1->data.F64[normIndex], kernels->solution1->data.F64[bgIndex]);
    436469#if 0
    437470        for (int i = 0; i < kernels->num; i++) {
    … …  
    455488        psFitsClose(fits);
    456489
     490#if 0
     491        {
     492            psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     493            psArray *covariances = psArrayAlloc(CONV_NUM);
     494            psVector *factors = psVectorAlloc(CONV_NUM, PS_TYPE_F32);
     495            double mean = 0.0;
     496            for (int i = 0; i < CONV_NUM; i++) {
     497                float x, y;
     498                p_pmSubtractionPolynomialNormCoords(&x, &y, psRandomUniform(rng) * SKY_SIZE,
     499                                                    psRandomUniform(rng) * SKY_SIZE,
     500                                                    kernels->xMin, kernels->xMax,
     501                                                    kernels->yMin, kernels->yMax);
     502                psKernel *kernel = pmSubtractionKernel(kernels, x, y, false);
     503                psKernel *covar = covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
     504                psFree(kernel);
     505                mean += factors->data.F32[i] = psImageCovarianceFactor(covar);
     506            }
     507            psFree(rng);
     508            psFree(covariances);
     509
     510            mean /= WARP_NUM;
     511
     512            double stdev = 0.0;
     513            for (int i = 0; i < CONV_NUM; i++) {
     514                stdev += PS_SQR(factors->data.F32[i] - mean);
     515            }
     516            stdev = sqrt(stdev/(CONV_NUM-1));
     517            fprintf(stderr, "Conv covariance mean: %f stdev: %f\n", mean, stdev);
     518            psFree(factors);
     519        }
     520#endif
     521
    457522        pmReadout *conv = pmReadoutAlloc(NULL);
    458523#if 1
    … …  
    460525        conv->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
    461526        conv->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
    462         if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.01,
     527        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 2 * kernels->size + 1, 0.0, 0.01,
    463528                                       0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
    464529            psErrorStackPrint(stderr, "Error:");
    … …  
    494559        }
    495560
    496         fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n",
     561        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
    497562                i,
    498563                signoise(conv->image, conv->mask, conv->variance, conv->covariance),
    499564                psImageCovarianceFactor(conv->covariance),
    500                 meanVar(conv->variance, conv->mask, conv->covariance));
     565                meanVar(conv->variance, conv->mask, conv->covariance),
     566                covarianceSum(conv->covariance));
    501567
    502568        phot(conv->image, conv->mask, conv->variance, conv->covariance);
    503569        readouts->data[i] = conv;
    504570
     571        //        exit(1);
    505572    }
    506573
    … …  
    624691                meanVar(diffVariance, diffMask, diffCovar));
    625692
    626         phot(diffImage, diffMask, diffVariance, diffCovar);
    627 
    628693        writeImage(diffImage, "ssdiff.image.fits");
    629694        writeImage(diffMask, "ssdiff.mask.fits");
Note: See TracChangeset for help on using the changeset viewer.