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-10009: Ensure masks are valid from ImageMapReduceTask #63

Merged
merged 1 commit into from Jun 16, 2017

Conversation

djreiss
Copy link
Contributor

@djreiss djreiss commented Jun 5, 2017

Ensure all pixels that are not touched by the reducer are flagged in masks

These pixels are likely to be the result of a bug in the mapper.

Also revent divide by zero when averaging results in weights of 0.

@@ -275,15 +276,36 @@ def run(self, mapperResults, exposure, **kwargs):
wtsView = afwImage.ImageI(weights, item.getBBox())
wtsView.getArray()[isNotNan] += 1

isNan = (np.isnan(newMI.getImage().getArray()) |
np.isnan(newMI.getVariance().getArray()))
Copy link
Contributor

Choose a reason for hiding this comment

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

Since nan * anything is also a nan, you could probably just do np.isnan(newMI.getImage().getArray() * newMI.getVariance().getArra()), and save one loop over all the pixels.

@@ -275,15 +276,36 @@ def run(self, mapperResults, exposure, **kwargs):
wtsView = afwImage.ImageI(weights, item.getBBox())
wtsView.getArray()[isNotNan] += 1

isNan = (np.isnan(newMI.getImage().getArray()) |
np.isnan(newMI.getVariance().getArray()))
if np.sum(isNan) > 0:
Copy link
Contributor

Choose a reason for hiding this comment

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

more over, since you are looping over the whole bool mask array here an again below, with [isNan], by don't you change the above to:
tmp = newMI.getImage().getArray() * newMI.getVariance().getArray()
isNan = np.where(tmp != tmp) # nan never equals anything

That is only 3 loops over the data, and then if if statement can be:
if len(isNan[0]) > 0: # len will fetch the precomputed property,
and you can use isNan[0] and isNan[1] as indicies to mask.getArray(), which will be many fewer indices to loop over as it can do direct array random access.

self.log.info('AVERAGE: Minimum overlap: %f', wts.min())
wtsZero = wts == 0.
self.log.info('AVERAGE: Number of zero pixels: %f (%f%%)', np.sum(wtsZero),
np.mean(wtsZero) * 100.)
Copy link
Contributor

Choose a reason for hiding this comment

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

here you are needlessly looping twice, consider doing:
wtsZeroSum = np.sum(wtsZero)
self.log.info('AVERAGE: Number of zero pixels: %f (%f%%)', wtsZeroSum,
wtsZeroSum/wtsZeroSum.size * 100.)

Copy link
Contributor

Choose a reason for hiding this comment

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

Also just consider using afw statistics directly, I think you can have it build and return multiple statistics all at once.

wtsZero = wts == 0.
newMI.getImage().getArray()[wtsZero] = newMI.getVariance().getArray()[wtsZero] = np.nan
# TBD: set mask to something for pixels where wts == 0. Shouldn't happen. (DM-10009)
if np.sum(wtsZero) > 0:
Copy link
Contributor

Choose a reason for hiding this comment

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

consider moving this if statement higher, so some of the above does not need to get done if it is not true, also reuse the variable I suggested.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

thought about moving it higher, but the logic doesn't work. But your second comment is valid.

# New mask plane - for debugging map-reduced images
mask.addMaskPlane('INVALID_MAPREDUCE')
bad = mask.getPlaneBitMask(['INVALID_MAPREDUCE', 'BAD', 'NO_DATA'])
mask.getArray()[wtsZero] |= bad
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 combine the steps to produce a single mask, you would only need to add a mask plane and populate it once. More over I think this code will fail if np.sum(isNan) >0 and if np.sum(wtsZero) >0, as the mask plane will already exists and an exception will be thrown when trying to add it a second time.

wtsZero = wts == 0.
self.log.info('AVERAGE: Number of zero pixels: %f (%f%%)', np.sum(wtsZero),
np.mean(wtsZero) * 100.)
wts[wtsZero] = 1e-29 # avoid division by zero warnings -- will set to nan and flag them below
Copy link
Contributor

Choose a reason for hiding this comment

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

This line is not needed and loops over pixels needlessly, see below

wtsZero = wts == 0.
self.log.info('AVERAGE: Number of zero pixels: %f (%f%%)', np.sum(wtsZero),
np.mean(wtsZero) * 100.)
wts[wtsZero] = 1e-29 # avoid division by zero warnings -- will set to nan and flag them below
newMI.getImage().getArray()[:, :] /= wts
newMI.getVariance().getArray()[:, :] /= wts
Copy link
Contributor

Choose a reason for hiding this comment

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

numpy ufuncs save us here, you can just do:
np.divide(newMI.getImage().getArray(), wts, out=newMI.getImage().getArray(), where=~wtsZero)
This will do inplace operation (output is the same as input) and only loop over and do a calculation where wtsZero is not true, all in one step.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Only works as:

tmp = newMI.getImage().getArray()
np.divide(tmp, wts, out=tmp, where=~wtsZero)

will do it this way.

newMI.getImage().getArray()[wtsZero] = newMI.getVariance().getArray()[wtsZero] = np.nan
# TBD: set mask to something for pixels where wts == 0. Shouldn't happen. (DM-10009)
if np.sum(wtsZero) > 0:
newMI.getImage().getArray()[wtsZero] = newMI.getVariance().getArray()[wtsZero] = np.nan
Copy link
Contributor

Choose a reason for hiding this comment

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

why are you setting these values to nan here? is that not redundant info with the mask plane bits set? Do we ever want nans in an image? Doesn't that make things like afw statistics need to do more work?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Not sure about your last question. This was for the case where I used wts[wtsZero] = 1e-29 . Since you helped me get rid of that, this line can go now too.

mask = newMI.getMask() # Now check the mask
self.assertGreater(mask.getMaskPlane('INVALID_MAPREDUCE'), 0)
maskBit = mask.getPlaneBitMask('INVALID_MAPREDUCE')
nMasked = np.sum((mask.getArray()[:, :] & maskBit) != 0)
Copy link
Contributor

Choose a reason for hiding this comment

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

Im pretty sure here the operation you want is numpy.bitwise_and, you can avoid the [:,:] in that case.

@@ -261,6 +261,7 @@ def run(self, mapperResults, exposure, **kwargs):
patchMI = item.getMaskedImage()
isNotNan = ~(np.isnan(patchMI.getImage().getArray()) |
np.isnan(patchMI.getVariance().getArray()))
Copy link
Contributor

Choose a reason for hiding this comment

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

See my comment below about multiplying the arrays together to prevent needless loops. If you used np.isfinite instead of np.isnan you can save the inverting loop (if inf is also not a valid value).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

At risk of premature optimization, I will implement everything you've recommended here.

These pixels are likely to be the result of a bug in the mapper.

Prevent divide by zero when averaging results in weights of 0.
Add unit test to check mask
@djreiss djreiss merged commit 9916df4 into master Jun 16, 2017
@ktlim ktlim deleted the tickets/DM-10009 branch August 25, 2018 06:44
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

2 participants