Skip to content

Commit

Permalink
Merge pull request #575 from VChristiaens/master
Browse files Browse the repository at this point in the history 
Minor bug fix for PACO + added options for correct_ann_outliers function
  • Loading branch information
@VChristiaens
VChristiaens authored Feb 10, 2023
2 parents c7ed181 + c0ae8b2 commit 0ac5797
Show file tree
Hide file tree
Showing 5 changed files with 119 additions and 93 deletions.
2 changes: 2 additions & 0 deletions README.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@ VIP - Vortex Image Processing package
.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.7499314.svg
:target: https://doi.org/10.5281/zenodo.7499314

.. image:: https://img.shields.io/badge/EMAC-2207--116-blue
:target: https://emac.gsfc.nasa.gov/?cid=2207-116

::

Expand Down
2 changes: 1 addition & 1 deletion vip_hci/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
__version__ = "1.3.6"
__version__ = "1.4.0"


from . import preproc
Expand Down
160 changes: 89 additions & 71 deletions vip_hci/invprob/paco.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,9 @@ def __init__(self,

@abstractmethod
def PACOCalc(self,
phi0s : np.ndarray,
use_subpixel_psf_astrometry : Optional[bool] = True,
cpu : Optional[int] = 1) -> None:
phi0s: np.ndarray,
use_subpixel_psf_astrometry: Optional[bool] = True,
cpu: Optional[int] = 1) -> None:
"""
This function is algorithm dependant, and sets up the actual
calculation process.
Expand Down Expand Up @@ -194,11 +194,11 @@ def PACOCalc(self,
"""

