IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 28, 2013, 4:24:43 PM (13 years ago)
Author:
watersc1
Message:

Commit before switch to apparently-already-available-convolution-code.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r35717 r35726  
    3030                              const pmReadout *ro2,
    3131                              const psArray *sources,
    32                               int size
     32                              int size,
     33                              psImageMaskType maskVal,
     34                              psImageMaskType maskBad,
     35                              psImageMaskType maskPoor
    3336                              ) {
    3437  //
    … …  
    3841  float sigma1,sigma2,sigmaKern;
    3942  int convolution_direction = 0;
     43  psImage *image1;
     44  psImage *mask1;
     45  psImage *var1;
     46
     47  psImage *image2;
     48  psImage *mask2;
     49  psImage *var2;
     50
     51  psImage *imageC1;
     52  psImage *maskC1;
     53  psImage *varC1;
     54
     55  psImage *imageC2;
     56  psImage *maskC2;
     57  psImage *varC2;
    4058
    4159  // Allocate images, as this is usually done by subtractionMatchAlloc after this function is called. 
    42   psImageMaskType maskBad = 1;
    4360  int numCols = ro1->image->numCols;
    4461  int numRows = ro1->image->numRows;
    45  
    46   if (!conv1->image) {
    47     conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    48   }
    49   psImageInit(conv1->image, NAN);
    50   if (ro1->variance) {
    51     if (!conv1->variance) {
    52       conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    53     }
    54     psImageInit(conv1->variance, NAN);
    55   }
    56   if (!conv1->mask) {
    57     conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    58   }
    59   psImageInit(conv1->mask, maskBad);
    60   if (!conv2->image) {
    61     conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    62   }
    63   psImageInit(conv2->image, NAN);
    64   if (ro2->variance) {
    65     if (!conv2->variance) {
    66       conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    67     }
    68     psImageInit(conv2->variance, NAN);
    69   }
    70   if (!conv2->mask) {
    71     conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    72   }
    73   psImageInit(conv2->mask, maskBad);
    74 
    75  
    76  
    77  
    78   psImage *image1 = ro1->image;
    79   psImage *mask1  = ro1->mask;
    80   psImage *var1   = ro1->variance;
    81 
    82   psImage *image2 = ro2->image;
    83   psImage *mask2  = ro2->mask;
    84   psImage *var2   = ro2->variance;
    85 
    86   psImage *imageC1 = conv1->image;
    87   psImage *maskC1  = conv1->mask;
    88   psImage *varC1   = conv1->variance;
    89 
    90   psImage *imageC2 = conv2->image;
    91   psImage *maskC2  = conv2->mask;
    92   psImage *varC2   = conv2->variance;
    93  
     62
     63  if (conv1) {
     64    if (!conv1->image) {
     65      conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     66    }
     67    psImageInit(conv1->image, NAN);
     68    if (ro1->variance) {
     69      if (!conv1->variance) {
     70        conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     71      }
     72      psImageInit(conv1->variance, NAN);
     73    }
     74    if (!conv1->mask) {
     75      conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     76    }
     77    psImageInit(conv1->mask, maskBad);
     78  }
     79  if (conv2) {
     80    if (!conv2->image) {
     81      conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     82    }
     83    psImageInit(conv2->image, NAN);
     84    if (ro2->variance) {
     85      if (!conv2->variance) {
     86        conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     87      }
     88      psImageInit(conv2->variance, NAN);
     89    }
     90    if (!conv2->mask) {
     91      conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     92    }
     93    psImageInit(conv2->mask, maskBad);
     94  }
     95 
     96  // Assign local aliases to images
     97  image1 = ro1->image;
     98  mask1  = ro1->mask;
     99  var1   = ro1->variance;
     100
     101  image2 = ro2->image;
     102  mask2  = ro2->mask;
     103  var2   = ro2->variance;
     104
     105  if (conv1) {
     106    imageC1 = conv1->image;
     107    maskC1  = conv1->mask;
     108    varC1   = conv1->variance;
     109  }
     110 
     111  if (conv2) {
     112    imageC2 = conv2->image;
     113    maskC2  = conv2->mask;
     114    varC2   = conv2->variance;
     115  }
     116 
    94117  //
    95118  // Determine Gaussian widths
    … …  
    105128    sigmaKern = sqrt(PS_SQR(sigma2) - PS_SQR(sigma1));
    106129  }
     130  if (!conv1) {
     131    convolution_direction = 2;
     132  }
     133  if (!conv2) {
     134    convolution_direction = 1;
     135  }
    107136
    108137  //
    … …  
    110139  psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity.
    111140
    112   float normalization = 1.0;
    113 
    114   for (int i = 0; i < 2 * size + 1; i++) {
    115     kernelVec->data.F32[i] *= sqrt(normalization); // Needs to be this because we apply this vector twice.
    116   }
    117141 
    118142  //
    119143  // Do convolutions
    120144  psImage *temp = psImageAlloc(numCols,numRows,PS_TYPE_F32);
    121   for (int y = 0; y < image1->numRows; y++) {
    122     for (int x = 0; x < image1->numCols; x++) {
    123 
     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++) {
    124150      temp->data.F32[y][x] = 0.0;
    125      
     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
    126155      if (convolution_direction == 1) {
    127         imageC1->data.F32[y][x] = 0.0;
    128         maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
    129         varC1->data.F32[y][x] = 0.0;
    130 
    131         imageC2->data.F32[y][x] = image2->data.F32[y][x];
    132         maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask2->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
    133         varC2->data.F32[y][x] = var2->data.F32[y][x];
     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        }
    134166      }
    135167      else if (convolution_direction == 2) {
    136         imageC2->data.F32[y][x] = 0.0;
    137         maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
    138         varC2->data.F32[y][x] = 0.0;
    139 
    140         imageC1->data.F32[y][x] = image1->data.F32[y][x];
    141         maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask1->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
    142         varC1->data.F32[y][x] = var1->data.F32[y][x];
    143       }
     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
    144181      for (int v = -size; v <= size; v++) {
    145182        if ((y-v >= 0)&&(y-v < numRows)) {
    146183          if (convolution_direction == 1) {
     184            // Insert mask/finite check here.
    147185            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];
    148187          }
    149188          else if (convolution_direction == 2) {
     189            // And here.
    150190            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];
    151192          }
    152193        }
    … …  
    154195    }
    155196  }
    156   for (int y = 0; y < image1->numRows; y++) {
    157     for (int x = 0; x < image1->numCols; x++) {
     197
     198  // Do x-direction convolution
     199  for (int y = 0; y < numRows; y++) {
     200    for (int x = 0; x < numCols; x++) {
    158201        for (int u = -size; u <= size; u++) {
    159202          if ((x-u >= 0)&&(x-u < numCols)) {
    160203            if (convolution_direction == 1) {
     204              // And here
    161205              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];
    162207            }
    163208            else if (convolution_direction == 2) {
     209              // And here
    164210              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];
    165212            }
    166213          }
    … …  
    171218  psFree(kernelVec);
    172219
     220  //
     221  // Do normalization
     222  float normalization = 1.0;
     223
     224  // Something with the source list here
     225
     226  // Apply normalization
     227  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      }
     237    }
     238  }
     239 
     240
     241  //
    173242  // Make a fake pmSubtractionKernels element so we can add it appropriately.
    174243  // I call it fake because we've successfully done everything at this point
    … …  
    199268  kernels->solution1->data.F32[1] = 0.0;
    200269  kernels->solution1->data.F32[2] = 0.0;
    201 /*   kernels->solution2 = psVectorAlloc(3,PS_TYPE_F64); */
    202 /*   kernels->solution2->data.F32[0] = 1.0; */
    203270  kernels->solution1err = psVectorAlloc(3,PS_TYPE_F32);
    204271  kernels->solution1err->data.F32[0] = 0.0;
    205272  kernels->solution1err->data.F32[1] = 0.0;
    206273  kernels->solution1err->data.F32[2] = 0.0;
    207 /*   kernels->solution2err = psVectorAlloc(3,PS_TYPE_F32); */
    208 /*   kernels->solution2err->data.F32[0] = 0.0; */
    209 /*   kernels->solution2err->data.F32[1] = 0.0; */
    210 /*   kernels->solution2err->data.F32[2] = 0.0; */
    211  
    212 
     274 
     275  //
    213276  // Actually add it to the headers
    214  
    215277  psMetadata *outAnalysis = psMetadataAlloc();
    216278  psMetadata *outHeader   = psMetadataAlloc();
    … …  
    248310  psFree(outAnalysis);
    249311  psFree(outHeader);
     312
     313  //
    250314  // Return
    251 
    252315  return(true);
    253316}
Note: See TracChangeset for help on using the changeset viewer.