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.

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-22221: Add subtask to find and mask satellite trails #370

Merged
merged 1 commit into from Aug 28, 2020

Conversation

cmsaunders
Copy link
Contributor

This commit adds a new subtask in findSatellites.py and a call to
that task in CompareWarpAssembleCoaddTask, where it is used to
detect satellite trails or other linear features in difference
images. This will only run if the config option
doFilterMorphological is set to True.

Copy link
Contributor

@yalsayyad yalsayyad left a comment

Choose a reason for hiding this comment

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

Just the unit test


def setUp(self):
self.fst = FindSatellitesTask()
self.fst.config.dChi2Tolerance = 1e-10
Copy link
Contributor

Choose a reason for hiding this comment

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

We don't change the state of a config after the Task has been instantiated. In fact, when a cmdlineTask is launched it calls Config.freeze() to make config to make it immutable: https://github.com/lsst/pipe_base/blob/master/python/lsst/pipe/base/argumentParser.py#L725 This isn't enforced when using as a regular task, which is why this works. But the standard way of doing what you want to do is to edit the config and then instantiate the Task

self.config = FindSatellitesTask.ConfigClass()
self.config.dChi2Tolerance = 1e-10
self.fst = FindSatellitesTask(config=self.config)

*In this example I'm saving the config as an instance variable because you'll reuse it later in line 43

self.testy = 600
self.exposure = lsst.afw.image.ExposureF(self.testy, self.testx)
self.exposure.maskedImage.image.array = np.random.randn(self.testx, self.testy).astype(np.float32)
self.exposure.maskedImage.variance.array = np.ones((self.testx, self.testy)).astype(np.float32)
Copy link
Contributor

Choose a reason for hiding this comment

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

self.exposure.variance.set(1)

self.testx = 500
self.testy = 600
self.exposure = lsst.afw.image.ExposureF(self.testy, self.testx)
self.exposure.maskedImage.image.array = np.random.randn(self.testx, self.testy).astype(np.float32)
Copy link
Contributor

Choose a reason for hiding this comment

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

A stack-y way of doing this is:

rand = lsst.afw.math.Random()
lsst.afw.math.randomGaussianImage(exposure.image, rand)

which also has the benefit of making your test deterministic, because there is a default seed in lsst.afw.math.random.Random(algorithm: str, seed: int=1)


reshapeBinning = self.fst.processMaskedImage(self.exposure)
with self.assertWarns(Warning):
scipyBinning = self.fst.processMaskedImage(self.exposure, checkMethod=True)
Copy link
Contributor

Choose a reason for hiding this comment

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

By calling this parameter checkMethod you are communicating that this option only exist to support these unit tests? Is that true? By reading your inline comments in the task, I'd think that it's a slower, more flexible option.

@@ -0,0 +1,94 @@
import unittest
Copy link
Contributor

Choose a reason for hiding this comment

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

Need the preamble. See https://developer.lsst.io/stack/license-and-copyright.html?#python-preamble on how to get one from @sqrbot-jr

scipyBinning = self.fst.processMaskedImage(self.exposure, checkMethod=True)
self.assertAlmostEqual(reshapeBinning.tolist(), scipyBinning.tolist())

self.fst.config.imageBinning = None
Copy link
Contributor

Choose a reason for hiding this comment

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

People look to unit tests as examples, since they're guaranteed to work. Whereas this works, it sets a bad example. Please instantiate a new task rather than modifying a config.

self.config.imageBinning = None
task = FindSatelliteTask(config=self.config)
nobinImage = task.processMaskedImage(self.exposure)


result = self.fst.run(testExposure, binaryImage=binaryImage)
self.assertEqual(len(result.lines), 1)
self.assertLess(abs(input['rho'] - result.lines[0]['rho']), 0.01)
Copy link
Contributor

Choose a reason for hiding this comment

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

self.assertAlmostEqual(input['rho'], result.lines[0]['rho'], delta=0.01) would prob be more readable,


# Make image with line and check that input line is recovered:
testExposure = self.exposure.clone()
input = np.array((150, 45, 3), dtype=[('rho', float), ('theta', float), ('sigma', float)])
Copy link
Contributor

