kammerje / spaceKLIP

Pipeline for reducing JWST high-contrast imaging data. Published in Kammerer et al. 2022 and Carter et al. 2022.
https://ui.adsabs.harvard.edu/abs/2022SPIE12180E..3NK/abstract
MIT License
16 stars 9 forks source link

Add sqrt_stretch and custom mask to align steps #183

Closed wbalmer closed 4 weeks ago

wbalmer commented 1 month ago

See discussion in #114 . For data with more complex speckle patterns (e.g. NIRCam bar NARROW coronagraphy), it can be necessary to emphasize the speckle pattern itself rather than the central leakage term. I've found a combination of two strategies reproduces the commanded shifts from the SGD: masking the central leakage using a custom mask during the align_frames step, and taking the square root of the images before cross correlation. This PR introduces non-intrusive toggles to do both of these things.

wbalmer commented 1 month ago

Bug with take_sqrt option in find_nircam_centers (now fixed)

AarynnCarter commented 1 month ago

@wbalmer, looks great! I did have one suggestion that could be a nice adjustment. Instead of a binary do/do not do a square root. Could we instead provide an exponent to scale the data by? The square root worked in your case, but others might find in the future that something like a cube root could be more useful, or who knows, even an exponent >1. With this, we can set the default exponent==1, and can cut down on some of the added if statements.

Let me know what you think!

wbalmer commented 1 month ago

Hi Aarynn. I see your point about generalizing/simplifying the statement. I suppose the problem is that np.power, the function I'd use to generalize this, returns nans when negative values (which there will be in median subtracted data) are raised to non integer exponents (like 1/2). np.sqrt and np.cbrt provides an interface for returning the principle root. So to implement something like you describe, I feel like it would end up being about as many if/else statements as there are currently (a case for exp=1/2 or 1/3 and then a generalized case for exp>=1). I can do this, if you'd still like to see the generalization

AarynnCarter commented 1 month ago

@wbalmer Okay gotcha, in that case maybe we can strike a middle ground. For the nominal case (exp=1), we could use your current else statement

else:
    img1 = datasub* masksub
    img2 = model_psf* masksub

But for anything else, we can adjust to the following:

if take_sqrt:
    img1 = np.sqrt(np.abs(datasub))* masksub
    img2 = np.sqrt(np.abs(model_psf))* masksub

to

if scale_exponent != 1:  #Or some other sensible name
    img1 = np.power(np.abs(datasub), scale_exponent)* masksub
    img2 = np.power(np.abs(model_psf), scale_exponent)* masksub

I feel like this keeps us to the current number of if statements, but more general, and avoids using the absolute function for the nominal case. For exponents >1, it won't matter that we've taken the absolute value as the result will be the same. Do you think that would make sense?

wbalmer commented 1 month ago

Sounds good! I'll test this now.

wbalmer commented 1 month ago

Okay, tested in both the shft_exp = 1 and shft_exp = 0.5 cases and expected behavior was reproduced. Should be good to go

wbalmer commented 4 weeks ago

Should be fixed!