IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35771 for trunk/psModules


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.

Location:
trunk/psModules/src/imcombine
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/Makefile.am

    r30622 r35771  
    2020        pmSubtractionStamps.c   \
    2121        pmSubtractionThreads.c  \
     22        pmSubtractionSimple.c   \
    2223        pmPSFEnvelope.c         \
    2324        pmSubtractionVisual.c   \
    … …  
    4243        pmSubtractionStamps.h   \
    4344        pmSubtractionThreads.h  \
     45        pmSubtractionSimple.h   \
    4446        pmPSFEnvelope.h         \
    4547        pmSubtractionVisual.h   \

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

    r35455 r35771  
    167167              break;
    168168          }
     169        case PM_SUBTRACTION_KERNEL_SIMPLE: {
     170              solvedKernelPreCalc(kernel, kernels, 1.0, i);
     171              break;
     172        }         
     173         
    169174          case PM_SUBTRACTION_KERNEL_ISIS:
    170175          case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:

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

    r30622 r35771  
    6565    return result;
    6666}
     67
     68// Generate 1D convolution kernel for SIMPLE
     69psVector *pmSubtractionKernelSIMPLE(float sigma, // Gaussian width
     70                                    int order,   // Unused polynomial order
     71                                    int size     // Kernel half-size
     72                                    )
     73{
     74  int fullSize = 2 * size + 1;    // Full size of the kernel
     75  psVector *kernel = psVectorAlloc(fullSize, PS_TYPE_F32); // Kernel to return
     76  float expNorm = -0.5 / PS_SQR(sigma); // Normalization for exponential
     77  float norm    = 1.0 / sqrtf(2.0 * M_PI * sigma * sigma); // Correct Normalization for Gaussian
     78  if (sigma < 0.1) {
     79    kernel->data.F32[size] = 1.0;
     80    return(kernel);
     81  }
     82  for (int i = 0, x = -size; x <= size; i++,x++) {
     83    kernel->data.F32[i] = norm * expf(expNorm * PS_SQR(x));
     84  }
     85  return(kernel);
     86}
     87
    6788
    6889// Generate 1D convolution kernel for ISIS
    … …  
    771792        preCalc->poly    = NULL;
    772793        break;
     794    case PM_SUBTRACTION_KERNEL_SIMPLE:
     795      preCalc->xKernel = pmSubtractionKernelSIMPLE(sigma,uOrder,size);
     796      preCalc->yKernel = pmSubtractionKernelSIMPLE(sigma,vOrder,size);
     797      preCalc->uCoords = NULL;
     798      preCalc->vCoords = NULL;
     799      preCalc->poly    = NULL;
     800      break;
    773801      case PM_SUBTRACTION_KERNEL_RINGS:
    774802        // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
    … …  
    12301258      case PM_SUBTRACTION_KERNEL_RINGS:
    12311259        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode);
     1260    case PM_SUBTRACTION_KERNEL_SIMPLE:
     1261      return pmSubtractionKernelsISIS(size,spatialOrder,fwhms,orders,penalty,bounds,mode);
    12321262      default:
    12331263        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
    … …  
    12981328      case PM_SUBTRACTION_KERNEL_HERM:
    12991329      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
     1330    case PM_SUBTRACTION_KERNEL_SIMPLE:
    13001331        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
    13011332
    … …  
    13791410        return PM_SUBTRACTION_KERNEL_RINGS;
    13801411    }
    1381 
     1412    if (strncasecmp(type, "SIMPLE", nameLength) == 0) {
     1413      return PM_SUBTRACTION_KERNEL_SIMPLE;
     1414    }
    13821415    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", type);
    13831416    return PM_SUBTRACTION_KERNEL_NONE;
    … …  
    14361469        break;
    14371470
     1471    case PM_SUBTRACTION_KERNEL_SIMPLE:
     1472      psStringAppend(&kernels->description, "SIMPLE(%d,%s)",kernels->size,params);
     1473      break;
    14381474      default:
    14391475        psAbort("unknown kernel");

  • trunk/psModules/src/imcombine/pmSubtractionKernels.h

    r30622 r35771  
    2222    } \
    2323    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_ISIS_RADIAL) { \
     24        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
     25        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
     26        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
     27    } \
     28    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_SIMPLE) { \
    2429        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
    2530        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
    … …  
    7176}
    7277
     78// Generate 1D convolution kernel for SIMPLE
     79psVector *pmSubtractionKernelSIMPLE(float sigma, // Gaussian width
     80                                    int order,   // Unused polynomial order
     81                                    int size     // Kernel half-size
     82                                    );
     83
    7384// Generate 1D convolution kernel for ISIS
    7485psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width

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

    r33424 r35771  
    2121#include "pmSubtractionThreads.h"
    2222#include "pmSubtractionVisual.h"
     23#include "pmSubtractionSimple.h"
    2324#include "pmErrorCodes.h"
    2425
    … …  
    537538    }
    538539
     540    // Bail here if we're doing the simple matching
     541    if (type == PM_SUBTRACTION_KERNEL_SIMPLE) {
     542      if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor)) {
     543        return false;
     544      }
     545      return(true);
     546    }
     547
    539548    // Where does our variance map come from?
    540549    // Getting the variance exactly right is not necessary --- it's just used for weighting.
    … …  
    593602        goto MATCH_ERROR;
    594603    }
     604
    595605   
    596606    memCheck("start");
    … …  
    648658        psFree(bg);
    649659    }
    650 
     660   
    651661    subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, subMode, numCols, numRows);
    652 
     662   
    653663    // Iterate over iso-kernel regions
    654664    for (int j = 0; j < yRegions; j++) {

  • 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

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

    r35716 r35771  
    1313                              const pmReadout *ro2,
    1414                              const psArray *sources,
    15                               int size
     15                              int size,
     16                              psImageMaskType maskVal,
     17                              psImageMaskType maskBad,
     18                              psImageMaskType maskPoor
    1619                              );
    1720

  • trunk/psModules/src/imcombine/pmSubtractionTypes.h

    r31671 r35771  
    4242    PM_SUBTRACTION_KERNEL_GUNK,         ///< Grid United with Normal Kernel --- POIS and ISIS hybrid
    4343    PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
     44    PM_SUBTRACTION_KERNEL_SIMPLE,       ///< Simple Gaussian kernel to avoid complications
    4445} pmSubtractionKernelsType;
    4546
Note: See TracChangeset for help on using the changeset viewer.