Index: branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionMatch.c (revision 36069) +++ branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionMatch.c (revision 36075) @@ -540,5 +540,5 @@ // Bail here if we're doing the simple matching if (type == PM_SUBTRACTION_KERNEL_SIMPLE) { - if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor)) { + if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor,optThreshold)) { return false; } Index: branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.c =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.c (revision 36069) +++ branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.c (revision 36075) @@ -63,4 +63,39 @@ } +bool simple_apply_mask(psImage *image, psImage *weight, psImage *mask, + psImageMaskType maskVal) { + for (int y = 0; y < mask->numRows; y++) { + for (int x = 0; x < mask->numCols; x++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) { + image->data.F32[y][x] = NAN; + if (weight) { + weight->data.F32[y][x] = NAN; + } + } + } + } + return(true); +} + + +// Copied from pmSubtraction +static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated + const pmSubtractionKernels *kernels, // Kernel basis functions + float value, // Normalisation value for basis function + int index // Index of basis function of interest + ) +{ + int size = kernels->size; // Kernel half-size + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + kernel->kernel[v][u] += value * preCalc->kernel->kernel[v][u]; + } + } + + return; +} +//End copy + bool pmSubtractionSimpleMatch(pmReadout *conv1, pmReadout *conv2, @@ -71,29 +106,32 @@ psImageMaskType maskVal, psImageMaskType maskBad, - psImageMaskType maskPoor + psImageMaskType maskPoor, + float deconvolveThreshold ) { // // We've already validated the input values at this level float sig2fwhm = 2.0 * sqrt(2.0 * log(2.0)); - float fwhm1,fwhm2; - float sigma1,sigma2,sigmaKern; + float fwhm1 = 0,fwhm2 = 0; + float sigma1 = 0,sigma2 = 0,sigmaKern = 0; float chisq = 1.0; int convolution_direction = 0; - psImage *image1; - psImage *mask1; - psImage *var1; - - psImage *image2; - psImage *mask2; - psImage *var2; - - psImage *imageC1; - psImage *maskC1; - psImage *varC1; - - psImage *imageC2; - psImage *maskC2; - psImage *varC2; - + psImage *image1 = NULL; + psImage *mask1 = NULL; + psImage *var1 = NULL; + + psImage *image2 = NULL; + psImage *mask2 = NULL; + psImage *var2 = NULL; + + psImage *imageC1 = NULL; + psImage *maskC1 = NULL; + psImage *varC1 = NULL; + + psImage *imageC2 = NULL; + psImage *maskC2 = NULL; + psImage *varC2 = NULL; + + psImage *maskTemp = NULL; + // Allocate images, as this is usually done by subtractionMatchAlloc after this function is called. int numCols = ro1->image->numCols; @@ -101,5 +139,5 @@ if (conv1) { - conv1->covariance = psMemIncrRefCounter(ro1->covariance); + // conv1->covariance = psMemIncrRefCounter(ro1->covariance); if (!conv1->image) { conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); @@ -118,5 +156,5 @@ } if (conv2) { - conv2->covariance = psMemIncrRefCounter(ro2->covariance); + // if (!conv2->image) { conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); @@ -171,13 +209,17 @@ if (!conv1) { if (convolution_direction == 1) { - chisq = 100; - } - convolution_direction = 2; + if (sigma1 - sigma2 > deconvolveThreshold) { + chisq = 100; + } + } + // convolution_direction = 2; } if (!conv2) { if (convolution_direction == 2) { - chisq = 100; - } - convolution_direction = 1; + if (sigma2 - sigma1 > deconvolveThreshold) { + chisq = 100; + } + } + // convolution_direction = 1; } @@ -186,114 +228,9 @@ int maskBox = (int) ceil(sigmaKern * 1.1774); // diameter is 1/2 FWHM int maskBlank = 8; // I should be able to get this from a reference, right? - - // - // Construct required kernel. No longer needed as we can direct convolve - // psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity. - // psFree(kernelVec); - - // - // Do convolutions - if (convolution_direction == 1) { - psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6); - psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6); - maskC1 = psImageConvolveMask(maskC1,mask1,maskVal,maskBad, - -maskBox,maskBox,-maskBox,maskBox); - if (conv2) { - imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32); - varC2 = psImageCopy(varC2,var2,PS_TYPE_F32); - maskC2 = psImageCopy(maskC2,mask2,PS_TYPE_IMAGE_MASK); - } - pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank); - pmSubtractionMaskApply(imageC1,varC1,maskC1,PM_SUBTRACTION_MODE_1); - } - else if (convolution_direction == 2) { - psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6); - psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6); - maskC2 = psImageConvolveMask(maskC2,mask2,maskVal,maskBad, - -maskBox,maskBox,-maskBox,maskBox); - if (conv1) { - imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32); - varC1 = psImageCopy(varC1,var1,PS_TYPE_F32); - maskC1 = psImageCopy(maskC1,mask1,PS_TYPE_IMAGE_MASK); - } - pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank); - pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2); - } - - // - // Do normalization - float normalization = 1.0; - - // Scan source list, do box photometry on peaks, and then solve the linear relation. - int photRadius = (int) floor(PS_MAX(sigma1,sigma2) * 2.0 * sqrt(2.0 * log(2.0))); // Go out a FWHM diameter from the center. - psVector *logFluxDifferences = psVectorAlloc(sources->n,PS_TYPE_F32); - psVector *fitMask = psVectorAlloc(sources->n,PS_TYPE_VECTOR_MASK); - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - int nPix1,nPix2; - float flux1,flux2; - - if (convolution_direction == 1) { - simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius); - if (conv2) { - simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius); - } - else { - simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius); - } - } - else if (convolution_direction == 2) { - simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius); - if (conv1) { - simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius); - } - else { - simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius); - } - } - logFluxDifferences->data.F32[i] = flux2 - flux1; - fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; - if ((PS_MIN(nPix1,nPix2) <= 0.75 * PS_MAX(nPix1,nPix2))|| - (!isfinite(flux1))||(!isfinite(flux2))) { - fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; - } - - // fprintf(stderr,"SOURCES: %d %g %g %g -> %d %d %g %g %d %g\n",i,source->peak->xf,source->peak->yf,source->psfMag, - // nPix1,nPix2,flux1,flux2,fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i],logFluxDifferences->data.F32[i]); - - } - - // Given the differences in log-flux space, the normalization factor is just the exponential of the median difference - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - if (!psVectorStats(stats,logFluxDifferences,NULL,fitMask,0xff)) { - // This should complain. - normalization = 1.0; - } - - normalization = pow(10,stats->robustMedian); - // fprintf(stderr,"NORM: %g+/-%g\n",stats->robustMedian,stats->robustStdev); - - psFree(stats); - psFree(logFluxDifferences); - psFree(fitMask); - - // Apply normalization - if (normalization != 1.0) { - if ((conv1)&&((convolution_direction == 1))) { - psBinaryOp(imageC1,imageC1,"*",psScalarAlloc((float) normalization, PS_TYPE_F32)); - psBinaryOp(varC1,varC1,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32)); - } - else if ((conv2)&&(convolution_direction == 2)) { - normalization = 1.0 / normalization; // Because we fit one way, but are using it in the other. - psBinaryOp(imageC2,imageC2,"*",psScalarAlloc((float) normalization, PS_TYPE_F32)); - psBinaryOp(varC2,varC2,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32)); - } - } - - + int maskPoorVal = 16384; // Another value that should be found elsewhere. // // Make a fake pmSubtractionKernels element so we can add it appropriately. - // I call it fake because we've successfully done everything at this point - // without having to define these things. + // Defining everything here is a bit clunky, but it's necessary to get the covariance + // correct. psVector *fwhms = psVectorAlloc(1,PS_TYPE_F32); fwhms->data.F32[0] = sigmaKern * sig2fwhm; @@ -329,4 +266,126 @@ kernels->numStamps = sources->n; + psKernel *kernel = psKernelAlloc(-size,size,-size,size); + solvedKernelPreCalc(kernel,kernels,1.0,0); + + // + // Do convolutions + if (convolution_direction == 1) { + if (conv1) { + psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6); + psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6); + + maskTemp = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + maskTemp = psImageConvolveMask(maskTemp,mask1,maskVal,maskBad, // Mask bad values + -maskBox,maskBox,-maskBox,maskBox); + maskC1 = psImageConvolveMask(maskC1,maskTemp,maskPoorVal,maskPoor, // Mask poor values + -maskBox,maskBox,-maskBox,maskBox); + psFree(maskTemp); + + conv1->covariance = psImageCovarianceCalculate(kernel,ro1->covariance); + pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank); + simple_apply_mask(imageC1,varC1,maskC1,maskBad); + } + if (conv2) { + imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32); + varC2 = psImageCopy(varC2,var2,PS_TYPE_F32); + maskC2 = psImageCopy(maskC2,mask2,PS_TYPE_IMAGE_MASK); + conv2->covariance = psMemIncrRefCounter(ro2->covariance); + } + } + else if (convolution_direction == 2) { + if (conv2) { + psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6); + psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6); + + maskTemp = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + maskTemp = psImageConvolveMask(maskTemp,mask2,maskVal,maskBad, // Mask bad values + -maskBox,maskBox,-maskBox,maskBox); + maskC2 = psImageConvolveMask(maskC2,maskTemp,maskPoorVal,maskPoor, // Mask poor values + -maskBox,maskBox,-maskBox,maskBox); + psFree(maskTemp); + + conv2->covariance = psImageCovarianceCalculate(kernel,ro2->covariance); + pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank); + simple_apply_mask(imageC2,varC2,maskC2,maskBad); + } + if (conv1) { + imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32); + varC1 = psImageCopy(varC1,var1,PS_TYPE_F32); + maskC1 = psImageCopy(maskC1,mask1,PS_TYPE_IMAGE_MASK); + conv1->covariance = psMemIncrRefCounter(ro1->covariance); + } + } + + psFree(kernel); // No longer needed after doing covariance calculation + + // + // Do normalization + float normalization = 1.0; + + // Scan source list, do box photometry on peaks, and then solve the linear relation. + int photRadius = (int) floor(PS_MAX(sigma1,sigma2) * 2.0 * sqrt(2.0 * log(2.0))); // Go out a FWHM diameter from the center. + psVector *logFluxDifferences = psVectorAlloc(sources->n,PS_TYPE_F32); + psVector *fitMask = psVectorAlloc(sources->n,PS_TYPE_VECTOR_MASK); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + int nPix1,nPix2; + float flux1,flux2; + + if (conv1) { + simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius); + } + else { + simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius); + } + + if (conv2) { + simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius); + } + else { + simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius); + } + + logFluxDifferences->data.F32[i] = flux2 - flux1; + fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; + if ((PS_MIN(nPix1,nPix2) <= 0.75 * PS_MAX(nPix1,nPix2))|| + (!isfinite(flux1))||(!isfinite(flux2))) { + fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; + } + + // fprintf(stderr,"SOURCES: %d %g %g %g -> %d %d %g %g %d %g\n",i,source->peak->xf,source->peak->yf,source->psfMag, + // nPix1,nPix2,flux1,flux2,fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i],logFluxDifferences->data.F32[i]); + } + + // Given the differences in log-flux space, the normalization factor is just the exponential of the median difference + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + if (!psVectorStats(stats,logFluxDifferences,NULL,fitMask,0xff)) { + // This should complain. + normalization = 1.0; + } + + normalization = pow(10,stats->robustMedian); + // fprintf(stderr,"NORM: %g+/-%g\n",stats->robustMedian,stats->robustStdev); + + psFree(stats); + psFree(logFluxDifferences); + psFree(fitMask); + + // Apply normalization + if (normalization != 1.0) { + if ((conv1)&&((convolution_direction == 1))) { + psBinaryOp(imageC1,imageC1,"*",psScalarAlloc((float) normalization, PS_TYPE_F32)); + psBinaryOp(varC1,varC1,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32)); + } + else if ((conv2)&&(convolution_direction == 2)) { + normalization = 1.0 / normalization; // Because we fit one way, but are using it in the other. + psBinaryOp(imageC2,imageC2,"*",psScalarAlloc((float) normalization, PS_TYPE_F32)); + psBinaryOp(varC2,varC2,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32)); + } + } + + + // + // // Actually add it to the headers Index: branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.h =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.h (revision 36069) +++ branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.h (revision 36075) @@ -16,5 +16,6 @@ psImageMaskType maskVal, psImageMaskType maskBad, - psImageMaskType maskPoor + psImageMaskType maskPoor, + float deconvolveThreshold ); Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/objects/pmPSFtryFitPSF.c (revision 36069) +++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmPSFtryFitPSF.c (revision 36075) @@ -75,5 +75,5 @@ psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y); - return false; + continue; }