Choose a reason for hiding this comment

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

tests/test_findSatellites.py Outdated Show resolved Hide resolved
Copy link
Contributor

@natelust natelust left a comment

Choose a reason for hiding this comment

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

I'm going to leave these comments now, so everyone can look over them. I may read over this again soon to see if there is anything else I see.

#
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
import lsst.newhoughtransform
Copy link
Contributor

Choose a reason for hiding this comment

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

In general a space between lsst imports and standard imports is nice

Copy link
Member

Choose a reason for hiding this comment

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

Is lsst.newhoughtransform a new LSST supported package that wraps the new third party hough transform package?

Copy link
Contributor

Choose a reason for hiding this comment

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

RFC-680. This code is just undergoing first round review while that RFC is implemented.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

lsst.newhoughtransform is a placeholder until we get the third-party hough transform package implemented, but the functionality should be the same.

from sklearn.cluster import KMeans
import warnings

__all__ = ["FindSatellitesConfig", "FindSatellitesTask"]
Copy link
Contributor

Choose a reason for hiding this comment

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

I know we are not great at following this convention, but I think all should be above the imports

Copy link
Contributor

Choose a reason for hiding this comment

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

a line following the given coordinates.
"""
def __init__(self, data, weights, mask, invSigmaInit, line=None):
"""
Copy link
Contributor

Choose a reason for hiding this comment

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

its weird, but the parameters documentation of init should be with the class level documentation not in init

2d array with mask
invSigmaInit : float
Initial guess for inverse of Moffat sigma parameter
line : np.ndarray (optional)
Copy link
Contributor

Choose a reason for hiding this comment

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

indicate somewhere that it will default to None

Copy link
Contributor

Choose a reason for hiding this comment

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

Change the data type to reflect what it really is, a numpy.recarray.

This part of the comment is a little program design heavy, and you might know it already, but if not I hope it helps. I would strongly suggest you consider using a different data structure that you define. Record arrays can be useful, but often they are an anti pattern. If you created a simple data structure using namedtuple (we talked about this in bookclub) or a dataclass (we will talk about this soon in book club) it would make it easier for other people (or future you) to work with this code.

Data structures you create will have their own types which would make function signatures easier to understand, i.e. know exactly what needs to be passed. These data structures have their own schema that people will be able to see, and will be required to populate when the object is created.

For instance if you had a datastructure like

# There is a version of this using typing.NamedTuple, but that involves type annotations
# that we have not talked about yet, so using what book club covered.
LineOutline = namedtuple("LineOutline", ("rho", "theta", "invSigma"))

you could be sure that the line argument passed to this function would always has theta defined. This is also useful for the user, because they know what "schema" the object might have. Because recarrays are basically arbitrary containers, a user would have no idea without reading your code either here or where one is produced what this object should have in it.

Defined datastructures also have the advantage that you can check their type if you really wanted to validate that something is exactly the type you want. They are also more reusable in other locations, and provide a single place where you can change a definition rather than relying on identical implementations of recarray creation at multiple places in a code base. This will make refactoring code later much easier.

Ok long comment over, this is an optional suggestion, but I hope you find it useful to think about.

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 have implemented a Line dataclass and a LineSet class as discussed. These should remove all of the cases where there are recarrays or separate variables for rho, theta, and sigma, so this also addresses a lot of other comments on findSatellites.py.

yrange = np.arange(ymax) - ymax / 2.
self.rhoMax = ((0.5 * ymax)**2 + (0.5 * xmax)**2)**0.5
self.xmesh, self.ymesh = np.meshgrid(xrange, yrange)
self.mask = (weights != 0)
Copy link
Contributor

Choose a reason for hiding this comment

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

The () are not needed here

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 just do this because I think it is more readable. I will remove the parentheses if it is against the convention.

# change more than the allowed bin in rho or theta:
if ((abs(fit.rho - line.rho) > 2 * self.config.rhoBin) or
(abs(fit.theta - line.theta) > 2 * self.config.thetaBin)):
fitSuccess = False
Copy link
Contributor

Choose a reason for hiding this comment

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

I would probably restructure this if statement as