def run(self,
cpu : Optional[int] = 1,
imlib : Optional[str] = 'vip-fft',
interpolation : Optional[str]= 'lanczos4',
keep_center : Optional[bool] = True,
use_subpixel_psf_astrometry : Optional[bool] = True) -> Tuple[np.ndarray, np.ndarray]:
cpu: Optional[int] = 1,
imlib: Optional[str] = 'vip-fft',
interpolation: Optional[str] = 'lanczos4',
keep_center: Optional[bool] = True,
use_subpixel_psf_astrometry: Optional[bool] = True) -> Tuple[np.ndarray, np.ndarray]:
"""
Run method of the PACO class. This function wraps up the PACO
calculation steps, and returns the snr and flux estimates as
Expand Down Expand Up @@ -587,21 +587,24 @@ def flux_estimate(self,
# Create arrays needed for storage
# Store for each image pixel, for each temporal frame an image
# for patches: for each time, we need to store a column of patches
normalised_psf,norm,fwhm = normalize_psf(self.psf,
fwhm='fit',
size=None,
threshold=None,
mask_core=None,
model='airy',
imlib='vip-fft',
interpolation='lanczos4',
force_odd=False,
full_output=True,
verbose=self.verbose,
debug=False)

psf_mask = create_boolean_circular_mask(normalised_psf.shape, radius=self.fwhm)
hoff = np.zeros((self.num_frames,self.num_frames, self.patch_area_pixels)) # The off axis PSF at each point
normalised_psf, norm, fwhm = normalize_psf(self.psf,
fwhm='fit',
size=None,
threshold=None,
mask_core=None,
model='airy',
imlib='vip-fft',
interpolation='lanczos4',
force_odd=False,
full_output=True,
verbose=self.verbose,
debug=False)

psf_mask = create_boolean_circular_mask(
normalised_psf.shape, radius=self.fwhm)
# The off axis PSF at each point
hoff = np.zeros(
(self.num_frames, self.num_frames, self.patch_area_pixels))
# Create arrays needed for storage
# Store for each image pixel, for each temporal frame an image
# for patches: for each time, we need to store a column of patches
Expand All @@ -611,8 +614,9 @@ def flux_estimate(self,
ests = []
stds = []
for i, p0 in enumerate(phi0s):
p0 = (p0[1],p0[0])
angles_px = np.array(get_rotated_pixel_coords(x, y, p0, self.angles))
p0 = (p0[1], p0[0])
angles_px = np.array(
get_rotated_pixel_coords(x, y, p0, self.angles))
hon = []
for l, ang in enumerate(angles_px):
# Get the column of patches at this point
Expand All @@ -622,10 +626,11 @@ def flux_estimate(self,
imlib='vip-fft',
interpolation='lanczos4',
border_mode='reflect')[psf_mask]
hoff[l,l] = offax
hoff[l, l] = offax
hon.append(offax)

Cinv, m, patches = self.compute_statistics(np.array(angles_px).astype(int))
Cinv, m, patches = self.compute_statistics(
np.array(angles_px).astype(int))
# Get Angles
# Ensure within image bounds
# Extract relevant patches and statistics
Expand All @@ -635,28 +640,29 @@ def flux_estimate(self,
for l, ang in enumerate(angles_px):
Cinlst.append(Cinv[int(ang[0]), int(ang[1])])
mlst.append(m[int(ang[0]), int(ang[1])])
patch.append(patches[int(ang[0]), int(ang[1]),l])
a = self.al(hon,Cinlst)
b = self.bl(hon, Cinlst,patch, mlst)
patch.append(patches[int(ang[0]), int(ang[1]), l])
a = self.al(hon, Cinlst)
b = self.bl(hon, Cinlst, patch, mlst)
print(b/a)
# Fill patches and signal template


# Unbiased flux estimation
ahat = initial_est[i]
aprev = 1e10 # Arbitrary large value so that the loop will run
aprev = 1e10 # Arbitrary large value so that the loop will run
while np.abs(ahat - aprev) > np.abs(ahat * eps):
a = 0.0
b = 0.0
# the mean of a temporal column of patches at each pixel
m = np.zeros((self.num_frames, self.patch_area_pixels))
# the inverse covariance matrix at each point
Cinv = np.zeros((self.num_frames, self.patch_area_pixels, self.patch_area_pixels))
Cinv = np.zeros(
(self.num_frames, self.patch_area_pixels, self.patch_area_pixels))

# Patches here are columns in time
for l,ang in enumerate(angles_px):
for l, ang in enumerate(angles_px):
apatch = self.get_patch(ang.astype(int))
m[l], Cinv[l] = self.iterate_flux_calc(ahat, apatch, hoff[l])
m[l], Cinv[l] = self.iterate_flux_calc(
ahat, apatch, hoff[l])
# Patches here are where the planet is expected to be
a = self.al(hon, Cinv)
b = self.bl(hon, Cinv, patch, m)
Expand Down Expand Up @@ -698,8 +704,9 @@ def iterate_flux_calc(self, est: float, patch: np.ndarray,
if patch is None:
return None, None

unbiased = np.array([apatch - est*model[l] for l,apatch in enumerate(patch)])
m,Cinv = compute_statistics_at_pixel(unbiased)
unbiased = np.array([apatch - est*model[l]
for l, apatch in enumerate(patch)])
m, Cinv = compute_statistics_at_pixel(unbiased)
return m, Cinv

def subpixel_threshold_detect(self, snr_map: np.ndarray,
Expand All @@ -719,7 +726,7 @@ def subpixel_threshold_detect(self, snr_map: np.ndarray,
be defined as a region of an image in which some properties are constant
or vary within a prescribed range of values. See ``Notes`` below to read
about the algorithm details.
Parameters
----------
snr_map : numpy ndarray, 2d
Expand All @@ -746,7 +753,7 @@ def subpixel_threshold_detect(self, snr_map: np.ndarray,
The number of processes for running the ``snrmap`` function.
verbose : bool, optional
Whether to print to stdout information about found blobs.
Returns
-------
peaks : np.ndarray
Expand All @@ -764,7 +771,7 @@ def subpixel_threshold_detect(self, snr_map: np.ndarray,
nproc=cpu,
plot=False,
debug=False,
full_output=full_output,
full_output=True,
verbose=self.verbose)
return peaks.T

Expand Down Expand Up @@ -802,7 +809,7 @@ def pixel_threshold_detection(
y.append(y_center)
return np.array(list(zip(x, y)))

def compute_statistics(self, phi0s : np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
def compute_statistics(self, phi0s: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
This function computes the mean and inverse covariance matrix for
each patch in the image stack in Serial. Used by FastPACO and flux
Expand Down Expand Up @@ -832,23 +839,27 @@ def compute_statistics(self, phi0s : np.ndarray) -> Tuple[np.ndarray, np.ndarray

# Store for each image pixel, for each temporal frame an image
# for patches: for each time, we need to store a column of patches
patch = np.zeros((self.width, self.height, self.num_frames, self.patch_area_pixels))
patch = np.zeros(
(self.width, self.height, self.num_frames, self.patch_area_pixels))

# the mean of a temporal column of patches centered at each pixel
m = np.zeros((self.height, self.width, self.patch_area_pixels))
# the inverse covariance matrix at each point
Cinv = np.zeros((self.height, self.width, self.patch_area_pixels, self.patch_area_pixels))
Cinv = np.zeros(
(self.height, self.width, self.patch_area_pixels, self.patch_area_pixels))

# *** SERIAL ***
# Loop over all pixels
# i is the same as theta_k in the PACO paper
for p0 in phi0s:
apatch = self.get_patch(p0)
# For some black magic reason this needs to be inverted here.
m[p0[1]][p0[0]], Cinv[p0[1]][p0[0]] = compute_statistics_at_pixel(apatch)
m[p0[1]][p0[0]], Cinv[p0[1]][p0[0]
] = compute_statistics_at_pixel(apatch)
patch[p0[1]][p0[0]] = apatch
return Cinv, m, patch


"""
**************************************************
* *
Expand All @@ -864,9 +875,9 @@ class FastPACO(PACO):
"""

def PACOCalc(self,
phi0s : np.ndarray,
use_subpixel_psf_astrometry : Optional[bool] = True,
cpu : Optional[int] = 1) -> None:
phi0s: np.ndarray,
use_subpixel_psf_astrometry: Optional[bool] = True,
cpu: Optional[int] = 1) -> None:
"""
PACOCalc
Expand Down Expand Up @@ -899,7 +910,7 @@ def PACOCalc(self,

a = np.zeros(npx) # Setup output arrays
b = np.zeros(npx)
phi0s = np.array([phi0s[:,1],phi0s[:,0]]).T
phi0s = np.array([phi0s[:, 1], phi0s[:, 0]]).T

if cpu == 1:
Cinv, m, patches = self.compute_statistics(phi0s)
Expand All @@ -917,7 +928,8 @@ def PACOCalc(self,
full_output=False,
verbose=self.verbose,
debug=False)
psf_mask = create_boolean_circular_mask(normalised_psf.shape, radius=self.fwhm)
psf_mask = create_boolean_circular_mask(
normalised_psf.shape, radius=self.fwhm)

# Create arrays needed for storage
# Store for each image pixel, for each temporal frame an image
Expand Down Expand Up @@ -999,13 +1011,13 @@ def compute_statistics_parallel(self,
"""

if self.verbose:
print("Precomputing Statistics using %d Processes..."%cpu)
print("Precomputing Statistics using %d Processes..." % cpu)
# the mean of a temporal column of patches at each pixel
m = np.zeros((self.height*self.width*self.patch_area_pixels))
# the inverse covariance matrix at each point
Cinv = np.zeros((self.height*
self.width*
self.patch_area_pixels*
Cinv = np.zeros((self.height *
self.width *
self.patch_area_pixels *
self.patch_area_pixels))

# *** Parallel Processing ***
Expand Down Expand Up @@ -1046,9 +1058,9 @@ def compute_statistics_parallel(self,
self.width,
self.patch_area_pixels,
self.patch_area_pixels))
patches = np.swapaxes(patches,0,1)
m = np.swapaxes(m,0,1)
Cinv = np.swapaxes(Cinv,0,1)
patches = np.swapaxes(patches, 0, 1)
m = np.swapaxes(m, 0, 1)
Cinv = np.swapaxes(Cinv, 0, 1)

return Cinv, m, patches

Expand All @@ -1068,9 +1080,9 @@ class FullPACO(PACO):
"""

def PACOCalc(self,
phi0s : np.ndarray,
use_subpixel_psf_astrometry : Optional[bool] = True,
cpu : Optional[int] = 1) -> None:
phi0s: np.ndarray,
use_subpixel_psf_astrometry: Optional[bool] = True,
cpu: Optional[int] = 1) -> None:
"""
PACOCalc
Expand Down Expand Up @@ -1116,7 +1128,8 @@ def PACOCalc(self,
full_output=False,
verbose=self.verbose,
debug=False)
psf_mask = create_boolean_circular_mask(normalised_psf.shape, radius=self.fwhm)
psf_mask = create_boolean_circular_mask(
normalised_psf.shape, radius=self.fwhm)

if self.verbose:
print("Running Full PACO...")
Expand All @@ -1128,17 +1141,20 @@ def PACOCalc(self,
print("Multiprocessing for full PACO is not yet implemented!")

# Store intermediate results
patch = np.zeros((self.width, self.height, self.num_frames, self.patch_area_pixels))
patch = np.zeros(
(self.width, self.height, self.num_frames, self.patch_area_pixels))
# the mean of a temporal column of patches centered at each pixel
m = np.zeros((self.height, self.width, self.patch_area_pixels))
# the inverse covariance matrix at each point
Cinv = np.zeros((self.height, self.width, self.patch_area_pixels, self.patch_area_pixels))
Cinv = np.zeros(
(self.height, self.width, self.patch_area_pixels, self.patch_area_pixels))

# Loop over all pixels
# i is the same as theta_k in the PACO paper
for i, p0 in enumerate(phi0s):
# Get list of pixels for each rotation angle
angles_px = get_rotated_pixel_coords(x, y, (p0[1],p0[0]), self.angles)
angles_px = get_rotated_pixel_coords(
x, y, (p0[1], p0[0]), self.angles)

# Ensure within image bounds
if(int(np.max(angles_px.flatten())) >= self.width or
Expand All @@ -1155,16 +1171,17 @@ def PACOCalc(self,
clst = []
for l, ang in enumerate(angles_px):
# Get the column of patches at this point
if np.max(patch[int(ang[0]),int(ang[1])]) == 0:
apatch = self.get_patch((int(ang[1]),int(ang[0])))
patch[int(ang[0]),int(ang[1])] = apatch
m[int(ang[0]),int(ang[1])], Cinv[int(ang[0]),int(ang[1])] = compute_statistics_at_pixel(apatch)
if np.max(patch[int(ang[0]), int(ang[1])]) == 0:
apatch = self.get_patch((int(ang[1]), int(ang[0])))
patch[int(ang[0]), int(ang[1])] = apatch
m[int(ang[0]), int(ang[1])], Cinv[int(ang[0]), int(
ang[1])] = compute_statistics_at_pixel(apatch)
else:
apatch = patch[int(ang[0]),int(ang[1])]
apatch = patch[int(ang[0]), int(ang[1])]
if apatch is None:
continue
mlst.append(m[int(ang[0]),int(ang[1])])
clst.append(Cinv[int(ang[0]),int(ang[1])])
mlst.append(m[int(ang[0]), int(ang[1])])
clst.append(Cinv[int(ang[0]), int(ang[1])])

current_patch.append(apatch)
if use_subpixel_psf_astrometry:
Expand All @@ -1179,7 +1196,8 @@ def PACOCalc(self,

h.append(offax)
current_patch = np.array(current_patch)
patches = np.array([current_patch[l,l] for l in range(len(angles_px))])
patches = np.array([current_patch[l, l]
for l in range(len(angles_px))])
h = np.array(h)
mlst = np.array(mlst)
clst = np.array(clst)
Expand Down
Loading

0 comments on commit 0ac5797

Please sign in to comment.