Index: trunk/archive/noise_model/simulate.c =================================================================== --- trunk/archive/noise_model/simulate.c (revision 28006) +++ trunk/archive/noise_model/simulate.c (revision 28159) @@ -10,9 +10,10 @@ #define SCALE 0.654321 #define ROT M_PI_2 -#define INTERPOLATION PS_INTERPOLATE_LANCZOS3 +#define INTERPOLATION PS_INTERPOLATE_LANCZOS4 #define OFFSET 16 #define SMOOTH_SIGMA 6.54321 #define SMOOTH_N_SIGMA 2.0 #define DUAL_KERNEL "sub.subkernel" +#define WARP_NUM 10000 static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 }; @@ -94,21 +95,80 @@ } } - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV); - psImageBackground(stats, NULL, sn, mask, 0xFF, rng); - float noise = stats->robustStdev; + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_STDEV); + psImageStats(stats, sn, mask, 0xFF); + float noise = stats->sampleStdev; psFree(stats); - psFree(rng); + + psHistogram *hist = psHistogramAlloc(-5, +5, 1001); + psVector *data = psVectorAlloc(sn->numCols * sn->numRows, PS_TYPE_F32); + psVector *dataMask = psVectorAlloc(sn->numCols * sn->numRows, PS_TYPE_VECTOR_MASK); + psVectorInit(dataMask, 0); + long num = 0; + for (int y = 0, i = 0; y < sn->numRows; y++) { + for (int x = 0; x < sn->numCols; x++, i++) { + data->data.F32[i] = sn->data.F32[y][x]; + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) { + dataMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF; + } else { + num++; + } + } + } + psVectorHistogram(hist, data, NULL, dataMask, 0xFF); + psFree(data); + psFree(dataMask); psFree(sn); + + static int number = 0; + psString name = NULL; + fprintf(stderr, "Writing histogram %d\n", number); + psStringAppend(&name, "hist_%d.dat", number++); + FILE *file = fopen(name, "w"); + psFree(name); + fprintf(stderr, "# Sig Frac\n"); + for (int i = 0; i < hist->bounds->n - 1; i++) { + fprintf(file, "%f %f\n", + 0.5 * (hist->bounds->data.F32[i] + hist->bounds->data.F32[i+1]), + hist->nums->data.F32[i] / num); + } + fclose(file); + psFree(hist); + return noise; } +void convolve(psImage **smoothImage, psImage **smoothMask, psImage **smoothVariance, psKernel **smoothCovar, + const psImage *image, const psImage *mask, const psImage *variance, const psKernel *covar, + const psKernel *kernel) +{ + psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax); + double sum = 0.0, sum2 = 0.0; + for (int y = kernel->yMin; y <= kernel->yMax; y++) { + for (int x = kernel->xMin; x <= kernel->xMax; x++) { + sum += kernel->kernel[y][x]; + sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]); + } + } + for (int y = kernel->yMin; y <= kernel->yMax; y++) { + for (int x = kernel->xMin; x <= kernel->xMax; x++) { + kernel2->kernel[y][x] /= sum2; + } + } + *smoothImage = psImageConvolveDirect(NULL, image, kernel); + *smoothVariance = psImageConvolveDirect(NULL, variance, kernel2); + *smoothCovar = psImageCovarianceCalculate(kernel, covar); + *smoothMask = psImageCopy(NULL, mask, PS_TYPE_IMAGE_MASK); // Cheating + psFree(kernel2); +} + + void phot(psImage *image, psImage *mask, psImage *variance, psKernel *covar) { +#if 1 psImage *smoothImage = psImageCopy(NULL, image, PS_TYPE_F32); + psImage *smoothVariance = psImageCopy(NULL, variance, PS_TYPE_F32); psImageSmoothMask(smoothImage, smoothImage, mask, 0xFF, SMOOTH_SIGMA, SMOOTH_N_SIGMA, 0.1); - psImage *smoothVariance = psImageCopy(NULL, variance, PS_TYPE_F32); psImageSmoothMask(smoothVariance, smoothVariance, mask, 0xFF, - SMOOTH_SIGMA * M_SQRT1_2, SMOOTH_N_SIGMA, 0.1); + SMOOTH_SIGMA * M_SQRT1_2, SMOOTH_N_SIGMA / M_SQRT1_2, 0.1); int extent = SMOOTH_SIGMA * SMOOTH_N_SIGMA + 0.5; psImage *smoothMask = psImageConvolveMask(NULL, mask, 0xFF, 0xFF, -extent, extent, -extent, extent); @@ -118,4 +178,25 @@ psFree(kernel); psImageCovarianceTransfer(smoothVariance, smoothCovar); +#else + psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing + psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax); + double sum = 0.0, sum2 = 0.0; + for (int y = kernel->yMin; y <= kernel->yMax; y++) { + for (int x = kernel->xMin; x <= kernel->xMax; x++) { + sum += kernel->kernel[y][x]; + sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]); + } + } + fprintf(stderr, "Kernel sum: %f %f\n", sum, sum2); + for (int y = kernel->yMin; y <= kernel->yMax; y++) { + for (int x = kernel->xMin; x <= kernel->xMax; x++) { + kernel2->kernel[y][x] /= sum2; + } + } + psImage *smoothImage = psImageConvolveDirect(NULL, image, kernel); + psImage *smoothVariance = psImageConvolveDirect(NULL, variance, kernel2); + psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar); + psImage *smoothMask = psImageCopy(NULL, mask, PS_TYPE_IMAGE_MASK); +#endif fprintf(stderr, "Phot: S/N: %f Covar: %f Var: %f\n", @@ -200,31 +281,52 @@ float xOffset = psRandomUniform(rng) * 2 * OFFSET - OFFSET; float yOffset = psRandomUniform(rng) * 2 * OFFSET - OFFSET; - for (int y = 0; y < SKY_SIZE; y++) { - for (int x = 0; x < SKY_SIZE; x++) { - float xIn = SCALE * ((x - SKY_SIZE/2) * cos(ROT) + (y - SKY_SIZE/2) * sin(ROT)) + DET_SIZE / 2 + xOffset; - float yIn = SCALE * ((x - SKY_SIZE/2) * -sin(ROT) + (y - SKY_SIZE/2) * cos(ROT)) + DET_SIZE / 2 + yOffset; - - double img, var; - psImageMaskType msk; + for (int y = 0, index = 0; y < SKY_SIZE; y++) { + float dy = y - SKY_SIZE/2; + for (int x = 0; x < SKY_SIZE; x++, index++) { + float dx = x - SKY_SIZE/2; + float xIn = SCALE * (dx * cos(ROT) + dy * sin(ROT)) + DET_SIZE / 2 + xOffset; + float yIn = SCALE * (dx * -sin(ROT) + dy * cos(ROT)) + DET_SIZE / 2 + yOffset; + +#if 0 + xIn += 0.12345e-6 * PS_SQR(dx); + yIn += 0.12345e-6 * PS_SQR(dy); +#endif + + double img = NAN, var = NAN; + psImageMaskType msk = 0; psImageInterpolate(&img, &var, &msk, xIn, yIn, interp); (*outImage)->data.F32[y][x] = img; (*outVariance)->data.F32[y][x] = var; - (*outMask)->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = isfinite(img) || isfinite(var) ? msk : 0xFF; - } - } - - psArray *covariances = psArrayAlloc(1000); - for (int i = 0; i < covariances->n; i++) { + (*outMask)->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = isfinite(img) && isfinite(var) ? msk : 0xFF; + } + } + + psArray *covariances = psArrayAlloc(WARP_NUM); + psVector *factors = psVectorAlloc(WARP_NUM, PS_TYPE_F32); + double mean = 0.0; + for (int i = 0; i < WARP_NUM; i++) { psKernel *kernel = psImageInterpolationKernel(psRandomUniform(rng), psRandomUniform(rng), INTERPOLATION); - covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar); + psKernel *covar = covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar); psFree(kernel); + mean += factors->data.F32[i] = psImageCovarianceFactor(covar); } psFree(rng); psKernel *avgCovar = psImageCovarianceAverage(covariances); psFree(covariances); + *outCovar = psImageCovarianceScale(avgCovar, SCALE); psFree(avgCovar); + + mean /= WARP_NUM; + + double stdev = 0.0; + for (int i = 0; i < WARP_NUM; i++) { + stdev += PS_SQR(factors->data.F32[i] - mean); + } + stdev = sqrt(stdev/(WARP_NUM-1)); + fprintf(stderr, "Warp covariance mean: %f stdev: %f\n", mean, stdev); + psFree(factors); } @@ -265,5 +367,18 @@ #endif - psImageCovarianceTransfer(inVariance, inCovar); + // psImageCovarianceTransfer(inVariance, inCovar); + + { + psString imageName = NULL, maskName = NULL, varName = NULL; + psStringAppend(&imageName, "input.image.%d.fits", i); + psStringAppend(&maskName, "input.mask.%d.fits", i); + psStringAppend(&varName, "input.var.%d.fits", i); + writeImage(inImage, imageName); + writeImage(inMask, maskName); + writeImage(inVariance, varName); + psFree(imageName); + psFree(maskName); + psFree(varName); + } fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n", @@ -272,4 +387,6 @@ psImageCovarianceFactor(inCovar), meanVar(inVariance, inMask, inCovar)); + + phot(inImage, inMask, inVariance, inCovar); psImage *warpImage = NULL, *warpMask = NULL, *warpVariance = NULL; @@ -281,5 +398,21 @@ psFree(inCovar); - psImageCovarianceTransfer(warpVariance, warpCovar); + // psImageCovarianceTransfer(warpVariance, warpCovar); + + { + psString imageName = NULL, maskName = NULL, varName = NULL, covarName = NULL; + psStringAppend(&imageName, "warp.image.%d.fits", i); + psStringAppend(&maskName, "warp.mask.%d.fits", i); + psStringAppend(&varName, "warp.var.%d.fits", i); + psStringAppend(&covarName, "warp.covar.%d.fits", i); + writeImage(warpImage, imageName); + writeImage(warpMask, maskName); + writeImage(warpVariance, varName); + writeImage(warpCovar->image, covarName); + psFree(imageName); + psFree(maskName); + psFree(varName); + psFree(covarName); + } fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n", @@ -300,4 +433,21 @@ pmReadoutReadSubtractionKernels(ro, fits); pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); +#if 0 + for (int i = 0; i < kernels->num; i++) { + for (int y = 0, index = i; y < kernels->spatialOrder; y++) { + for (int x = 0; x < kernels->spatialOrder - y; x++, index++) { + if (x != 0 && y != 0) { + if (kernels->solution1) { + kernels->solution1->data.F64[index] = 0.0; + } + if (kernels->solution2) { + kernels->solution2->data.F64[index] = 0.0; + } + } + } + } + } +#endif + kernels->xMax = SKY_SIZE; kernels->yMax = SKY_SIZE; @@ -305,15 +455,41 @@ pmReadout *conv = pmReadoutAlloc(NULL); +#if 1 conv->image = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32); conv->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK); conv->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32); - if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.001, + if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.01, 0xFF, 0xF0, 0x0F, 0.1, 1.0)) { psErrorStackPrint(stderr, "Error:"); exit(1); } +#else + { + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); + convolve(&conv->image, &conv->mask, &conv->variance, &conv->covariance, + ro->image, ro->mask, ro->variance, ro->covariance, + kernel); + psFree(kernel); + } +#endif psFree(ro); - psImageCovarianceTransfer(conv->variance, conv->covariance); + // psImageCovarianceTransfer(conv->variance, conv->covariance); + + { + psString imageName = NULL, maskName = NULL, varName = NULL, covarName = NULL; + psStringAppend(&imageName, "conv.image.%d.fits", i); + psStringAppend(&maskName, "conv.mask.%d.fits", i); + psStringAppend(&varName, "conv.var.%d.fits", i); + psStringAppend(&covarName, "conv.covar.%d.fits", i); + writeImage(conv->image, imageName); + writeImage(conv->mask, maskName); + writeImage(conv->variance, varName); + writeImage(conv->covariance->image, covarName); + psFree(imageName); + psFree(maskName); + psFree(varName); + psFree(covarName); + } fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n", @@ -325,4 +501,5 @@ phot(conv->image, conv->mask, conv->variance, conv->covariance); readouts->data[i] = conv; + }