if not fitSuccess or
    abs(fit.rho - line.rho) > 2 * self.config.rhoBin or
    abs(fit.theta - line.theta) > 2 * self.config.thetaBin:
       continue

and drop the second conditional below

Copy link
Contributor

Choose a reason for hiding this comment

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

(ignore the weird formatting in the above, use standard pep8 spacing)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

fitSuccess originally comes from a few lines further up. I wanted to only have on continue statement for all fit failures, so that is why I have formatted it this way, if that makes sense.

lineModel.setLineMask(fit)
finalModel = lineModel.makeProfile(fit.rho, fit.theta, fit.invSigma)
# Take absolute value, as trails are allowed to be negative
finalModelMax = (abs(finalModel)).max()
Copy link
Contributor

Choose a reason for hiding this comment

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

you should not need the extra ()

lineFits.append([fit.rho, fit.theta, sigma, chi2, finalModelMax])
finalLineMasks.append(finalLineMask)

if len(lineFits) == 0:
Copy link
Contributor

Choose a reason for hiding this comment

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

See comment about about a container class for all of this. The constructor of that type can be setup such that you would not need the if else in this case.

line : np.ndarray (optional)
Guess for position of line, with datatype names `rho`, `theta`, and
`sigma`
"""
Copy link
Contributor

Choose a reason for hiding this comment

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

When you build a class it is important to think about what story, its interface, that you want told to the rest of the world. If a user gets passed an instance of your class what are they supposed to know about it, and how should they use it. This is part of documentation yes, but also in what methods/members your class presents to them. Should a user look at an instance of your class and do instance.sigma to get sigma info? Is instance.mask how they should be getting the mask, or should it be the result of a method call? Anything that would be an implementation detail should have its name start with an underscore.

For instance would a user ever call setLineMask on some instance they were passed, or is this a method you used in your class while doing some calculation the user does care about. If the later, the name should be _setLineMask. Same is true for the attributes. I am pretty sure users should not be looking at initLine so it should be _initLine.

Setting up your names like this does a few things.

  • it provides a smaller api that you are reasonably sure will not change, so you only have to keep a guarantee on those interfaces, and you are free to refactor the rest without fearing you will break something someone was using.

  • It is self documenting, as in when a user does dir they are presented with only the interfaces that are important for them, helping them in discoverability

  • It helps other developers know what you intended when you wrote the code. If someone comes along and sees 2 public methods that can be sure those are the meat and potatoes and to focus on them when developing against your library, or writing a new compatible one.

Im sure if you talked to others there would be more things on this list, but these are the ones that sprang to mind now. In LSST we have not been particularly good about this in the past, but it would be really great if people coming up to speed started doing this more, so that our newest code, that others are likely to see and use will already be cleaner and more maintainable.

fitResult = self.findSatellites.run(warpDiffExp, slateIm)
satMask = afwImage.Mask(coaddBBox)
satMask.array[fitResult.mask] = 1
spanSetSatellite = afwGeom.SpanSet.fromMask(satMask)
Copy link
Contributor

Choose a reason for hiding this comment

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

I am pretty sure what you want to do here is to call the split method on the spanset that is returned from fromMask. fromMask returns one giant SpanSet that basically just maps exactly onto the mask. But I think what you want here is a bunch of individual SpanSets with all the rest dropped. This is basically a bunch of footprints of all the satellites that are to be ignored. @yalsayyad may know better than me about what this algorithm expects though.

Also, if possible, could your findSatellites.run just use afwMasks? Or return one maybe?

Copy link
Contributor

Choose a reason for hiding this comment

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

It probably doesn't matter if we treat one satellite at a time or all-the-satellites as one.... This does explain why list.append(list) below worked...Because it wasn't a list!

Unrelated, can you move this doFilterMorphological block above the doPrefilterArtifacts block. I think I want to prefilter artifacts contained by your satellites too eventually.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

For moving the doFilterMorphological block above the doPrefilterArtifacts block, there are a couple issues. First, the satellite code uses slateIm, which is created from the prefiltered artifacts. I think interpolated regions are prefiltered here (as SUSPECT, but maybe I am wrong), and it is important that these are removed before running findSatellites.
Second, I don't think prefilterArtifacts works well for the satellite trails. I tried what I think you are suggesting when I was figuring out the implementation, but ran into problems. Often the satellite trail intersects with another large footprint, like a bright star. Then, the "bad" fraction of the total footprint ends up being less than self.prefilterArtifactsRatio, so the whole footprint is kept, not filtered out and the satellite trail still goes into the final image.

# Initial estimate should be quite close: fit is deemed unsuccessful if rho or theta
# change more than the allowed bin in rho or theta:
if ((abs(fit.rho - line.rho) > 2 * self.config.rhoBin) or
(abs(fit.theta - line.theta) > 2 * self.config.thetaBin)):
Copy link
Contributor

Choose a reason for hiding this comment

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

E128: indent one more

@@ -2232,6 +2256,8 @@ def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList):
if spanSetList:
filteredSpanSetList = self.filterArtifacts(spanSetList, epochCountImage, nImage,
templateFootprints)
if self.config.doFilterMorphological and (spanSetBadMorphoList[i].getArea() != 0):
filteredSpanSetList.append(spanSetBadMorphoList[i])
Copy link
Contributor

Choose a reason for hiding this comment

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

You mean += instead of append?

>>> a = [1,2,3]
>>> b = [4,5,6]
>>> a.append(b)
>>> a
[1, 2, 3, [4, 5, 6]]
>>> a = [1,2,3]
>>> a += b
>>> a
[1, 2, 3, 4, 5, 6]

I think this whole block should get its own loop, in the scientifically impossible but programmatically possible chance that there are no artifacts but there is a satellite trail.

        if self.config.doFilterMorphological:
            for i, (satellites, artifacts) in enumerate(zip(spanSetBadMorphoList, spanSetArtifactList)):
                artifacts += satellites

fitResult = self.findSatellites.run(warpDiffExp, slateIm)
satMask = afwImage.Mask(coaddBBox)
satMask.array[fitResult.mask] = 1
spanSetSatellite = afwGeom.SpanSet.fromMask(satMask)
Copy link
Contributor

Choose a reason for hiding this comment

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

It probably doesn't matter if we treat one satellite at a time or all-the-satellites as one.... This does explain why list.append(list) below worked...Because it wasn't a list!

Unrelated, can you move this doFilterMorphological block above the doPrefilterArtifacts block. I think I want to prefilter artifacts contained by your satellites too eventually.

"""

