Skip to content

Commit

Permalink
Merge pull request #31 from lsst/tickets/DM-7376
Browse files Browse the repository at this point in the history 
Fix NaiveDipoleCentroid logic
  • Loading branch information
@djreiss
djreiss committed Aug 31, 2016
2 parents c789682 + 24a03cd commit ad0d2d3
Show file tree
Hide file tree
Showing 3 changed files with 61 additions and 20 deletions.
2 changes: 1 addition & 1 deletion include/lsst/ip/diffim/DipoleAlgorithms.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -288,7 +288,7 @@ class NaiveDipoleCentroid : public DipoleCentroidAlgorithm {
afw::image::Exposure<float> const & exposure
) const;

void mergeCentroids(afw::table::SourceRecord & source) const;
void mergeCentroids(afw::table::SourceRecord & source, double posValue, double negValue) const;

void fail(
afw::table::SourceRecord & measRecord,
Expand Down
41 changes: 35 additions & 6 deletions python/lsst/ip/diffim/dipoleFitTask.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,11 @@ class DipoleFitPluginConfig(measBase.SingleFramePluginConfig):
"""!Configuration for DipoleFitPlugin
"""

fitAllDiaSources = pexConfig.Field(
dtype=float, default=False,
doc="""Attempte dipole fit of all diaSources (otherwise just the ones consisting of overlapping
positive and negative footprints)""")

maxSeparation = pexConfig.Field(
dtype=float, default=5.,
doc="Assume dipole is not separated by more than maxSeparation * psfSigma")
Expand Down Expand Up @@ -102,10 +107,15 @@ def setDefaults(self):
"base_PixelFlags",
"base_SkyCoord",
"base_PsfFlux",
"base_GaussianCentroid",
"base_PeakLikelihoodFlux",
"base_PeakCentroid",
"base_SdssCentroid",
"base_NaiveCentroid",
"ip_diffim_NaiveDipoleCentroid",
"ip_diffim_NaiveDipoleFlux",
"ip_diffim_PsfDipoleFlux",
"ip_diffim_ClassificationDipole"
"ip_diffim_ClassificationDipole",
]

self.slots.calibFlux = None
Expand Down Expand Up @@ -538,10 +548,11 @@ def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None
cenNeg = cenPos = fpCentroid

pks = fp.getPeaks()

if len(pks) >= 1:
cenPos = pks[0].getF() # if individual (merged) peaks were detected, use those
if len(pks) >= 2:
cenNeg = pks[1].getF()
if len(pks) >= 2: # peaks are already sorted by centroid flux so take the most negative one
cenNeg = pks[-1].getF()

# For close/faint dipoles the starting locs (min/max) might be way off, let's help them a bit.
# First assume dipole is not separated by more than 5*psfSigma.
Expand Down Expand Up @@ -919,6 +930,10 @@ def _setupSchema(self, config, name, schema, metadata):
schema.join(name, "flag", "classification"), type="Flag",
doc="Flag indicating diaSource is classified as a dipole")

self.classificationAttemptedFlagKey = schema.addField(
schema.join(name, "flag", "classificationAttempted"), type="Flag",
doc="Flag indicating diaSource was attempted to be classified as a dipole")

