In [4]:
#everything we will need for this demo
from kbmodpy import kbmod as kb
import numpy as np
import matplotlib.pyplot as plt
import os
path = "../data/demo/"
2D Point Spread Function Array
2D Image array
A Complete image represented as 3 raw_image layers (science, mask, variance)
Stack of layered_images, intended to be the same frame captured at different times
Searches an image_stack for a moving psf
Stores an object's position and motion through an image_stack
In [5]:
p = kb.psf(1.0)
The psf can be cast into a numpy array
In [6]:
np.array(p)
Out[6]:
A psf can also be initialized or set from a numpy array, but the array must be square and have odd dimensions
In [7]:
arr = np.linspace(0.0, 1.0, 9).reshape(3,3)
p2 = kb.psf(arr) # initialized from array
arr = np.square(arr)
p2.set_array(arr) # set from array
np.array(p2)
Out[7]:
There are several methods that get information about its properties
In [8]:
p.get_dim() # dimension of kernel width and height
p.get_radius() # distance from center of kernel to edge
p.get_size() # total number of pixels in the kernel
p.get_sum() # total sum of all pixels in the kernel,
#should be close to 1.0 for a normalized kernel
Out[8]:
In [9]:
#im = kb.layered_image(path+"example.fits")
B. Generate a new image from scratch
In [10]:
im = kb.layered_image("image2", 100, 100, 5.0, 25.0, 0.0)
# name, width, height, background_noise_sigma, variance, capture_time
Artificial objects can easily be added into a layered_image
In [11]:
im.add_object(20.0, 35.0, 2500.0, p)
# x, y, flux, psf
The image pixels can be retrieved as a 2D numpy array
In [17]:
pixels = im.science()
pixels
Out[17]:
The image can mask itself by providing a bitmask of flags (note: masked pixels are set to -9999 so they can be distinguished later from 0.0 pixles)
In [18]:
flags = ~0
flag_exceptions = [32,39]
# mask all of pixels with flags except those with specifed combiniations
im.apply_mask_flags( flags, flag_exceptions )
The image can be convolved with a psf kernel
In [19]:
im.convolve(p)
# note: This function is called interally by stack_search and doesn't need to be
# used directy. It is only exposed because it happens to be a fast
# implementation of a generally useful function
The image at any point can be saved to a file
In [20]:
#im.save_layers(path+"/out") # file will use original name
A layered_image can have its layers set from any numpy array
In [33]:
raw = kb.raw_image( np.ones_like(pixels) )
In [35]:
im.set_science( raw )
im.set_variance( raw )
im.science()
Out[35]:
Get properties
In [16]:
im.get_width()
im.get_height()
im.get_time()
im.get_ppi() # pixels per image, width*height
Out[16]:
In [14]:
count = 10
imlist = [ kb.layered_image("img"+str(n), 50, 50, 10.0, 5.0, n/count) for n in range(count) ]
stack = kb.image_stack( imlist )
# this creates a stack with 10 50x50 images, and times ranging from 0 to 1
Manually set the times the images in the stack were taken
In [15]:
stack.set_times( [0,2,3,4.5,5,6,7,10,11,14] )
A shortcut is provided to initialize a stack automatically from a list of files. If 'MJD' is in the header for each image, the stack will automatically load the times as well. If not, you can set them as above.
In [16]:
import os
files = os.listdir(path)
files = [path+f for f in files if '.fits' in f]
files.sort()
files
Out[16]:
In [17]:
stack = kb.image_stack(files)
A master mask can be generated and applied to the stack
In [18]:
flags = ~0 # mask pixels with any flags
flag_exceptions = [32,39] # unless it has one of these special combinations of flags
master_flags = int('100111', 2) # mask any pixels which have any of
# these flags in more than two images
Most features of the layered_image can be used on the whole stack
In [19]:
stack.apply_mask_flags(flags, flag_exceptions)
stack.apply_master_mask(master_flags, 2)
stack.convolve(p)
stack.get_width()
stack.get_height()
stack.get_ppi()
stack.get_images() # retrieves list of layered_images back from the stack
stack.get_times()
Out[19]:
Here, we will create a very bright object and add it to the images and create a new image stack with the new object.
In [20]:
im_list = stack.get_images()
In [21]:
new_im_list = []
for im, time in zip(im_list, stack.get_times()):
im.add_object(20.0+(time*800.), 35.0+(time*0.), 25000.0, p)
new_im_list.append(im)
In [22]:
stack = kb.image_stack(new_im_list)
In [23]:
search = kb.stack_search( stack, p )
To save psi and images, a directory with "psi" and "phi" folders must be specified.
In [24]:
if os.path.exists(os.path.join(path,'out/psi')) is False:
os.mkdir(os.path.join(path,'out/psi'))
if os.path.exists(os.path.join(path,'out/phi')) is False:
os.mkdir(os.path.join(path,'out/phi'))
search.save_psi_phi(os.path.join(path, 'out'))
Launch a search
In [25]:
#search.gpu(100, 100, 0.2, 0.4, 750, 850, 2)
search.gpu(10, 10, -0.1, 0.1, 750, 850, 2)
# angle_steps, velocity_steps, min_angle, max_angle, min_velocity, max_velocity, min_observations
Save the results to a files
note: format is {x, y, xv, yv, likelihood, flux}
In [26]:
search.save_results(path+"results.txt", 0.05)
# path, fraction of total results to save in file
Trajectories can be retrieved directly from search without writing and reading to file.
However, this is not recommended for a large number of trajectories, as it is not returned as a numpy array, but as a list of the trajectory objects described below
In [27]:
top_results = search.get_results(0, 100)
# start, count
In [28]:
best = top_results[0]
In [29]:
# these numbers are wild because mask flags and search parameters above were chosen randomly
best.flux
best.lh
best.x
best.y
best.x_v
best.y_v
Out[29]:
tests/test_search.py shows a simple example of how to generate a set of images, add an artificial source, and recover it with search
In [30]:
# These top_results are all be duplicating searches on the same bright object we added.
top_results[:20]
Out[30]:
In [ ]: