Consider / design "Rebin Strategy" Interface

[calcuttj]

We should have some interface that allows us to be consistent with Rebinning.

I.e. if we have to downsample at some point and then upsample later, we want to make sure that happens consistently and doesn't need to be configured in multiple places

Need to develop thoughts/sketch this out

[brettviren]

Here are my (somewhat exhaustive) thoughts on this.

Prior work

WCT aux has a DFP graph node to resample IFrame

https://github.com/WireCell/wire-cell-toolkit/blob/master/aux/inc/WireCellAux/Resampler.h

This (non torch) WCT node applies LMN resampling with code directly implemented in the node.

I like the idea to factor out the core resampling algorithm into a "service" type interface.

Proposal for an ITorchResample interface class

Omitting some needed boilerplate I suggest this interface class:

class ITorchResample {
public:
    virtual torch::Tensor resample(torch::Tensor ten,
                                   double input_period, 
                                   double output_period, 
                                   size_t axis=-1) const = 0;
};

The axis determines which axis of tensor ten to resample. All other axes are considered as "batch". Default of axis=-1 would be the time dimension using SPNG's convention. The last section describes how this single-axis interface limits possible resampling methods.

An implementation of this interface determines whether it does upsample or downsample by the ratio of periods.

At first, I thought we needed separate upsample/downsample interfaces thinking some resampling methods have no clear bidirectionality. However, I am not able to think of any. I did some chatting with Gemini and while it gave some examples none were applicable.

If some future algorithm lacks "one way" then it can raise an exception when given periods that indicate going "the other way".

Here are the list of 3 resampling algorithms that I expect may be of interest to implement:

Interval-domain interpolation

A linear or cubic interpolation can be formed and resampling is then done simply by evaluating the interpolation at the new period (old/new ratio).

Rebin and unrebin

Rebin downsampling works only if D = output_period/input_period is integer. It combines D consecutive input samples with a configured an operations (eg "add" or "maxpool" or "decimate") determine how to combine D samples including a normalization.

Unrebin upsampling is the "moral inverse" of "rebin" and only for U = input_period/output_period integer. Simplest implementation is to replicate each input sample to U consecutive output samples and apply a configurable operation to the U output samples ("average", "copy") to determine how to produce the U separate samples.

LMN

LMN works only if R = input_period/output_period is a ratio of integers (or is approximated to one). Upsampling and downsampling are almost identical except up adds frequency bins, down removes them.

Configuration issues

Each implementation needs some kind of normalization directives that are expressed in an algorithm-specific way. Eg, rebin/unrebin needs the "operators". LMN takes an "sample interpretation" enum (eg instantaneous vs integral) and produces corresponding normalization.

Some related nodes/services

A few operations are related to resampling in general and are best implemented separately and not combined with any given resampling operation itself. Ideally, we defined "service" type interfaces for these to be used from a DFP node (possibly combined with the resampling "service").

Truncating/padding

Some methods can require interval-domain truncating or padding in order to hit either the integer ratio (rebin/unrebin) or ratio-of-integer (LMN) requirement.

A good implementation will configure the truncation/padding target in terms of this requirement. A less general one might be configured with a hard-wired size.

Both truncating and padding requires a choice to operate at the start of the tensor, the end or both.

Padding requires a choice for the values in the padded samples. Reasonable options include "zero", some scalar value or some decay (linear, exponential, sinusoidal). The "right" choice depends on the resampling method and the meaning of the tensor.

Filtering

Both downsampling and upsampling may required filters to combat artifacts.

Generally, downsampling will lead to aliasing when there is any power above the new (lower) Nyquist frequency. A low-pass ("antialiasing") filter MUST be applied prior to downsampling to avoid this. This is of particular problem for interpolation and rebin downsampling. LMN downsampling applies a "perfect" antialiasing filter in that removes the frequency bins above the new Nyquist frequency.

LMN upsampling and downsampling will produce interval-domain high-frequency ringing when the resampled spectrum has a sharp cutoff. This can be removed by a low-pass filter applied either before or after LMN resampling. (Note, for the case of 512ns->500ns LMN upsampling of ADC waveforms, the low-pass filter in SP is more than enough to remove this.)

Tensor data types

Not all resampling algorithms are relevant for all tensor data types. In particular, interpolation (LMN frequency-domain and the interval-domain variants) are not appropriate for resampling "ROI" binary tensors. It is up to the user writing the configuration to choose wisely.

LMN in particular might allow complex input tensors and assume they are already in Fourier space. Or have an option to leave the result in Fourier space. This allows pre/post spectral filters to be applied with fewer FFTs.

Dimensionality caveat

As proposed, bilinear/bicubic interval-domain interpolation resampling is excluded as it requires 2 axes and the input/output periods must go from scalar to tensor.

The rebin/unrebin algorithm can easily be extended to operate on N-dimensions. These extensions are also excluded from this interface.

Sometimes I do hear people consider "rebin" along the time AND channel dimension. With this current interface, they must do these in a serial manner. This is probably a bit slower than a 2D rebin/unrebin.