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

DM-23781: Improve Sky Object Placement #312

Merged
merged 1 commit into from
Apr 13, 2023
Merged

DM-23781: Improve Sky Object Placement #312

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.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 18 additions & 10 deletions python/lsst/meas/algorithms/skyObjects.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

__all__ = ["SkyObjectsConfig", "SkyObjectsTask", "generateSkyObjects"]

from scipy.stats import qmc

from lsst.pex.config import Config, Field, ListField
from lsst.pipe.base import Task

Expand Down Expand Up @@ -33,9 +35,11 @@ def generateSkyObjects(mask, seed, config):
through the provided `mask` (in which objects are typically flagged
as `DETECTED`).

The algorithm for determining sky objects is random trial and error:
we try up to `nTrialSkySources` random positions to find `nSources`
sky objects.
Sky objects are positioned using a quasi-random Halton sequence number
generator. This is a deterministic sequence that mimics a random trial and
error approach whilst acting to minimize clustering of points for a given
field of view. Up to `nTrialSources` points are generated, returning the
first `nSources` that do not overlap with the mask.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

RE:

"acting to minimize clustering of points"

does this algorithm guarantee no overlapping positions? If not, you might consider adding:

# Add doubled-in-size sky object spanSet to the avoid mask.
avoid = avoid.union(spans.dilated(int(skySourceRadius)))

after the line:

skyFootprints.append(fp)

at line 93 below in order to have a cumulative avoid mask. It worked nicely using the old trial 'n error algorithm (but may not be necessary now?)

Copy link
Contributor Author

@leeskelvin leeskelvin Apr 12, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's nothing explicit in the SciPy Halton Sequence implementation that prevents overlapping positions as N tends to infinity, however, it is (by design) exceptionally unlikely.

However, I do like your additional logic to prevent this and I don't think it adds too much overhead. As SkyObjectsConfig already has a growMask field, I'm tempted to recommend using this information rather than hard-coding a double-in-size condition?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Following testing, I think I convinced myself that I prefer your initial idea of growing the sky object footprints by their own radius. I experimented with growing the growMask parameter, and while it works, it occurred to me that that parameter was intended to mask the flux in the wings of astrophysical sources. We want to grow sky object footprints for a different reason altogether (to avoid clustering for stats purposes). To that end, I've adopted your initial recommendation above as-is.


Parameters
----------
Expand Down Expand Up @@ -71,15 +75,14 @@ def generateSkyObjects(mask, seed, config):
if config.growMask > 0:
avoid = avoid.dilated(config.growMask)

rng = lsst.afw.math.Random(seed=seed)
sampler = qmc.Halton(d=2, seed=seed).random(nTrialSkySources)
sample = qmc.scale(sampler, [xMin, yMin], [xMax, yMax])

skyFootprints = []
for _ in range(nTrialSkySources):
for x, y in zip(sample[:, 0].astype(int), sample[:, 1].astype(int)):
if len(skyFootprints) == nSkySources:
break

x = int(rng.flat(xMin, xMax))
y = int(rng.flat(yMin, yMax))
spans = lsst.afw.geom.SpanSet.fromShape(int(skySourceRadius), offset=(x, y))
if spans.overlaps(avoid):
continue
Expand All @@ -88,6 +91,9 @@ def generateSkyObjects(mask, seed, config):
fp.addPeak(x, y, 0)
skyFootprints.append(fp)

# Add doubled-in-size sky object spanSet to the avoid mask.
avoid = avoid.union(spans.dilated(int(skySourceRadius)))

return skyFootprints


Expand All @@ -103,9 +109,11 @@ def run(self, mask, seed):
through the provided `mask` (in which objects are typically flagged
as `DETECTED`).

The algorithm for determining sky objects is random trial and error:
we try up to `nTrialSkySources` random positions to find `nSources`
sky objects.
Sky objects are positioned using a quasi-random Halton sequence
number generator. This is a deterministic sequence that mimics a random
trial and error approach whilst acting to minimize clustering of points
for a given field of view. Up to `nTrialSources` points are generated,
returning the first `nSources` that do not overlap with the mask.

Parameters
----------
Expand Down
7 changes: 6 additions & 1 deletion tests/test_dynamicDetection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,12 @@ def setUp(self):

# Relative tolerance for tweak factor
# Not sure why this isn't smaller; maybe due to use of Poisson instead of Gaussian noise?
self.rtol = 0.1
# It seems as if some sky objects are being placed in the extra
# background region, which is causing the offset between the expected
# factor and the measured factor to be larger than otherwise expected.
# This relative tolerance was increased from 0.1 to 0.15 on DM-23781 to
# account for this.
self.rtol = 0.15
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be worth a comment in the commit message about this change (not that we can really explain it at the moment...😞).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought about adding one, but the initial in-line that's already there seems to capture everything I would want to say anyway. In the spirit of adding something, I've added this:

# It seems as if some sky objects are being placed in the extra
# background region, which is causing the offset between the expected
# factor and the measured factor to be larger than otherwise expected.
# This relative tolerance was increased from 0.1 to 0.15 on DM-23781 to
# account for this.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good (I did just mean adding something to the commit message, but this works too!)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, the commit message! For some reason I didn't parse that the first time around. Commit message also updated with some background too.


def tearDown(self):
del self.exposure
Expand Down