IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 10, 2010, 6:28:51 PM (16 years ago)
Author:
watersc1
Message:

Skycell Summary and stack Association stuff should be finished. I'll merge and do final tests on monday.

Location:
branches/czw_branch/20100519
Files:
3 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20100519

  • branches/czw_branch/20100519/archive/noise_model

    • Property svn:ignore set to
      *.fits
      hist_*.dat
      *.ps

  • branches/czw_branch/20100519/archive/noise_model/simulate.c

    r28006 r28304  
    1010#define SCALE 0.654321
    1111#define ROT M_PI_2
    12 #define INTERPOLATION PS_INTERPOLATE_LANCZOS3
     12#define INTERPOLATION PS_INTERPOLATE_LANCZOS4
    1313#define OFFSET 16
    1414#define SMOOTH_SIGMA 6.54321
    1515#define SMOOTH_N_SIGMA 2.0
    1616#define DUAL_KERNEL "sub.subkernel"
     17#define WARP_NUM 10000
    1718
    1819static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 };
    … …  
    9495        }
    9596    }
    96     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    97     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV);
    98     psImageBackground(stats, NULL, sn, mask, 0xFF, rng);
    99     float noise = stats->robustStdev;
     97    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     98    psImageStats(stats, sn, mask, 0xFF);
     99    float noise = stats->sampleStdev;
     100    float offset = stats->sampleMean;
    100101    psFree(stats);
    101     psFree(rng);
     102
     103    psHistogram *hist = psHistogramAlloc(-5, +5, 101);
     104    psVector *data = psVectorAlloc(sn->numCols * sn->numRows, PS_TYPE_F32);
     105    psVector *dataMask = psVectorAlloc(sn->numCols * sn->numRows, PS_TYPE_VECTOR_MASK);
     106    psVectorInit(dataMask, 0);
     107    long num = 0;
     108    for (int y = 0, i = 0; y < sn->numRows; y++) {
     109        for (int x = 0; x < sn->numCols; x++, i++) {
     110            data->data.F32[i] = sn->data.F32[y][x] - offset;
     111            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) {
     112                dataMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF;
     113            } else {
     114                num++;
     115            }
     116        }
     117    }
     118    psVectorHistogram(hist, data, NULL, dataMask, 0xFF);
     119    psFree(data);
     120    psFree(dataMask);
    102121    psFree(sn);
     122
     123    static int number = 0;
     124    psString name = NULL;
     125    fprintf(stderr, "Writing histogram %d\n", number);
     126    psStringAppend(&name, "hist_%d.dat", number++);
     127    FILE *file = fopen(name, "w");
     128    psFree(name);
     129    fprintf(file, "# Sig Frac\n");
     130    for (int i = 0; i < hist->bounds->n - 1; i++) {
     131        fprintf(file, "%f %f\n",
     132                0.5 * (hist->bounds->data.F32[i] + hist->bounds->data.F32[i+1]),
     133                hist->nums->data.F32[i] / num);
     134    }
     135    fclose(file);
     136    psFree(hist);
     137
    103138    return noise;
    104139}
    105140
     141void convolve(psImage **smoothImage, psImage **smoothMask, psImage **smoothVariance, psKernel **smoothCovar,
     142              const psImage *image, const psImage *mask, const psImage *variance, const psKernel *covar,
     143              const psKernel *kernel)
     144{
     145    psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
     146    double sum = 0.0, sum2 = 0.0;
     147    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
     148        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;
     156        }
     157    }
     158    *smoothImage = psImageConvolveDirect(NULL, image, kernel);
     159    *smoothVariance = psImageConvolveDirect(NULL, variance, kernel2);
     160    *smoothCovar = psImageCovarianceCalculate(kernel, covar);
     161    *smoothMask = psImageCopy(NULL, mask, PS_TYPE_IMAGE_MASK); // Cheating
     162    psFree(kernel2);
     163}
     164
     165
    106166void phot(psImage *image, psImage *mask, psImage *variance, psKernel *covar)
    107167{
     168#if 1
    108169    psImage *smoothImage = psImageCopy(NULL, image, PS_TYPE_F32);
     170    psImage *smoothVariance = psImageCopy(NULL, variance, PS_TYPE_F32);
    109171    psImageSmoothMask(smoothImage, smoothImage, mask, 0xFF, SMOOTH_SIGMA, SMOOTH_N_SIGMA, 0.1);
    110     psImage *smoothVariance = psImageCopy(NULL, variance, PS_TYPE_F32);
    111172    psImageSmoothMask(smoothVariance, smoothVariance, mask, 0xFF,
    112                       SMOOTH_SIGMA * M_SQRT1_2, SMOOTH_N_SIGMA, 0.1);
     173                      SMOOTH_SIGMA * M_SQRT1_2, SMOOTH_N_SIGMA / M_SQRT1_2, 0.1);
    113174    int extent = SMOOTH_SIGMA * SMOOTH_N_SIGMA + 0.5;
    114175    psImage *smoothMask = psImageConvolveMask(NULL, mask, 0xFF, 0xFF, -extent, extent, -extent, extent);
    … …  
    118179    psFree(kernel);
    119180    psImageCovarianceTransfer(smoothVariance, smoothCovar);
     181#else
     182    psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
     183    psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
     184    double sum = 0.0, sum2 = 0.0;
     185    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
     186        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;
     195        }
     196    }
     197    psImage *smoothImage = psImageConvolveDirect(NULL, image, kernel);
     198    psImage *smoothVariance = psImageConvolveDirect(NULL, variance, kernel2);
     199    psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar);
     200    psImage *smoothMask = psImageCopy(NULL, mask, PS_TYPE_IMAGE_MASK);
     201#endif
    120202
    121203    fprintf(stderr, "Phot: S/N: %f Covar: %f Var: %f\n",
    … …  
    200282    float xOffset = psRandomUniform(rng) * 2 * OFFSET - OFFSET;
    201283    float yOffset = psRandomUniform(rng) * 2 * OFFSET - OFFSET;
    202     for (int y = 0; y < SKY_SIZE; y++) {
    203         for (int x = 0; x < SKY_SIZE; x++) {
    204             float xIn = SCALE * ((x - SKY_SIZE/2) * cos(ROT) + (y - SKY_SIZE/2) * sin(ROT)) + DET_SIZE / 2 + xOffset;
    205             float yIn = SCALE * ((x - SKY_SIZE/2) * -sin(ROT) + (y - SKY_SIZE/2) * cos(ROT)) + DET_SIZE / 2 + yOffset;
    206 
    207             double img, var;
    208             psImageMaskType msk;
     284    for (int y = 0, index = 0; y < SKY_SIZE; y++) {
     285        float dy = y - SKY_SIZE/2;
     286        for (int x = 0; x < SKY_SIZE; x++, index++) {
     287            float dx = x - SKY_SIZE/2;
     288            float xIn = SCALE * (dx * cos(ROT) + dy * sin(ROT)) + DET_SIZE / 2 + xOffset;
     289            float yIn = SCALE * (dx * -sin(ROT) + dy * cos(ROT)) + DET_SIZE / 2 + yOffset;
     290
     291#if 0
     292            xIn += 0.12345e-6 * PS_SQR(dx);
     293            yIn += 0.12345e-6 * PS_SQR(dy);
     294#endif
     295
     296            double img = NAN, var = NAN;
     297            psImageMaskType msk = 0;
    209298            psImageInterpolate(&img, &var, &msk, xIn, yIn, interp);
    210299
    211300            (*outImage)->data.F32[y][x] = img;
    212301            (*outVariance)->data.F32[y][x] = var;
    213             (*outMask)->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = isfinite(img) || isfinite(var) ? msk : 0xFF;
    214         }
    215     }
    216 
    217     psArray *covariances = psArrayAlloc(1000);
    218     for (int i = 0; i < covariances->n; i++) {
     302            (*outMask)->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = isfinite(img) && isfinite(var) ? msk : 0xFF;
     303        }
     304    }
     305
     306    psArray *covariances = psArrayAlloc(WARP_NUM);
     307    psVector *factors = psVectorAlloc(WARP_NUM, PS_TYPE_F32);
     308    double mean = 0.0;
     309    for (int i = 0; i < WARP_NUM; i++) {
    219310        psKernel *kernel = psImageInterpolationKernel(psRandomUniform(rng), psRandomUniform(rng),
    220311                                                      INTERPOLATION);
    221         covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
     312        psKernel *covar = covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
    222313        psFree(kernel);
     314        mean += factors->data.F32[i] = psImageCovarianceFactor(covar);
    223315    }
    224316    psFree(rng);
    225317    psKernel *avgCovar = psImageCovarianceAverage(covariances);
    226318    psFree(covariances);
     319
    227320    *outCovar = psImageCovarianceScale(avgCovar, SCALE);
    228321    psFree(avgCovar);
     322
     323    mean /= WARP_NUM;
     324
     325    double stdev = 0.0;
     326    for (int i = 0; i < WARP_NUM; i++) {
     327        stdev += PS_SQR(factors->data.F32[i] - mean);
     328    }
     329    stdev = sqrt(stdev/(WARP_NUM-1));
     330    fprintf(stderr, "Warp covariance mean: %f stdev: %f\n", mean, stdev);
     331    psFree(factors);
    229332}
    230333
    … …  
    265368#endif
    266369
    267         psImageCovarianceTransfer(inVariance, inCovar);
     370        //        psImageCovarianceTransfer(inVariance, inCovar);
     371
     372        {
     373            psString imageName = NULL, maskName = NULL, varName = NULL;
     374            psStringAppend(&imageName, "input.image.%d.fits", i);
     375            psStringAppend(&maskName, "input.mask.%d.fits", i);
     376            psStringAppend(&varName, "input.var.%d.fits", i);
     377            writeImage(inImage, imageName);
     378            writeImage(inMask, maskName);
     379            writeImage(inVariance, varName);
     380            psFree(imageName);
     381            psFree(maskName);
     382            psFree(varName);
     383        }
    268384
    269385        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n",
    … …  
    272388                psImageCovarianceFactor(inCovar),
    273389                meanVar(inVariance, inMask, inCovar));
     390
     391        phot(inImage, inMask, inVariance, inCovar);
    274392
    275393        psImage *warpImage = NULL, *warpMask = NULL, *warpVariance = NULL;
    … …  
    281399        psFree(inCovar);
    282400
    283         psImageCovarianceTransfer(warpVariance, warpCovar);
     401        //        psImageCovarianceTransfer(warpVariance, warpCovar);
     402
     403        {
     404            psString imageName = NULL, maskName = NULL, varName = NULL, covarName = NULL;
     405            psStringAppend(&imageName, "warp.image.%d.fits", i);
     406            psStringAppend(&maskName, "warp.mask.%d.fits", i);
     407            psStringAppend(&varName, "warp.var.%d.fits", i);
     408            psStringAppend(&covarName, "warp.covar.%d.fits", i);
     409            writeImage(warpImage, imageName);
     410            writeImage(warpMask, maskName);
     411            writeImage(warpVariance, varName);
     412            writeImage(warpCovar->image, covarName);
     413            psFree(imageName);
     414            psFree(maskName);
     415            psFree(varName);
     416            psFree(covarName);
     417        }
    284418
    285419        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n",
    … …  
    300434        pmReadoutReadSubtractionKernels(ro, fits);
    301435        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
     436#if 0
     437        for (int i = 0; i < kernels->num; i++) {
     438            for (int y = 0, index = i; y < kernels->spatialOrder; y++) {
     439                for (int x = 0; x < kernels->spatialOrder - y; x++, index++) {
     440                    if (x != 0 && y != 0) {
     441                        if (kernels->solution1) {
     442                            kernels->solution1->data.F64[index] = 0.0;
     443                        }
     444                        if (kernels->solution2) {
     445                            kernels->solution2->data.F64[index] = 0.0;
     446                        }
     447                    }
     448                }
     449            }
     450        }
     451#endif
     452
    302453        kernels->xMax = SKY_SIZE;
    303454        kernels->yMax = SKY_SIZE;
    … …  
    305456
    306457        pmReadout *conv = pmReadoutAlloc(NULL);
     458#if 1
    307459        conv->image = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
    308460        conv->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
    309461        conv->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
    310         if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.001,
     462        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.01,
    311463                                       0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
    312464            psErrorStackPrint(stderr, "Error:");
    313465            exit(1);
    314466        }
     467#else
     468        {
     469            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false);
     470            convolve(&conv->image, &conv->mask, &conv->variance, &conv->covariance,
     471                     ro->image, ro->mask, ro->variance, ro->covariance,
     472                     kernel);
     473            psFree(kernel);
     474        }
     475#endif
    315476        psFree(ro);
    316477
    317         psImageCovarianceTransfer(conv->variance, conv->covariance);
     478        //        psImageCovarianceTransfer(conv->variance, conv->covariance);
     479
     480        {
     481            psString imageName = NULL, maskName = NULL, varName = NULL, covarName = NULL;
     482            psStringAppend(&imageName, "conv.image.%d.fits", i);
     483            psStringAppend(&maskName, "conv.mask.%d.fits", i);
     484            psStringAppend(&varName, "conv.var.%d.fits", i);
     485            psStringAppend(&covarName, "conv.covar.%d.fits", i);
     486            writeImage(conv->image, imageName);
     487            writeImage(conv->mask, maskName);
     488            writeImage(conv->variance, varName);
     489            writeImage(conv->covariance->image, covarName);
     490            psFree(imageName);
     491            psFree(maskName);
     492            psFree(varName);
     493            psFree(covarName);
     494        }
    318495
    319496        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n",
    … …  
    325502        phot(conv->image, conv->mask, conv->variance, conv->covariance);
    326503        readouts->data[i] = conv;
     504
    327505    }
    328506
    … …  
    340518                psImageCovarianceFactor(diffCovar),
    341519                meanVar(diffVariance, diffMask, diffCovar));
     520
     521        phot(diffImage, diffMask, diffVariance, diffCovar);
    342522
    343523        writeImage(diffImage, "wwdiff.image.fits");
    … …  
    444624                meanVar(diffVariance, diffMask, diffCovar));
    445625
     626        phot(diffImage, diffMask, diffVariance, diffCovar);
     627
    446628        writeImage(diffImage, "ssdiff.image.fits");
    447629        writeImage(diffMask, "ssdiff.mask.fits");
    … …  
    536718                meanVar(diffVariance, diffMask, diffCovar));
    537719
     720        phot(diffImage, diffMask, diffVariance, diffCovar);
     721
    538722        writeImage(diffImage, "wsdiff.image.fits");
    539723        writeImage(diffMask, "wsdiff.mask.fits");
    … …  
    541725        writeImage(diffCovar->image, "wsdiff.covar.fits");
    542726
    543         phot(diffImage, diffMask, diffVariance, diffCovar);
    544 
    545727        psFree(diffImage);
    546728        psFree(diffMask);
Note: See TracChangeset for help on using the changeset viewer.