ConfigClass = FindSatellitesConfig
_DefaultName = "findSatellitesConfig"
Copy link
Contributor

Choose a reason for hiding this comment

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

_DefaultName = "findSatellites"

dtype=float,
default=10.**-1,
)
imageBinning = pexConfig.Field(
Copy link
Contributor

Choose a reason for hiding this comment

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

We call this binSize in SubtractBackgroundTask and other parts of the stack. imageBin would keep it consistent with your thetaBin etc.. I'd support rhoBinSize, imageBinSize

dtype=float,
default=2,
)
minImageSignaltoNoise = pexConfig.Field(
Copy link
Contributor

Choose a reason for hiding this comment

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

essentially the same thing as thresholdValue in SourceDetectionTask.
https://github.com/lsst/meas_algorithms/blob/master/python/lsst/meas/algorithms/detection.py#L66 Your processMaskedImage = SourceDetectionTask (with thresholdType='pixel_stdev') + binning.

Copy link
Contributor

Choose a reason for hiding this comment

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

I still think minImageSignaltoNoise is not going to register as the detection threshold for stack users. And since I just called it the detection threshold without thinking about it, I'd recommend detectionThreshold.

imageBinning = pexConfig.Field(
doc="Number of pixels by which to bin image",
dtype=int,
default=2,
Copy link
Contributor

Choose a reason for hiding this comment

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

Out of curiosity, why bin the image? Does it make the kht faster?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It does make the kht faster, but it also helps make satellite trails sharper.

"""Configuration parameters for `FindSatellitesTask`
"""
minimumKernelHeight = pexConfig.Field(
doc="minimum height of the satellite finding kernel relative to the tallest kernel",
Copy link
Contributor

Choose a reason for hiding this comment

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

"satellite-finding". Start with a Capital M. Units? When I hear "relative," I think fractional, rather than additive offset.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It is fractional.

mask = maskedImage.getMask()

detectionMask = ((mask.array & mask.getPlaneBitMask("DETECTED")))
badMask = ((mask.array & mask.getPlaneBitMask("NO_DATA")) |
Copy link
Contributor

Choose a reason for hiding this comment

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

For example, that would look like:

badPixelMask = mask.getPlaneBitMask(self.config.badMaskPlanes)
badMask = (mask.array & badPixelMask) >  0

badMask is already a bool.

fitWeights = np.copy(weights)
fitWeights[~fitMask] = 0

if binning is not None:
Copy link
Contributor

Choose a reason for hiding this comment

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

In the long run, this isn't the right place for a binning algorithm. Consider opening a ticket to move this lower into the stack. In the stack we have afwMath.binImage() which can only handle means. Most the image binning used by the stack is actually inafwMath.makeBackground which can do clipped means and medians. https://github.com/lsst/afw/blob/master/doc/lsst.afw.math/Background-example.rst

Also, why do you need inverse-variance weighted means when binning this image?

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 looked at using afwMath.binImage() here, but I think it would be more convoluted. Even without doing a weighted mean, you need the binned weights in order to get the S/N at the end. I think this means that, to use afwMath.binImage(), you would have make an afwImage ImageF or similar, put the weights in that, and run afwMath.binImage() on that.

Separately, doing a weighted mean takes care of removing bad pixels, though it is probably not crucial otherwise.

Copy link
Contributor

@yalsayyad yalsayyad left a comment

Choose a reason for hiding this comment

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

Adding comments I started a while back. 4 comments are from before the module got renamed to maskStreaks. I took a screenshot if they don't show up. We can talk about higher level stuff today.

2d array of data
weights : np.ndarray
2d array of weights
mask : np.ndarray
Copy link
Contributor

Choose a reason for hiding this comment

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

There's no mask argument in the constructor

dtype=float,
default=0.2,
)
nSigmas = pexConfig.Field(
Copy link
Contributor

Choose a reason for hiding this comment

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

nSigma would be easier to remember IMO.

dtype=float,
default=2,
)
minImageSignaltoNoise = pexConfig.Field(
Copy link
Contributor

Choose a reason for hiding this comment

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

I still think minImageSignaltoNoise is not going to register as the detection threshold for stack users. And since I just called it the detection threshold without thinking about it, I'd recommend detectionThreshold.

default=("NO_DATA", "INTRP", "BAD", "SAT", "EDGE")
)
detectedPlanes = pexConfig.ListField(
doc="Pixels that were detected above threshold in image",
Copy link
Contributor

Choose a reason for hiding this comment

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

Would you add what this is used for in the docstring?

)

def processMaskedImage(self, maskedImage, forceSlowBin=False):
"""Make binary image array from maskedImage object
Copy link
Contributor

Choose a reason for hiding this comment

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

This summary string describes the input and output but doesn't really describe what it does which is (I think) "Make a detection map" or "Bin and detect" The stack has various detection tasks and binning routines, and we need to make clear why we need another one and eventually move it to where the others live.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This function has been renamed "setDetectionMask" and downgraded to a helper function outside the class.


Parameters
----------
maskedImage : `lsst.afw.image.Exposure`
Copy link
Contributor

Choose a reason for hiding this comment

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

Calling a lsst.afw.image.Exposure maskedImage is confusing.

Returns
-------
out_data : `np.ndarray`
2-d binary image of pixels above the signal-to-noise threshold.
Copy link
Contributor

Choose a reason for hiding this comment

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

np.zeros makes an array of float64s which is what it returns now, right? Not a great return type for a boolean map.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This has been removed because the function no longer returns anything, it just sets a mask plane in the input image.

ind = binnedDenominator != 0
np.divide(binnedNumerator, binnedDenominator, out=binnedData, where=ind)
binnedWeight = binnedDenominator
binMask = (binnedData * binnedWeight**0.5) > self.config.minImageSignaltoNoise
Copy link
Contributor

Choose a reason for hiding this comment

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

binMask is a boolean array which you could return as-is


Returns
-------
#lines : `np.ndarray'
Copy link
Contributor

Choose a reason for hiding this comment

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

Forgot to delete these lines?

- ``originalLines``: lines identified by kernel hough transform
- ``lineClusters``: lines grouped into clusters in rho-theta space
- ``lines``: final result for lines after line-profile fit
- ``mask``: 2-d boolean mask where detected lines=0
Copy link
Contributor

Choose a reason for hiding this comment

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

I'd expect lines to be True/1 and empty space to be False/0

sigma: float = 0


class LineSet:
Copy link
Contributor

Choose a reason for hiding this comment

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

This class does not implement the set protocol, so it probably should not be named set. Maybe LineCollection?



class LineSet:
"""Set of `Line` objects.
Copy link
Contributor

Choose a reason for hiding this comment

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

If you change the name change this

else:
self.sigmas = np.zeros(len(self.rhos))

self.lines = [Line(rho, theta, sigma) for (rho, theta, sigma) in
Copy link
Contributor

Choose a reason for hiding this comment

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

Do you want to allow there to be duplicate identical lines in this object?

Copy link
Contributor

Choose a reason for hiding this comment

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

I dont think you are exposing the fact that there is a lines attribute outside of the class (since you have a __len__ and __getitem__) so you should make this self._lines


def __getitem__(self, index):
return self.lines[index]

Copy link
Contributor

Choose a reason for hiding this comment

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

You are implementing __len__ and __getitem__ which is enough for iteration, but it would be more efficient if you were to implement __iter__ as well:

def __iter__(self):
    return iter(self._lines)

return self.lines[index]

def __repr__(self):
return ", ".join(str(line) for line in self.lines)
Copy link
Contributor

Choose a reason for hiding this comment

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

This could get REALLY long if you print this out, consider wrapping this in a textwrap.shorten call with some width like:

return textwrap.shorten(<joined string>, width=80, placeholder="...")

This will be more friendly to people that print it out

"""
self.rhos *= scalingFactor
self.lines = [Line(rho, theta, sigma) for (rho, theta, sigma) in
zip(self.rhos, self.thetas, self.sigmas)]
Copy link
Contributor

Choose a reason for hiding this comment

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

Consider adding a rescaleRho to the Line class so you can just say

for line in self._lines:
    line.rescale(scalingFactor)

This will save you from creating the objects all over again

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 removed the need for rescaling here.

newLine : `Line`
`Line` to add to current set of lines
"""
self.lines.append(newLine)
Copy link
Contributor

Choose a reason for hiding this comment

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

Do you want to take ownership of the line you are adding, so if it changes on the outside it does not change inside your class? Or for performance do you want to want to rely on people to be good?

self.lines.append(copy.copy(newLine))

"""

# Renormalize variables so that expected standard deviation in a
# cluster is 1.
Copy link
Contributor

Choose a reason for hiding this comment

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

I'm not sure I get how what you are doing below follows from the comment you have here, but I may just be missing something.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This has been rewritten.

# change more than the allowed bin in rho or theta:
if ((abs(fit.rho - line.rho) > 2 * self.config.rhoBinSize)
or (abs(fit.theta - line.theta) > 2 * self.config.thetaBinSize)):
fitSuccess = False
Copy link
Contributor

Choose a reason for hiding this comment

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

If you changed this var name to fitFailed it would be more natural to say

if fitFailed:
    continue


Returns
-------
#lines : `np.ndarray'
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't think there should be # here

Copy link
Contributor

@natelust natelust left a comment

Choose a reason for hiding this comment

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

A few more minor things, but your almost there!


@property
def rhos(self):
return np.array([line.rho for line in self._lines])
Copy link
Contributor

Choose a reason for hiding this comment

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

change this to np.array(line.rho for line in self._lines). This will use a generator expression to build the array and saves you from creating an intermediate list that just gets thrown away.


@property
def thetas(self):
return np.array([line.theta for line in self._lines])
Copy link
Contributor

Choose a reason for hiding this comment

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

same as above

finalModel : np.ndarray
Model for line profile
"""
model, dmodel = self._makeMaskedProfile(line, fitFlux=fitFlux)
Copy link
Contributor

Choose a reason for hiding this comment

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

Im surprised flake8 did not warn you about this. if you are not using dmodel, you should assign it to _ instead:
model, _ = self._makeMaskedProfile(line, fitFlux=fitFlux)

dChi2 = oldChi2 - chi2
cholesky = scipy.linalg.cho_factor(A)
dx = scipy.linalg.cho_solve(cholesky, b)

Copy link
Contributor

Choose a reason for hiding this comment

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

You could save from defining a new line_search function on each interation of the while if you factored it outside the loop and had it take a second parameter dx and then pass dx in the brent call:

def line_search(c, dx):
    testx = x - c * dx
    testLine = Line(testx[0], testx[1], testx[2]**-1)
   return self._lineChi2(testLine, grad=False)

while abs(dChi2) > dChi2Tol:
    ....
    factor, fmin, _, _ = scipy.optimize.brent(line_search, args=(dx,), full_output=True, tol=0.05)

- ``mask``: 2-d boolean mask where detected lines are True
"""
mask = maskedImage.getMask()
detectionMask = ((mask.array & mask.getPlaneBitMask(self.config.detectedMaskPlane)))
Copy link
Contributor

Choose a reason for hiding this comment

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

two too many ()

Parameters
----------
maskedImage : `lsst.afw.image.maskedImage`
The image in which to search for streaks.
Copy link
Contributor

Choose a reason for hiding this comment

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

maybe make a note that the maskedImage mask detection plane is expected to be populated in a way required by your task

"""
filterData = image.astype(int)
cannyData = canny(filterData, low_threshold=0, high_threshold=1, sigma=0.1)
return cannyData
Copy link
Contributor

Choose a reason for hiding this comment

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

you can totally leave this as is, I just want to highlight you can save lines, typing, and readying by putting

return canny(filterData, low_threshold=0, high_threshold=1, sigma=0.1)

default=2,
)
invSigma = pexConfig.Field(
doc="Moffat sigma parameter (in units of pixels) describing the "
Copy link
Contributor

Choose a reason for hiding this comment

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

inverse of the sigma parameter?

# see <http://www.lsstcorp.org/LegalNotices/>.
#

__all__ = ["MaskStreaksConfig", "MaskStreaksTask"]
Copy link
Contributor

Choose a reason for hiding this comment

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

if you want setDetectionMask to be available for other code (like tests) it must be added to __all__, I dont know how the tests were working w/o this.

self.maskName = "STREAK"
self.detectedPlane = "DETECTED"

"""def test_input(self):
Copy link
Contributor

Choose a reason for hiding this comment

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

Something is up with this test, it is commented out? Should you put it back in? Should it be taken out? As an fyi you can use unittest.skip decorator to skip a test while you are testing and debugging. The final test output will indicate that one test was skipped so it is not forgotten like a comment.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

No, this test should have been removed. It is no longer relevant, since we don't allow a separate "binary image" anymore.

Defaults to None, in which case only data with `weights` = 0 is masked
out.
"""
def __init__(self, data, weights, line=None):
Copy link
Contributor

Choose a reason for hiding this comment

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

https://developer.lsst.io/python/numpydoc.html#docstrings-of-classes-should-be-followed-but-not-preceded-by-a-blank-line

Class docs are followed by blank line. Method docs are not. My linter picks this up as E301.

sigmas : np.ndarray (optional)
Array of `Line` sigma parameters
"""
def __init__(self, rhos, thetas, sigmas=None):
Copy link
Contributor

Choose a reason for hiding this comment

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

https://developer.lsst.io/python/numpydoc.html#docstrings-of-classes-should-be-followed-but-not-preceded-by-a-blank-line

Class docs are followed by blank line. Method docs are not. My linter picks this up as E301.

give the same result.
binning : int (optional)
Number of pixels by which to bin image
detectedPlanes : str (optional)
Copy link
Contributor

Choose a reason for hiding this comment

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

detectedPlane vs detectedPlanes

This commit adds a new subtask in findSatellites.py and a call to
that task in CompareWarpAssembleCoaddTask, where it is used to
detect satellite trails or other linear features in difference
images. This will only run if the config option
doFilterMorphological is set to True.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants