Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dev space resample #74

Open
wants to merge 29 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 19 commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
16ca19f
Comments and todo
oscarbates Oct 5, 2023
260068e
Anti-aliasing filter implementation
oscarbates Oct 5, 2023
a255c51
Autodetect mute and hann filter
oscarbates Oct 6, 2023
802c43c
resample pipline with filters
oscarbates Oct 10, 2023
a940c62
Identified an issue in the resampy library
oscarbates Oct 13, 2023
851bb43
Clean-up
oscarbates Oct 13, 2023
14cdac1
Merge branch 'master' into upd_time_resample
oscarbates Nov 12, 2023
67d5999
Filters and mutes
oscarbates Nov 21, 2023
db7e42f
Automatic muting diabled for now
oscarbates Nov 21, 2023
b7d5d83
Removing redundant breakpoint
oscarbates Nov 22, 2023
d326733
space_resample method into problem
oscarbates Nov 23, 2023
a78cede
Correct resample call, more helpful error message
oscarbates Dec 8, 2023
81b7295
Shape rounding in space resmaple with comment
oscarbates Dec 8, 2023
b4b584f
If mute does not exist, skip
oscarbates Dec 8, 2023
ca70009
Merge branch 'master' into dev_space_resample
oscarbates Dec 13, 2023
a6c49a9
Time resampling update
oscarbates Dec 15, 2023
fadd1bd
Tweaking space resampling
oscarbates Dec 18, 2023
8b99c60
Removing breakpoint
oscarbates Dec 18, 2023
9bc1304
Fixed evil shape bug
oscarbates Dec 18, 2023
9422582
Resolving comments in review
oscarbates Dec 19, 2023
2c95252
Flake8 corrections
oscarbates Dec 19, 2023
98484a0
Merge branch 'master' into dev_space_resample
oscarbates Dec 19, 2023
f91791b
Bug fix
oscarbates Dec 19, 2023
17fe08a
New new_shape calculation, now match ndimage zoom
oscarbates Dec 20, 2023
b21b5b2
Flake8 fix
oscarbates Dec 20, 2023
0fe16ed
Merge branch 'master' into dev_space_resample
oscarbates Feb 19, 2024
de3ac47
Changes for PR
oscarbates Mar 15, 2024
578c218
minor fix and flake8 corrections
oscarbates Mar 15, 2024
5f66d23
Default behaviour space resample
oscarbates Mar 15, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion stride/optimisation/pipelines/steps/filter_traces.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,9 @@ def _apply(self, traces, **kwargs):
filter_type = kwargs.pop('filter_type', self.filter_type or default_filter_type)

method_name = '%s_filter_%s' % (pass_type, filter_type)
method = getattr(filters, method_name)
method = getattr(filters, method_name, None)
if method is None:
raise Exception('Requested filter does not exist. Implemented filters are butterworth, fir & cos.')

filtered = method(traces.extended_data, *args, zero_phase=False, **kwargs)

Expand Down
18 changes: 12 additions & 6 deletions stride/problem/acquisitions.py
Original file line number Diff line number Diff line change
Expand Up @@ -121,9 +121,9 @@ def __init__(self, id, name=None, problem=None, **kwargs):
for receiver in receivers:
self._receivers[receiver.id] = receiver

self.wavelets = Traces(name='wavelets', transducer_ids=self.source_ids, grid=self.grid)
self.observed = Traces(name='observed', transducer_ids=self.receiver_ids, grid=self.grid)
self.delays = Traces(name='delays', transducer_ids=self.source_ids, shape=(len(sources), 1), grid=self.grid)
self.wavelets = self._traces(name='wavelets', transducer_ids=self.source_ids, grid=self.grid)
self.observed = self._traces(name='observed', transducer_ids=self.receiver_ids, grid=self.grid)
self.delays = self._traces(name='delays', transducer_ids=self.source_ids, shape=(len(sources), 1), grid=self.grid)

if delays is not None:
self.delays.data[:, 0] = delays
Expand Down Expand Up @@ -402,6 +402,12 @@ def append_observed(self, *args, **kwargs):
else:
self._acquisitions.dump(*args, **kwargs)

def _traces(self, name, transducer_ids, grid, shape=None):
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
if shape == None:
return Traces(name=name, transducer_ids=transducer_ids, grid=grid)
else:
return Traces(name=name, transducer_ids=transducer_ids, shape=shape, grid=grid)

def __get_desc__(self, **kwargs):
description = {
'id': self.id,
Expand Down Expand Up @@ -439,15 +445,15 @@ def __set_desc__(self, description):
receiver = None
self._receivers[receiver_id] = receiver

self.wavelets = Traces(name='wavelets', transducer_ids=self.source_ids, grid=self.grid)
self.wavelets = self._traces(name='wavelets', transducer_ids=self.source_ids, grid=self.grid)
if 'wavelets' in description:
self.wavelets.__set_desc__(description.wavelets)

self.observed = Traces(name='observed', transducer_ids=self.receiver_ids, grid=self.grid)
self.observed = self._traces(name='observed', transducer_ids=self.receiver_ids, grid=self.grid)
if 'observed' in description:
self.observed.__set_desc__(description.observed)

self.delays = Traces(name='delays', transducer_ids=self.source_ids, shape=(len(self.source_ids), 1), grid=self.grid)
self.delays = self._traces(name='delays', transducer_ids=self.source_ids, shape=(len(self.source_ids), 1), grid=self.grid)
if 'delays' in description:
self.delays.__set_desc__(description.delays)

Expand Down
2 changes: 1 addition & 1 deletion stride/problem/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ def slow_time(self):
return self._grid.slow_time

def resample(self, grid=None, space=None, time=None, slow_time=None):
raise NotImplementedError('Resampling has not been implemented yet.')
raise NotImplementedError('Resampling has not been implemented in this class yet. Alternatively, try to access space_resample or time_resample via problem.')


class Saved:
Expand Down
121 changes: 48 additions & 73 deletions stride/problem/data.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
from .base import GriddedSaved
from ..core import Variable
from .. import plotting
from ..utils import filters


__all__ = ['Data', 'StructuredData', 'Scalar', 'ScalarField', 'VectorField', 'Traces',
Expand Down Expand Up @@ -908,27 +909,18 @@ def stagger_data(self, data, stagger, method='nearest'):

return interp

def resample(self, space=None, **kwargs):
def _resample(self, old_spacing, new_spacing, **kwargs):
"""
Resample the internal (non-padded) data given some new space object.
oscarbates marked this conversation as resolved.
Show resolved Hide resolved

Parameters
----------
space : Space
New space.
old_spacing: float
The old spacing.
new_spacing: float
The new spacing.
order : int, optional
Order of the interplation, default is 3.
prefilter : bool, optional
Determines if the input array is prefiltered
before interpolation. If downsampling, this defaults to ``False`` as an anti-aliasing filter
will be applied instead. If upsampling, this defaults to ``True``.
anti_alias : bool, optional
Whether a Gaussian filter is applied to smooth the data before interpolation.
The default is ``True``. This is only applied when downsampling.
anti_alias_sigma : float or tuple of floats, optional
Gaussian filter standard deviations used for the anti-aliasing filter.
The default is (d - 1) / 2 where d is the downsampling factor and d > 1. When upsampling,
d < 1, and no anti-aliasing filter is applied.

Returns
-------
Expand All @@ -940,74 +932,23 @@ def resample(self, space=None, **kwargs):

interp = []
for t in range(data.shape[0]):
interp.append(self.resample_data(data[t], space, **kwargs))
interp.append(self._resample_data(data[t], old_spacing, new_spacing, **kwargs))

interp = np.stack(interp, axis=0)

else:
interp = self.resample_data(self.data, space, **kwargs)
interp = self._resample_data(self.data, old_spacing, new_spacing, **kwargs)

self.grid.space = space
self._init_shape()
self._init_shape() # NOTE this is an in-place operation, unlike time resample
self._data = self.pad_data(interp)

def resample_data(self, data, space, **kwargs):
"""
Resample the data given some new space object.

Parameters
----------
data : ndarray
Data to stagger.
space : Space
New space.
order : int, optional
Order of the interpolation, default is 3.
prefilter : bool, optional
Determines if the input array is prefiltered
before interpolation. If downsampling, this defaults to ``False`` as an anti-aliasing filter
will be applied instead. If upsampling, this defaults to ``True``.
anti_alias : bool, optional
Whether a Gaussian filter is applied to smooth the data before interpolation.
The default is ``True``. This is only applied when downsampling.
anti_alias_sigma : float or tuple of floats, optional
Gaussian filter standard deviations used for the anti-aliasing filter.
The default is (d - 1) / 2 where d is the downsampling factor and d > 1. When upsampling,
d < 1, and no anti-aliasing filter is applied.
def _resample_data(self, data, old_spacing, new_spacing, **kwargs):

Returns
-------
ndarray
Resampled data.

"""
order = kwargs.pop('order', 3)
prefilter = kwargs.pop('prefilter', True)

resampling_factors = np.array([dx_old/dx_new
for dx_old, dx_new in zip(self.space.spacing, space.spacing)])

# Anti-aliasing is only required for down-sampling interpolation
if any(factor < 1 for factor in resampling_factors):
anti_alias = kwargs.pop('anti_alias', True)

if anti_alias:
anti_alias_sigma = kwargs.pop('anti_alias_sigma', None)

if anti_alias_sigma is not None:
anti_alias_sigma = anti_alias_sigma * np.ones_like(resampling_factors)

if np.any(anti_alias_sigma < 0):
raise ValueError("Anti-alias standard dev. must be equal to or greater than zero")

# Estimate anti-alias standard deviations if none provided
else:
anti_alias_sigma = np.maximum(0, (1/resampling_factors - 1) / 2)

data = scipy.ndimage.gaussian_filter(data, anti_alias_sigma)

# Prefiltering is not necessary if anti-alias filter used
prefilter = False
for dx_old, dx_new in zip(old_spacing, new_spacing)])
ccuetom marked this conversation as resolved.
Show resolved Hide resolved

interp = scipy.ndimage.zoom(data, resampling_factors,
order=order, prefilter=prefilter)
Expand Down Expand Up @@ -1562,13 +1503,47 @@ def plot_one(self, id, **kwargs):

return axis

def _resample(self, factor, new_num, **kwargs):
def _resample(self, old_step, new_step, new_num, freq_niquist, len_hann=None, mute_start=None, **kwargs):
'''
In-place operation to resample a trace to a new time-spacing.
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
Sinc interpolation is used.

Parameters
----------
factor
The ratio dt_old/dt_new.
new_num
The length of the trace.
freq_niquist
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
min(fs_old, fs_new), where fs is the sampling frequency.
len_hann
The length of the Hann window applied around the edge of the mute.
mute_start, default: None
The length of the starting mute in num_samples.
filter_type_antialias, default: cos
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
The method used by the anti-aliasing filter. Options are butterworth, fir & cos.
filter_order_antialias, default: 8
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
The order of the anti-aliasing filter.
oscarbates marked this conversation as resolved.
Show resolved Hide resolved

Returns
-------
'''

filter_type_antialias = kwargs.pop('filter_type_antialias', 'cos')
filter_order_antialias = kwargs.pop('filter_order_antialias', 2)

factor = old_step/new_step
sr_orig = 1
sr_new = factor

if self.allocated:
data = resampy.resample(self.data, sr_orig, sr_new, axis=1) # resample
new_traces = Traces(name=self.name, grid=self.grid, transducer_ids=self._transducer_ids, data=data)
processed_data = self.data
# Resample
processed_data = resampy.resample(processed_data, sr_orig, sr_new, axis=1, parallel=True)

# Fill object
new_traces = Traces(name=self.name, grid=self.grid, transducer_ids=self._transducer_ids, data=processed_data)

else:
new_traces = Traces(name=self.name, grid=self.grid, transducer_ids=self._transducer_ids)

Expand Down
35 changes: 24 additions & 11 deletions stride/problem/domain.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,9 @@ class Space:
"""

def __init__(self, shape=None, spacing=None, extra=None, absorbing=None):
self._set_properties(shape=shape, spacing=spacing, extra=extra, absorbing=absorbing)

def _set_properties(self, shape, spacing, extra, absorbing):
if isinstance(spacing, float):
spacing = (spacing,)*len(shape)

Expand Down Expand Up @@ -74,21 +77,26 @@ def extended_size(self):
"""
return self.extended_limit

def resample(self, spacing, extra=None, absorbing=None):
if isinstance(spacing, float):
spacing = (spacing,)*self.dim
def resample(self, new_spacing, new_extra=None, new_absorbing=None):
'''
Resample space
'''

shape = tuple((np.array(self.size) / np.array(spacing) + 1).astype(int))
if isinstance(new_spacing, float):
new_spacing = (new_spacing,)*self.dim

if extra is None:
extra = tuple((np.array(self.spacing) * (np.array(self.extra) - 1) /
np.array(spacing) + 1).astype(int))
new_shape = tuple((np.round(np.array(self.size) / np.array(new_spacing)) + 1).astype(int))
# TODO there is possibly a rounding error

if absorbing is None:
absorbing = tuple((np.array(self.spacing) * (np.array(self.absorbing) - 1) /
np.array(spacing) + 1).astype(int))
if new_extra is None:
new_extra = tuple((np.round(np.array(self.spacing) * (np.array(self.extra) - 1) /
np.array(new_spacing)) + 1).astype(int))

return Space(shape=shape, spacing=spacing, extra=extra, absorbing=absorbing)
if new_absorbing is None:
new_absorbing = tuple((np.round(np.array(self.spacing) * (np.array(self.absorbing) - 1) /
np.array(new_spacing)) + 1).astype(int))

self._set_properties(shape=new_shape, spacing=new_spacing, extra=new_extra, absorbing=new_absorbing)
oscarbates marked this conversation as resolved.
Show resolved Hide resolved

@property
def inner(self):
Expand Down Expand Up @@ -286,8 +294,13 @@ def resample(self, new_step, new_num):
new_stop = interp_stop

self._set_properties(start=new_start, step=new_step, num=new_num) # Update time

try:
del self.__dict__['grid']
except:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you've creared a _clear_cache method below, why are you not using it?

pass

try:
del self.__dict__['extended_grid']
except:
pass
Expand Down
50 changes: 44 additions & 6 deletions stride/problem/problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -123,18 +123,42 @@ def load(self, *args, **kwargs):
if getattr(problem_property, grid_property) is None:
setattr(problem_property._grid, grid_property, getattr(self, grid_property))

def space_resample(self, new_spacing, new_extra=None, new_absorbing=None, **kwargs):
'''
In-place operation to resample the wavelets and data into a grid with new
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
time-spacing. Sinc interpolation is used.

Parameters
----------
new_spacing : float
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
The space spacing for the interpolated grid.
new_extra : int
The extra grid-points for the interpolated grid. Defaults to rescaling existing extra.
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
new_absorbing : float
The absorbing grid-points for the interpolated grid. Defaults to rescaling existing absorbing.

Returns
-------
'''
old_spacing = self.grid.space.spacing
self.grid.space.resample(new_spacing=new_spacing)
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
new_spacing = self.grid.space.spacing

self.medium.vp._resample(old_spacing, new_spacing, **kwargs) # NOTE this is in-place
oscarbates marked this conversation as resolved.
Show resolved Hide resolved
return self.medium.vp

def time_resample(self, new_step, new_num=None, **kwargs):
'''
In-place operatin to resample the wavelets and data into a grid with new
In-place operation to resample the wavelets and data into a grid with new
time-spacing. Sinc interpolation is used.

Parameters
----------
new_step : float
The time spacing for the interpolated grid
The time spacing for the interpolated grid.
new_num : int, optional
The number of time-points, default is calculated to match input pulse
length in [s]
length in [s].

Returns
-------
Expand All @@ -144,9 +168,23 @@ def time_resample(self, new_step, new_num=None, **kwargs):
old_num = self.time.num
self.grid.time.resample(new_step=new_step, new_num=new_num)
oscarbates marked this conversation as resolved.
Show resolved Hide resolved

for shot in self.acquisitions.shots:
shot.wavelets = shot.wavelets._resample(factor=old_step/new_step, new_num=new_num) # resample wavelet
shot.observed = shot.observed._resample(factor=old_step/new_step, new_num=new_num) # resample observed
freq_niquist_hz = min(1/old_step, 1/new_step) # anti-aliasing filter using max sampling frequency
oscarbates marked this conversation as resolved.
Show resolved Hide resolved

for shot in self.acquisitions.shots: # TODO, implement low-pass filtering with freq_max as pass band
oscarbates marked this conversation as resolved.
Show resolved Hide resolved

shot.wavelets = shot.wavelets._resample( # resample wavelet
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if resample is an in-place operation then you don't need to assign the result, only call the method

old_step=old_step,
new_step=new_step,
new_num=new_num,
freq_niquist=freq_niquist_hz*old_step, # Convert [Hz] to dimensionless frequency
**kwargs)

shot.observed = shot.observed._resample( # resample observed
old_step=old_step,
new_step=new_step,
new_num=new_num,
freq_niquist=freq_niquist_hz*old_step,
**kwargs)

def dump(self, *args, **kwargs):
"""
Expand Down
Loading
Loading