OttoStruve / muler

A Python package for working with pipeline-produced spectra from IGRINS, HPF, and Keck NIRSPEC
https://muler.readthedocs.io
MIT License
15 stars 9 forks source link

Add a cross correlate method #50

Open gully opened 3 years ago

gully commented 3 years ago

Cross correlation is a useful technique for template matching, and can be tricky to setup for new users (how to treat boundaries, resampling, etc). Here we propose to add a .cross_correlate(input_spec) method.

One interesting implementation question is what-object-to-return? A Spectrum1D seems not-quite-right because the cross correlation has a relative and not absolute spectral_axis, and it's not clear that specutils can handle relative coordinates?

gully commented 2 years ago

Specutils already has a cross correlation functionality, so implementhing this should be pretty straightforward.

The hiccup comes with ensuring that the pre-processing steps (continuum normalization, masking) are all completed adequately. But we have some good tutorials on that, so maybe we're not too far off from adding this capability.

kfkaplan commented 1 year ago

I have had success using spicy.ndimage.zoom to enlarge an array using linear interpolation and then scipy.signal.fftconvolve to cross correlate exposures to correct for flexure. It is quite efficient. Obviously, this is different than cross-correlation a template to get RVs but perhaps the code could be adapted for that. I am also not sure how different it would be to the cross-correlation in specutils. Below is a generalized example of my code for shifts in pixel space like so you can get an idea of how it might be implemented.

from scipy.ndimage import zoom
from scipy.signal import fftconvolve
import dumpy as np

#Data and template are two equal sized arrays you are trying to cross-correlate

zoom_amount = 100 #How many times to supersample arrays
maximum_pixel_search = 20 #How far (in pixels) to search for the cross-correlation maximum, used to ignore incorrect far off maxima
zoomed_template = zoom(template, zoom_amount, order=1) #Enlarge arrays to supersample them
zoomed_data = zoom(data, zoom_amount, order=1)
zoomed_template[np.isnan(zoomed_template)] = 0 #Zero out nans
zoomed_data[np.isnan(zoomed_data)] = 0
fft_result = fftconvolve(zoomed_template, zoomed_data[::-1]) #Do FFT cross- correlation
delta_sub_pixels = 0.5*maximum_pixel_search #Apply maximum pixel search cut for finding maxima 
fft_result_length = int(len(fft_result))
x1 = int(fft_result_length/2) - delta_sub_pixels
x2 = int(fft_result_length/2) + delta_sub_pixels
fft_sub_result = fft_result[x1:x2])
fft_sub_result_length =  int(len(fft_sub_result))
find_shift_from_maximum = np.unravel_index(np.argmax(fft_sub_result), fft_sub_result_length) #Find maximum
dx = (find_shift_from_maximum - (fft_sub_result_length/2)) / zoom_amount #This is the actual shift in pixels
gully commented 1 year ago

@kfkaplan I like that use of spicy.ndimage.zoom, sweet!

My inclination is to adopt the specutils technique when available, which makes the code easier to maintain in the long term. Have you taken a look at the specutils implementation? Your code could potentially be incorporated in the upstream specutils codebase if it were shown to be much more performant than their implementation.