IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 3, 2013, 4:21:05 PM (13 years ago)
Author:
watersc1
Message:

Working 1D-Gaussian convolution code. Recipes/reductions that use this have _1DG added.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtractionSimple.c

    r35726 r35771  
    2222#include "pmSubtractionVisual.h"
    2323#include "pmErrorCodes.h"
     24#include "pmSpan.h"
     25#include "pmFootprintSpans.h"
     26#include "pmFootprint.h"
     27#include "pmPeaks.h"
     28#include "pmMoments.h"
     29#include "pmModelFuncs.h"
     30
     31#include "pmSourceMasks.h"
     32#include "pmSourceExtendedPars.h"
     33#include "pmSourceDiffStats.h"
     34#include "pmSourceSatstar.h"
     35
     36#include "pmSource.h"
    2437
    2538#include "pmSubtractionSimple.h"
     39
     40
     41bool simple_do_boxphot(int *nPix,
     42                       float *flux,
     43                       pmSource *source,
     44                       psImage *image,
     45                       psImage *mask,
     46                       psImageMaskType maskVal,
     47                       int size) {
     48  *nPix = 0;
     49  *flux = 0.0;
     50  for (int y = source->peak->yf - size;y <= source->peak->yf + size;y++) {
     51    for (int x = source->peak->xf - size; x <= source->peak->xf + size; x++) {
     52      if ((y > 0)&&(y < image->numRows)&&
     53          (x > 0)&&(x < image->numCols)) {
     54        if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
     55          *nPix += 1;
     56          *flux += image->data.F32[y][x];
     57        }
     58      }
     59    }
     60  }
     61  *flux = log10(*flux);
     62  return(true);
     63} 
    2664
    2765bool pmSubtractionSimpleMatch(pmReadout *conv1,
    … …  
    4078  float fwhm1,fwhm2;
    4179  float sigma1,sigma2,sigmaKern;
     80  float chisq = 1.0;
    4281  int convolution_direction = 0;
    4382  psImage *image1;
    … …  
    62101
    63102  if (conv1) {
     103    conv1->covariance = psMemIncrRefCounter(ro1->covariance);
    64104    if (!conv1->image) {
    65105      conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    … …  
    78118  }
    79119  if (conv2) {
     120    conv2->covariance = psMemIncrRefCounter(ro2->covariance);
    80121    if (!conv2->image) {
    81122      conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    … …  
    129170  }
    130171  if (!conv1) {
     172    if (convolution_direction == 1) {
     173      chisq = 100;
     174    }
    131175    convolution_direction = 2;
    132176  }
    133177  if (!conv2) {
     178    if (convolution_direction == 2) {
     179      chisq = 100;
     180    }
    134181    convolution_direction = 1;
    135182  }
    136 
    137   //
    138   // Determine Normalization scaling
    139   psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity.
    140 
    141  
     183 
     184
     185  int maskBox = (int) ceil(sigmaKern * 1.1774); // diameter is 1/2 FWHM
     186  int maskBlank = 8;  // I should be able to get this from a reference, right?
     187
     188  //
     189  // Construct required kernel.  No longer needed as we can direct convolve
     190  //  psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity.
     191  //  psFree(kernelVec);
     192
    142193  //
    143194  // Do convolutions
    144   psImage *temp = psImageAlloc(numCols,numRows,PS_TYPE_F32);
    145   psImage *tempMask = psImageAlloc(numCols,numRows,PS_TYPE_IMAGE_MASK);
    146   psImage *tempVar  = psImageAlloc(numCols,numRows,PS_TYPE_F32);
    147 
    148   for (int y = 0; y < numRows; y++) {
    149     for (int x = 0; x < numCols; x++) {
    150       temp->data.F32[y][x] = 0.0;
    151       tempVar->data.F32[y][x] = 0.0;
    152       tempMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
    153 
    154       // Copy the image we're not convolving into the convolved output
    155       if (convolution_direction == 1) {
    156         if (conv1) {
    157           imageC1->data.F32[y][x] = 0.0;
    158           maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
    159           varC1->data.F32[y][x] = 0.0;
    160         }
    161         if (conv2) {
    162           imageC2->data.F32[y][x] = image2->data.F32[y][x];
    163           maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask2->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
    164           varC2->data.F32[y][x] = var2->data.F32[y][x];
    165         }
    166       }
    167       else if (convolution_direction == 2) {
    168         if (conv2) {
    169           imageC2->data.F32[y][x] = 0.0;
    170           maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
    171           varC2->data.F32[y][x] = 0.0;
    172         }
    173         if (conv1) {
    174           imageC1->data.F32[y][x] = image1->data.F32[y][x];
    175           maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask1->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
    176           varC1->data.F32[y][x] = var1->data.F32[y][x];
    177         }
    178       }
    179 
    180       // Do y-direction convolution
    181       for (int v = -size; v <= size; v++) {
    182         if ((y-v >= 0)&&(y-v < numRows)) {
    183           if (convolution_direction == 1) {
    184             // Insert mask/finite check here.
    185             temp->data.F32[y][x] += kernelVec->data.F32[v + size] * image1->data.F32[y - v][x];
    186             tempVar->data.F32[y][x] += kernelVec->data.F32[v + size] * var1->data.F32[y - v][x];
    187           }
    188           else if (convolution_direction == 2) {
    189             // And here.
    190             temp->data.F32[y][x] += kernelVec->data.F32[v + size] * image2->data.F32[y - v][x];
    191             tempVar->data.F32[y][x] += kernelVec->data.F32[v + size] * var2->data.F32[y - v][x];
    192           }
    193         }
    194       }
    195     }
    196   }
    197 
    198   // Do x-direction convolution
    199   for (int y = 0; y < numRows; y++) {
    200     for (int x = 0; x < numCols; x++) {
    201         for (int u = -size; u <= size; u++) {
    202           if ((x-u >= 0)&&(x-u < numCols)) {
    203             if (convolution_direction == 1) {
    204               // And here
    205               imageC1->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
    206               varC1->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
    207             }
    208             else if (convolution_direction == 2) {
    209               // And here
    210               imageC2->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
    211               varC2->data.F32[y][x] += kernelVec->data.F32[u + size] * temp->data.F32[y][x - u];
    212             }
    213           }
    214         }
    215     }
    216   }
    217   psFree(temp);
    218   psFree(kernelVec);
     195  if (convolution_direction == 1) {
     196    psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6);
     197    psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
     198    maskC1 = psImageConvolveMask(maskC1,mask1,maskVal,maskBad,
     199                                 -maskBox,maskBox,-maskBox,maskBox);
     200    if (conv2) {
     201      imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32);
     202      varC2   = psImageCopy(varC2,var2,PS_TYPE_F32);
     203      maskC2  = psImageCopy(maskC2,mask2,PS_TYPE_IMAGE_MASK);
     204    }
     205    pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank);
     206    pmSubtractionMaskApply(imageC1,varC1,maskC1,PM_SUBTRACTION_MODE_1);
     207  }
     208  else if (convolution_direction == 2) {
     209    psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6);
     210    psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
     211    maskC2 = psImageConvolveMask(maskC2,mask2,maskVal,maskBad,
     212                                 -maskBox,maskBox,-maskBox,maskBox);
     213    if (conv1) {
     214      imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32);
     215      varC1   = psImageCopy(varC1,var1,PS_TYPE_F32);
     216      maskC1  = psImageCopy(maskC1,mask1,PS_TYPE_IMAGE_MASK);
     217    }
     218    pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
     219    pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2);
     220  }   
    219221
    220222  //
    … …  
    222224  float normalization = 1.0;
    223225
    224   // Something with the source list here
     226  // Scan source list, do box photometry on peaks, and then solve the linear relation.
     227  int photRadius = (int) floor(PS_MAX(sigma1,sigma2) * 2.0 * sqrt(2.0 * log(2.0))); // Go out a FWHM diameter from the center.
     228  psVector *logFluxDifferences = psVectorAlloc(sources->n,PS_TYPE_F32);
     229  psVector *fitMask = psVectorAlloc(sources->n,PS_TYPE_VECTOR_MASK);
     230  for (int i = 0; i < sources->n; i++) {
     231    pmSource *source = sources->data[i];
     232    int nPix1,nPix2;
     233    float flux1,flux2;
     234
     235    if (convolution_direction == 1) {
     236      simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
     237      if (conv2) {
     238        simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
     239      }
     240      else {
     241        simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius);
     242      }
     243    }
     244    else if (convolution_direction == 2) {
     245      simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
     246      if (conv1) {
     247        simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
     248      }
     249      else {
     250        simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius);
     251      }
     252    }
     253    logFluxDifferences->data.F32[i] = flux2 - flux1;
     254    fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
     255    if ((PS_MIN(nPix1,nPix2) <= 0.75 * PS_MAX(nPix1,nPix2))||
     256        (!isfinite(flux1))||(!isfinite(flux2))) {
     257      fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
     258    }
     259
     260    //    fprintf(stderr,"SOURCES: %d %g %g %g -> %d %d %g %g %d %g\n",i,source->peak->xf,source->peak->yf,source->psfMag,
     261    //      nPix1,nPix2,flux1,flux2,fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i],logFluxDifferences->data.F32[i]);
     262   
     263  }
     264
     265  // Given the differences in log-flux space, the normalization factor is just the exponential of the median difference
     266  psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     267  if (!psVectorStats(stats,logFluxDifferences,NULL,fitMask,0xff)) {
     268    // This should complain.
     269    normalization = 1.0;
     270  }
     271
     272  normalization = pow(10,stats->robustMedian);
     273  // fprintf(stderr,"NORM: %g+/-%g\n",stats->robustMedian,stats->robustStdev);
     274 
     275  psFree(stats);
     276  psFree(logFluxDifferences);
     277  psFree(fitMask);
    225278
    226279  // Apply normalization
    227280  if (normalization != 1.0) {
    228     for (int y = 0; y < numRows; y++) {
    229       for (int x = 0; x < numCols; x++) {
    230         if ((conv1)&&((convolution_direction == 1))) {
    231           imageC1->data.F32[y][x] *= normalization;
    232         }
    233         else if ((conv2)&&(convolution_direction == 2)) {
    234           imageC2->data.F32[y][x] *= normalization;
    235         }
    236       }
     281    if ((conv1)&&((convolution_direction == 1))) {
     282      psBinaryOp(imageC1,imageC1,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
     283      psBinaryOp(varC1,varC1,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
     284    }
     285    else if ((conv2)&&(convolution_direction == 2)) {
     286      normalization = 1.0 / normalization; // Because we fit one way, but are using it in the other.
     287      psBinaryOp(imageC2,imageC2,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
     288      psBinaryOp(varC2,varC2,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
    237289    }
    238290  }
    … …  
    263315    psFree(kernels->preCalc->data[0]);
    264316  }
    265   kernels->preCalc->data[0] = preCalc;
     317  kernels->preCalc->data[0] = psMemIncrRefCounter(preCalc);
    266318  kernels->solution1 = psVectorAlloc(3,PS_TYPE_F64);
    267319  kernels->solution1->data.F32[0] = 1.0;
    … …  
    272324  kernels->solution1err->data.F32[1] = 0.0;
    273325  kernels->solution1err->data.F32[2] = 0.0;
     326  kernels->mean = 0.0;
     327  kernels->rms = chisq; // This is the chi^2 value that's passed to ppStack
     328  kernels->numStamps = sources->n;
    274329 
    275330  //
    … …  
    284339    psFree(kernels);
    285340  }
     341  // This is a hack.  Yes, I know.  pmSAnalysis doesn't get the normalization correct, so I just do it here.
     342  psMetadataRemoveKey(outAnalysis,PM_SUBTRACTION_ANALYSIS_NORM);
     343  psMetadataRemoveKey(outHeader,PM_SUBTRACTION_ANALYSIS_NORM);
     344  psMetadataAddF32(outAnalysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_REPLACE, "Normalisation", normalization);
     345  psMetadataAddF32(outHeader,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_REPLACE, "Normalisation", normalization);
     346
    286347  psFree(fwhms);
    287348  psFree(orders);
    … …  
    291352  if (conv1) {
    292353    conv1->analysis = psMetadataCopy(conv1->analysis, outAnalysis);
     354    conv1->data_exists = true;
     355    conv1->parent->data_exists = true;
     356    conv1->parent->parent->data_exists = true;
    293357  }
    294358  if (conv2) {
    295359    conv2->analysis = psMetadataCopy(conv2->analysis, outAnalysis);
     360    conv2->data_exists = true;
     361    conv2->parent->data_exists = true;
     362    conv2->parent->parent->data_exists = true;
    296363  }
    297364
Note: See TracChangeset for help on using the changeset viewer.