#239 Source detection and SNR computation code #301
Conversation
|
Hi @toshiyukimizuki, |
| from pyklip.kpp.stat.statPerPix_utils import get_image_stat_map_perPixMasking | ||
| from corgidrp.data import Dataset | ||
|
|
||
| def find_source(input_dataset: Dataset, psf=None, fwhm=3.5, nsigma_threshold=5.0): |
There was a problem hiding this comment.
We haven't so far been using type hints in the function calls, but I suppose we could. Should the type instead be corgidrp.data.Dataset?
There was a problem hiding this comment.
Actually I think it should be fine for this to accept an Image as input rather than a dataset since we don't expect the L4_to_TDA.py doesn't have to be scriptable with the walker, so the input type should be corgidrp.data.Image (and the variable name should reflect that too).
| Args: | ||
| input_dataset (Dataset): Input dataset containing images. | ||
| psf (np.ndarray, optional): 2D PSF array. If None, a Gaussian PSF is generated. | ||
| fwhm (float, optional): Full-width at half-maximum of the PSF. |
There was a problem hiding this comment.
specify the units here, I think it's "pixels", right?
| @@ -0,0 +1,89 @@ | |||
| import numpy as np | |||
There was a problem hiding this comment.
I think we should create a file called l4_to_TDA.py and put this there since it's an L4 to TDA processing step.
| sn_source, xy_source = [], [] | ||
|
|
||
| # Iteratively find sources above the SNR threshold | ||
| while np.nanmax(image_snmap) >= nsigma_threshold: |
There was a problem hiding this comment.
Could we modularize all these steps a bit more (e.g. the SNR detection, then then PSF scaling/subtraction)? I could imagine a case where we might want to use these steps with a bit more manual control. For example, the PSF's morphology will change close to the inner and outer edges of the field of view, so being able to adjust the input PSF could be helpful. I still think that we want it packaged nicely into a function like this, but having things be more modular could also allows us to carry out some of these steps manually.
|
|
||
| # Store detected sources in FITS header | ||
| for i in range(len(sn_source)): | ||
| input_dataset[0].pri_hdr[f'snyx{i:03d}'] = f'{sn_source[i]:5.1f},{xy_source[i][0]:4d},{xy_source[i][1]:4d}' |
There was a problem hiding this comment.
Let's add it to the extension header ext_hdr instead. The primary header will mostly be reserved for things set by the telescope.
| @@ -0,0 +1,173 @@ | |||
| import os, glob, copy | |||
There was a problem hiding this comment.
Please take a look at the structure of other test functions - we're using a pytest infrastructure which requires that we have functions start with test_* and make specific assertion tests.
| nondetection = np.vstack((nondetection, np.delete(np.column_stack((sn_rand, y_rand, x_rand)), idx1, axis=0))) | ||
| misdetection = np.vstack((misdetection, np.delete(snyx, idx2, axis=0))) | ||
|
|
||
| # Print summary of detection results |
There was a problem hiding this comment.
At the end here rather than printing things or plotting things we should make assertion statements to show that we've recovered the input x,y positions and the expected SNR
| nsigma_threshold = 5.0 | ||
|
|
||
| # Load the dataset from GPI data file | ||
| file_gpi = mockfilepath+'NoCompNoDisk-GPIcube.fits' |
There was a problem hiding this comment.
I actually quite like what you've done here, but in order to keep the size of the repository small we're trying to minimize the number of data files that we add. Rather than adding more data to the repository here we could just simulate a noise field. For this level of test I think its ok if the noise is just Gaussian, or if you wanted a radially changing Gaussian noise field that could be fine too. For the end to end tests we might have higher fidelity speckle residual fields.
| # Categorize results into detections, non-detections, and misdetections | ||
| detection = np.vstack((detection, np.hstack((np.column_stack((sn_rand, y_rand, x_rand))[idx1], snyx[idx2])))) | ||
| nondetection = np.vstack((nondetection, np.delete(np.column_stack((sn_rand, y_rand, x_rand)), idx1, axis=0))) | ||
| misdetection = np.vstack((misdetection, np.delete(snyx, idx2, axis=0))) |
There was a problem hiding this comment.
I like that you're injecting sources below the SNR threshold to make sure we don't detect them - that's a good test of the threshold.
| inputflux_rand*= 1.5 # somehow needed ? | ||
| ##### ##### ##### | ||
| x_rand = radius_rand * np.cos(np.radians(pa_rand)) + dataset_center[1] | ||
| y_rand = radius_rand * np.sin(np.radians(pa_rand)) + dataset_center[0] |
There was a problem hiding this comment.
I don't think we want to randomly inject point sources with random SNRs, we need to define them specifically so that we can show that we can recover them.
Describe your changes
Submission of code to create snmap and do point source detection, as well as code to test it
corgidrp/find_source.py
tests/test_find_source.py
tests/test_data/NoCompNoDisk-GPIcube.fits
Reference any relevant issues (don't forget the #)
#239
Checklist before requesting a review