self.flagKey = schema.addField(
schema.join(name, "flag"), type="Flag",
doc="General failure flag for dipole fit")
Expand Down Expand Up @@ -946,11 +961,21 @@ def measure(self, measRecord, exposure, posExp=None, negExp=None):
testing (@see `tests/testDipoleFitter.py`)
"""

result = None
pks = measRecord.getFootprint().getPeaks()
if len(pks) <= 1: # not a dipole for our analysis

# Check if the footprint consists of a putative dipole - else don't fit it.
if (
(len(pks) <= 1) or # one peak in the footprint - not a dipole
(len(pks) > 1 and (np.sign(pks[0].getPeakValue()) ==
np.sign(pks[-1].getPeakValue()))) # peaks are same sign - not a dipole
):
measRecord.set(self.classificationFlagKey, False)
measRecord.set(self.classificationAttemptedFlagKey, False)
self.fail(measRecord, measBase.MeasurementError('not a dipole', self.FAILURE_NOT_DIPOLE))
if not self.config.fitAllDiaSources:
return result

result = None
try:
alg = self.DipoleFitAlgorithmClass(exposure, posImage=posExp, negImage=negExp)
result, _ = alg.fitDipole(
Expand All @@ -966,6 +991,8 @@ def measure(self, measRecord, exposure, posExp=None, negExp=None):
self.fail(measRecord, measBase.MeasurementError('dipole fit failure', self.FAILURE_FIT))

if result is None:
measRecord.set(self.classificationFlagKey, False)
measRecord.set(self.classificationAttemptedFlagKey, False)
return result

self.log.log(self.log.DEBUG, "Dipole fit result: %d %s" % (measRecord.getId(), str(result)))
Expand All @@ -974,7 +1001,6 @@ def measure(self, measRecord, exposure, posExp=None, negExp=None):
self.fail(measRecord, measBase.MeasurementError('dipole fit failure', self.FAILURE_FIT))

# add chi2, coord/flux uncertainties (TBD), dipole classification

# Add the relevant values to the measRecord
measRecord[self.posFluxKey] = result.posFlux
measRecord[self.posFluxSigmaKey] = result.signalToNoise # to be changed to actual sigma!
Expand Down Expand Up @@ -1032,6 +1058,8 @@ def doClassify(self, measRecord, chi2val):
passesChi2 = significance < self.config.maxChi2DoF
allPass = allPass and passesChi2

measRecord.set(self.classificationAttemptedFlagKey, True)

if allPass: # Note cannot pass `allPass` into the `measRecord.set()` call below...?
measRecord.set(self.classificationFlagKey, True)
else:
Expand All @@ -1052,6 +1080,7 @@ def fail(self, measRecord, error=None):
if error.getFlagBit() == self.FAILURE_NOT_DIPOLE:
self.log.log(self.log.DEBUG, 'DipoleFitPlugin not run on record %d: %s' %
(measRecord.getId(), str(error)))
measRecord.set(self.classificationAttemptedFlagKey, False)
measRecord.set(self.flagKey, True)
else:
self.log.warn('DipoleFitPlugin failed on record %d' % measRecord.getId())
38 changes: 25 additions & 13 deletions src/DipoleAlgorithms.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,34 +137,46 @@ void NaiveDipoleCentroid::measure(
) const {
afw::detection::PeakCatalog const& peaks = source.getFootprint()->getPeaks();

naiveCentroid(source, exposure, peaks[0].getI(), (peaks[0].getPeakValue() >= 0 ?
getPositiveKeys() :
getNegativeKeys()));
if (peaks.size() > 1) {
naiveCentroid(source, exposure, peaks[1].getI(), (peaks[1].getPeakValue() >= 0 ?
getPositiveKeys() :
getNegativeKeys()));
int posInd = 0;
double posValue = peaks[posInd].getPeakValue(), negValue = 0;
if (posValue < 0.) { /* All peaks are negative so use the *most* negative value */
posInd = peaks.size() - 1;
posValue = peaks[posInd].getPeakValue();
}
naiveCentroid(source, exposure, peaks[posInd].getI(),
(posValue >= 0 ? getPositiveKeys() : getNegativeKeys()));

if (posValue > 0. && posInd == 0 && peaks.size() > 1) { /* See if there's also a negative peak */
int negInd = peaks.size() - 1;
negValue = peaks[negInd].getPeakValue();
if (posValue > 0. && negValue < 0.) {
naiveCentroid(source, exposure, peaks[negInd].getI(),
(negValue >= 0 ? getPositiveKeys() : getNegativeKeys()));
}
}

mergeCentroids(source);
mergeCentroids(source, posValue, negValue);

}

void NaiveDipoleCentroid::mergeCentroids(afw::table::SourceRecord & source) const {
void NaiveDipoleCentroid::mergeCentroids(afw::table::SourceRecord & source,
double posValue, double negValue) const {

float pos_x, pos_y;
float neg_x, neg_y;
double pos_x, pos_y, pos_f;
double neg_x, neg_y, neg_f;

pos_x = source.get(getPositiveKeys().getX());
pos_y = source.get(getPositiveKeys().getY());
pos_f = posValue;

neg_x = source.get(getNegativeKeys().getX());
neg_y = source.get(getNegativeKeys().getY());
neg_f = -negValue;

if(std::isfinite(pos_x) && std::isfinite(pos_y) &&
std::isfinite(neg_x) && std::isfinite(neg_y)) {
source.set(getCenterKeys().getX(), 0.5*(pos_x + neg_x));
source.set(getCenterKeys().getY(), 0.5*(pos_y + neg_y));
source.set(getCenterKeys().getX(), (pos_x * pos_f + neg_x * neg_f) / (pos_f + neg_f));
source.set(getCenterKeys().getY(), (pos_y * pos_f + neg_y * neg_f) / (pos_f + neg_f));
} else if (std::isfinite(pos_x) && std::isfinite(pos_y)) {
source.set(getCenterKeys().getX(), pos_x);
source.set(getCenterKeys().getY(), pos_y);
Expand Down

0 comments on commit ad0d2d3

Please sign in to comment.