Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix NaiveDipoleCentroid logic #31

Merged
merged 1 commit into from
Aug 31, 2016
Merged

Fix NaiveDipoleCentroid logic #31

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion include/lsst/ip/diffim/DipoleAlgorithms.h
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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