XMM-Newton Science Analysis System


asmooth (asmooth-2.32.1) [xmmsas_20230412_1735-21.0.0]

Algorithm

read parameters;
read --inset;

if (--withoutmaskset) {
  read --outmaskset;
} else {
  outMask = TRUE;
}

if (--withweightset) {
  read --weightset;
} else {
  weightImage = 1;
}

weightImage = weightImage / maxval(weightImage)

(xyLimits) = &findWidthOfBlackBorder(inImage);
inImage     = &cutBorder(inImage,     xyLimits);
outMask     = &cutBorder(outMask,     xyLimits);
weightImage = &cutBorder(weightImage, xyLimits);

if (--withexpimageset) {
  read --expimageset;
  expImage = &cutBorder(expImage, xyLimits);

  where(expImage>0) {
    inImage = inImage / expImage;
  } elsewhere {
    inImage = 0;
    weightImage = 0;
  }
}

# Choice of smoothing type:

if (--smoothstyle='simple') {
  if (--convolverstyle='gaussian') {
    convolvers(1) = &makeGaussianConvolver;
  } elsif(--convolverstyle='tophat') {
    convolvers(1) = &makeTopHatConvolver;
  } elsif(--convolverstyle='squarebox') {
    convolvers(1) = &makeBoxConvolver;
  }
} elsif(--smoothstyle='template') {
  (indexImage, convolvers) = &makeConvolversFromTemplate;
} elsif(--smoothstyle='withset') {
  (indexImage, convolvers) = &readConvolversFromSets;

} elsif(--smoothstyle='adaptive') {

  # Adaptive smoothing calculation of convolvers:

  if (--readvarianceset) {
    read --invarianceset;
    varianceImage = &cutBorder(varianceImage, xyLimits);
  } else {
    if (--withexpimageset) {
      varianceImage = inImage * expMapImage;
      # Because inImage has already been divided by expMapImage, and we
      # need to reverse that.
    } else {
      varianceImage = inImage;
    }
  }

  if (--withexpimageset) {
    if (--expmapuse='samesnr') {
      where(expImage>0) {
        varianceImage = varianceImage / expImage / expImage;
      } elsewhere {
        varianceImage = 0;
      }
    } else {
      where(expImage>0) {
        varianceImage = varianceImage / expImage / maxval(expImage);
      } elsewhere {
        varianceImage = 0;
      }
    }
  }

  convolvers = &calculateConvolverLibrary;

  where(outMask) {
    indexImage = 1;
  } elsewhere {
    indexImage = 0;
  }

  smoothedImage = 0.0;
  rmsSmoothedImage = 0.0;

  i = 1;
  smoothedImage    = smoothedImage         + &patchSmooth(convolvers(i)
    , weightImage, mask=(indexImage==i));
  rmsSmoothedImage = rmsSmoothedImage + sqrt(&patchSmooth(convolvers**2(i)
    , weightImage, mask=(indexImage==i)));

  where(indexImage == i) {
    where(! failureMask) {
      where(rmsSmoothedImage>0 && smoothedImage>0) {
        bestSnr = smoothedImage / rmsSmoothedImage;

        where(bestSnr > desiredSnr) {
          # This means that we are still
          # smoothing too hard, and need to try this pixel again with a
          # narrower gaussian. Therefore increment the index:
 
          indexImage = i + 1;
        }
      }
    }
  }

  foreach(i = 2, numConvolvers) {
    secondBestSnr = bestSnr;

    smoothedImage    = smoothedImage         + &patchSmooth(convolvers(i)
      , weightImage, mask=(indexImage==i));
    rmsSmoothedImage = rmsSmoothedImage + sqrt(&patchSmooth(convolvers**2(i)
      , weightImage, mask=(indexImage==i)));

    where(indexImage == i) {
      where(failureMask) {
        indexImage = i - 1; # This pixel didn't fail for the previous
        # value of index (=> next larger gaussian), otherwise it would
        # never have got here. Hence we will drop the index value for
        # this pixel back to this last known 'good' value.

      } elsewhere {
        where(rmsSmoothedImage>0 && smoothedImage>0) {
          bestSnr = smoothedImage / rmsSmoothedImage;
        } elsewhere {
          bestSnr = 0;
        }

        where(bestSnr > desiredSnr) {
          # This means that we are still
          # smoothing too hard, and need to try this pixel again with a
          # narrower gaussian. Therefore increment the index:
 
          indexImage = i + 1;
        } elsewhere { # go to the value that gave the best SNR:
          where(abs(secondBestSnr-desiredSnr)<abs(bestSnr-desiredSnr)) {
            indexImage = i - 1;
            # Otherwise, leave indexImage at i.
          }
        }
      }
    }
  }

  where(indexImage > numConvolvers) {indexImage = numConvolvers;}
}

# Do the smoothing:

smoothedImage = 0.0;
where(outMask) {smoothedImage = inImage;}
failureMask = FALSE;
foreach(i = 1, numConvolvers) {
  patchMask = outMask && (indexImage == i);

  # The actual convolution:
  foreach(xi = outStartX, outFinisX) {
    foreach(yi = outStartY, outFinisY) {
    next if (! outMask(xi, yi));

      summ = 0.0;
      weight = 0.0;
      foreach(cxi = 1, convolverXSize) {
        xxi = (cxi - 1 - halfConvolverXSize) - xi;
        foreach(cyi = 1, convolverYSize) {
          yyi = (cyi - 1 - halfConvolverYSize) - yi;
        next if (weightImage(xxi, yyi) <= 0.0);
          summ   = summ   + convolver(cxi, cyi) *     inImage(xxi, yyi);
          weight = weight + convolver(cxi, cyi) * weightImage(xxi, yyi);
        }
      }

      if (weight >= minAllowedWeight) {
        outImage(xi, yi) = summ * norm / weight;
      } else {
        outImage(xi, yi) = 0.0;
        failureMask(xi, yi) = TRUE; # indicates those pixels where the
        # weight is too small.
      }
    }
  }
}

# Outputs:

if (--writevarianceset) {
  outVarianceImage = &calculateOutVariance(inImage, weightImage, outMask
    , indexImage, convolvers);

  outVarianceImage = &restoreBlankBorder(outVarianceImage, xyLimits);
  write to --outvarianceset;
}

where(failureMask) {smoothedImage = 0;}

if (--writebadmaskset) {
  failureMask = &restoreBlankBorder(failureMask, xyLimits);
  write to --badmaskset;
}

if (--smoothstyle='adaptive') {
  if (--writetemplateset) {
    templateImage = &calculateTemplateImage;
    templateImage = &restoreBlankBorder(templateImage, xyLimits);
    write to --outtemplateset;
  }

  if (--writeconvolvers) {
    indexImage = &restoreBlankBorder(indexImage, xyLimits);
    write indexImage to --outindeximageset;
    write convolvers to --outconvolversset;
  }
}

if (--withexpimageset) {
  if (--remultiply) {
    smoothedImage = smoothedImage * expMapImage;
  }
}

smoothedImage = &restoreBlankBorder(smoothedImage, xyLimits);

write to --outset;