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

DM-38703: Adding stars that are missing from bright stars stamps to Bright Star Subtracting pipeline #339

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

DM-38703: Adding stars that are missing from bright stars stamps to Bright Star Subtracting pipeline #339

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fa4c7c5
DM-38703: Adding stars that are missing from bright stars stamps to B…
bazkiaei Jun 25, 2023
5473c6e
minor correction to pass testes
bazkiaei Jul 6, 2023
5404897
some typos corrected
bazkiaei Jul 20, 2023
0dcb3da
some modifications that are due to some changes in pipe_tasks are app…
bazkiaei Sep 4, 2023
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
165 changes: 116 additions & 49 deletions python/lsst/meas/algorithms/brightStarStamps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,8 @@ class BrightStarStamp(AbstractStamp):
archive_element: Persistable | None = None
annularFlux: float | None = None
minValidAnnulusFraction: float = 0.0
optimalInnerRadius: int | None = None
optimalOuterRadius: int | None = None

@classmethod
def factory(cls, stamp_im, metadata, idx, archive_element=None, minValidAnnulusFraction=0.0):
Expand Down Expand Up @@ -143,27 +145,26 @@ def measureAndNormalize(
Collection of mask planes to ignore when computing annularFlux.
"""
stampSize = self.stamp_im.getDimensions()
# create image: same pixel values within annulus, NO_DATA elsewhere
# Create image: science pixel values within annulus, NO_DATA elsewhere
maskPlaneDict = self.stamp_im.mask.getMaskPlaneDict()
annulusImage = MaskedImageF(stampSize, planeDict=maskPlaneDict)
annulusMask = annulusImage.mask
annulusMask.array[:] = 2 ** maskPlaneDict["NO_DATA"]
annulus.copyMaskedImage(self.stamp_im, annulusImage)
# set mask planes to be ignored
# Set mask planes to be ignored.
andMask = reduce(ior, (annulusMask.getPlaneBitMask(bm) for bm in badMaskPlanes))
statsControl.setAndMask(andMask)

annulusStat = makeStatistics(annulusImage, statsFlag, statsControl)
# determining the number of valid (unmasked) pixels within the annulus
# Determine the number of valid (unmasked) pixels within the annulus.
unMasked = annulusMask.array.size - np.count_nonzero(annulusMask.array)
# should we have some messages printed here?
logger.info(
"The Star's annulus contains %s valid pixels and the annulus itself contains %s pixels.",
unMasked,
annulus.getArea(),
)
if unMasked > (annulus.getArea() * self.minValidAnnulusFraction):
# compute annularFlux
# Compute annularFlux.
self.annularFlux = annulusStat.getValue()
logger.info("Annular flux is: %s", self.annularFlux)
else:
Expand All @@ -174,7 +175,7 @@ def measureAndNormalize(
raise RuntimeError("Annular flux computation failed, likely because there are no valid pixels.")
if self.annularFlux < 0:
raise RuntimeError("The annular flux is negative. The stamp can not be normalized!")
# normalize stamps
# Normalize stamps.
self.stamp_im.image.array /= self.annularFlux
return None

Expand Down Expand Up @@ -244,8 +245,7 @@ def __init__(
use_archive=False,
):
super().__init__(starStamps, metadata, use_mask, use_variance, use_archive)
# Ensure stamps contain a flux measurement if and only if they are
# already expected to be normalized
# Ensure stamps contain a flux measure if expected to be normalized.
self._checkNormalization(False, innerRadius, outerRadius)
self._innerRadius, self._outerRadius = innerRadius, outerRadius
if innerRadius is not None and outerRadius is not None:
Expand All @@ -267,6 +267,7 @@ def initAndNormalize(
use_archive=False,
imCenter=None,
discardNanFluxObjects=True,
forceFindFlux=False,
statsControl=StatisticsControl(),
statsFlag=stringToStatisticsProperty("MEAN"),
badMaskPlanes=("BAD", "SAT", "NO_DATA"),
Expand Down Expand Up @@ -331,14 +332,17 @@ def initAndNormalize(
stamps, stamps normalized with different annulus definitions, or if
stamps are to be normalized but annular radii were not provided.
"""
stampSize = starStamps[0].stamp_im.getDimensions()
if imCenter is None:
stampSize = starStamps[0].stamp_im.getDimensions()
imCenter = stampSize[0] // 2, stampSize[1] // 2
# Create SpanSet of annulus

# Create SpanSet of annulus.
outerCircle = SpanSet.fromShape(outerRadius, Stencil.CIRCLE, offset=imCenter)
innerCircle = SpanSet.fromShape(innerRadius, Stencil.CIRCLE, offset=imCenter)
annulusWidth = outerRadius - innerRadius
annulus = outerCircle.intersectNot(innerCircle)
# Initialize (unnormalized) brightStarStamps instance

# Initialize (unnormalized) brightStarStamps instance.
bss = cls(
starStamps,
innerRadius=None,
Expand All @@ -349,37 +353,99 @@ def initAndNormalize(
use_variance=use_variance,
use_archive=use_archive,
)
# Ensure no stamps had already been normalized.

# Ensure that no stamps have already been normalized.
bss._checkNormalization(True, innerRadius, outerRadius)
bss._innerRadius, bss._outerRadius = innerRadius, outerRadius
# Create a list to contain rejected stamps.

# Apply normalization.
rejects = []
# Apply normalization
badStamps = []
for stamp in bss._stamps:
try:
stamp.measureAndNormalize(
annulus, statsControl=statsControl, statsFlag=statsFlag, badMaskPlanes=badMaskPlanes
)
# Stars that are missing from the first steps results might
# still have a flux within the normalization annulus. The
# following two lines make sure that these stars are included
# in the subtraction process. Failing to assign the optimal
# radii values may result in an error in `createAnnulus`
# method of `SubtractBrightStarsTask` class. An alternative to
# handle this, is to create to types of stams that are missing
# from the input brightStarStamps. One for those that have
# flux within the normalization annulus and another for those
# that do not have a flux within the normalization annulus.
stamp.optimalOuterRadius = outerRadius
stamp.optimalInnerRadius = innerRadius
except RuntimeError as err:
logger.error(err)
# Optionally keep NaN flux objects, for bookkeeping purposes,
# and to avoid having to re-find and redo the preprocessing
# steps needed before bright stars can be subtracted.
if discardNanFluxObjects:
rejects.append(stamp)
elif forceFindFlux:
newInnerRadius = innerRadius
newOuterRadius = outerRadius
while True:
newOuterRadius += annulusWidth
newInnerRadius += annulusWidth
if newOuterRadius > min(imCenter):
logger.info("No flux found for the star with Gaia ID of %s", stamp.gaiaId)
stamp.annularFlux = None
badStamps.append(stamp)
break
newOuterCircle = SpanSet.fromShape(newOuterRadius, Stencil.CIRCLE, offset=imCenter)
newInnerCircle = SpanSet.fromShape(newInnerRadius, Stencil.CIRCLE, offset=imCenter)
newAnnulus = newOuterCircle.intersectNot(newInnerCircle)
try:
stamp.measureAndNormalize(
newAnnulus,
statsControl=statsControl,
statsFlag=statsFlag,
badMaskPlanes=badMaskPlanes,
)

except RuntimeError:
stamp.annularFlux = np.nan
logger.error(
"The annular flux was not found for radii %d and %d",
newInnerRadius,
newOuterRadius,
)
if stamp.annularFlux and stamp.annularFlux > 0:
logger.info("The flux is found within an optimized annulus.")
logger.info(
"The optimized annulus radii are %d and %d and the flux is %f",
newInnerRadius,
newOuterRadius,
stamp.annularFlux,
)
stamp.optimalOuterRadius = newOuterRadius
stamp.optimalInnerRadius = newInnerRadius
break
else:
stamp.annularFlux = np.nan

# Remove rejected stamps.
bss.normalized = True
if discardNanFluxObjects:
for reject in rejects:
bss._stamps.remove(reject)
bss.normalized = True
elif forceFindFlux:
for badStamp in badStamps:
bss._stamps.remove(badStamp)
bss._innerRadius, bss._outerRadius = None, None
return bss, badStamps
return bss

def _refresh_metadata(self):
"""Refresh metadata. Should be called before writing the object out."""
# add full list of positions, Gaia magnitudes, IDs and annularFlxes to
# shared metadata
"""Refresh metadata. Should be called before writing the object out.

This method adds full lists of positions, Gaia magnitudes, IDs and
annular fluxes to the shared metadata.
"""
self._metadata["G_MAGS"] = self.getMagnitudes()
self._metadata["GAIA_IDS"] = self.getGaiaIds()
positions = self.getPositions()
Expand All @@ -400,7 +466,7 @@ def readFits(cls, filename):
Parameters
----------
filename : `str`
Name of the file to read
Name of the file to read.
"""
return cls.readFitsWithOptions(filename, None)

Expand All @@ -411,9 +477,9 @@ def readFitsWithOptions(cls, filename, options):
Parameters
----------
filename : `str`
Name of the file to read
Name of the file to read.
options : `PropertyList`
Collection of metadata parameters
Collection of metadata parameters.
"""
stamps, metadata = readFitsWithOptions(filename, BrightStarStamp.factory, options)
nb90Rots = metadata["NB_90_ROTS"] if "NB_90_ROTS" in metadata else None
Expand Down Expand Up @@ -481,11 +547,12 @@ def extend(self, bss):
self._stamps += bss._stamps

def getMagnitudes(self):
"""Retrieve Gaia G magnitudes for each star.
"""Retrieve Gaia G-band magnitudes for each star.

Returns
-------
gaiaGMags : `list` [`float`]
Gaia G-band magnitudes for each star.
"""
return [stamp.gaiaGMag for stamp in self._stamps]

Expand All @@ -495,20 +562,23 @@ def getGaiaIds(self):
Returns
-------
gaiaIds : `list` [`int`]
Gaia IDs for each star.
"""
return [stamp.gaiaId for stamp in self._stamps]

def getAnnularFluxes(self):
"""Retrieve normalization factors for each star.
"""Retrieve normalization factor for each star.

These are computed by integrating the flux in annulus centered on the
bright star, far enough from center to be beyond most severe ghosts and
saturation. The inner and outer radii that define the annulus can be
recovered from the metadata.
saturation.
The inner and outer radii that define the annulus can be recovered from
the metadata.

Returns
-------
annularFluxes : `list` [`float`]
Annular fluxes which give the normalization factor for each star.
"""
return [stamp.annularFlux for stamp in self._stamps]

Expand All @@ -528,8 +598,7 @@ def selectByMag(self, magMin=None, magMax=None):
for stamp in self._stamps
if (magMin is None or stamp.gaiaGMag > magMin) and (magMax is None or stamp.gaiaGMag < magMax)
]
# This is an optimization to save looping over the init argument when
# it is already guaranteed to be the correct type
# This saves looping over init when guaranteed to be the correct type.
instance = BrightStarStamps(
(), innerRadius=self._innerRadius, outerRadius=self._outerRadius, metadata=self._metadata
)
Expand All @@ -542,8 +611,8 @@ def _checkRadius(self, innerRadius, outerRadius):
"""
if innerRadius != self._innerRadius or outerRadius != self._outerRadius:
raise AttributeError(
"Trying to mix stamps normalized with annulus radii "
f"{innerRadius, outerRadius} with those of BrightStarStamp instance\n"
f"Trying to mix stamps normalized with annulus radii {innerRadius, outerRadius} with those "
"of BrightStarStamp instance\n"
f"(computed with annular radii {self._innerRadius, self._outerRadius})."
)

Expand All @@ -555,42 +624,40 @@ def _checkNormalization(self, normalize, innerRadius, outerRadius):
nStamps = len(self)
nFluxVals = nStamps - noneFluxCount
if noneFluxCount and noneFluxCount < nStamps:
# at least one stamp contains an annularFlux value (i.e. has been
# normalized), but not all of them do
# At least one stamp contains an annularFlux value (i.e. has been
# normalized), but not all of them do.
raise AttributeError(
f"Only {nFluxVals} stamps contain an annularFlux value.\nAll stamps in a "
"BrightStarStamps instance must either be normalized with the same annulus "
"definition, or none of them can contain an annularFlux value."
f"Only {nFluxVals} stamps contain an annularFlux value.\nAll stamps in a BrightStarStamps "
"instance must either be normalized with the same annulus definition, or none of them can "
"contain an annularFlux value."
)
elif normalize:
# stamps are to be normalized; ensure annular radii are specified
# and they have no annularFlux
# Stamps are to be normalized; ensure annular radii are specified
# and they have no annularFlux.
if innerRadius is None or outerRadius is None:
raise AttributeError(
"For stamps to be normalized (normalize=True), please provide a valid "
"value (in pixels) for both innerRadius and outerRadius."
"For stamps to be normalized (normalize=True), please provide a valid value (in pixels) "
"for both innerRadius and outerRadius."
)
elif noneFluxCount < nStamps:
raise AttributeError(
f"{nFluxVals} stamps already contain an annularFlux value. For stamps to "
"be normalized, all their annularFlux must be None."
f"{nFluxVals} stamps already contain an annularFlux value. For stamps to be normalized, "
"all their annularFlux must be None."
)
elif innerRadius is not None and outerRadius is not None:
# Radii provided, but normalize=False; check that stamps
# already contain annularFluxes
# Radii provided, but normalize=False; check that stamps already
# contain annularFluxes.
if noneFluxCount:
raise AttributeError(
f"{noneFluxCount} stamps contain no annularFlux, but annular radius "
"values were provided and normalize=False.\nTo normalize stamps, set "
"normalize to True."
f"{noneFluxCount} stamps contain no annularFlux, but annular radius values were provided "
"and normalize=False.\nTo normalize stamps, set normalize to True."
)
else:
# At least one radius value is missing; ensure no stamps have
# already been normalized
# already been normalized.
if nFluxVals:
raise AttributeError(
f"{nFluxVals} stamps contain an annularFlux value. If stamps have "
"been normalized, the innerRadius and outerRadius values used must "
"be provided."
f"{nFluxVals} stamps contain an annularFlux value. If stamps have been normalized, the "
"innerRadius and outerRadius values used must be provided."
)
return None