Using cross-correlation of spatially shifted data to gain geometric information

Question

Cross-correlation is usually used to estimate e.g. lags between signals that are detected by different sensors at different times. In my case, the situation is slightly different in that I measure signals with sensors at different grid points and I want to recover geometric information about the grid. To be more specific, let's consider the following example.

I have sensors located on a grid $x_1$ from $x = -0.179$ to $x = 0.169$ m, with a spacing of 0.003 m. The sensors basically sample a function that is located somewhere on this grid; let's consider a parabolic function for simplicity, like $f(x) = -2.5 * (x - \delta x)^2 + 0.016$ that is only defined for positive values of $f(x)$, otherwise it is $0$, and $\delta x$ denotes a shift along the x-axis.* I measure this function twice on the grid above, and for the second measurement, the function is mirrored.** However, the issue is that I do not know where exactly it is mirrored and that is, in fact, my goal, i.e. I want to find some $x_0$ that, in terms of the function above, fulfills the condition $\delta x - x_0 \rightarrow -(\delta x - x_0)$, the caveat being that in reality, $\delta x$ is not known, either.

For that purpose, I intend to use the cross correlation of my measured functions $f_1$ and $f_2$, before and after mirroring. For the purposes of this post, let's say that the initial $\delta x$ is $-0.07$ and the (in reality unknown) mirroring position is at $-0.02$ so that the mirrored function's $\delta x$ is at $0.03$. From the step size, we can see that the functions are sampled at different function values, so there is a "non-integer" shift between them. What follows are plots of the functions and their cross correlation. (Cross-correlation is computed using Matlab's xcorr function.)

I can determine the shift between the two functions from the position of the maximum value of the cross correlation function, and interpolation can make this value more accurate. In this case, the maximum value lies at a lag value (or $\tau_{\text{discrete}}$) of -33. Multiplied by the step size I therefore get $-0.099$. This could be improved via interpolation of the cross-correlation function, by which we would estimate $\tau_{\text{grid}}$, and $\tau_{\text{continuous}} = \tau_{\text{discrete}} + \tau_{\text{grid}}$. Since the value is already very close to the true value of $-0.1$, though, I will leave this out for now. If one $\delta x$ was known, the center of rotation $x_0$ could easily be determined by adding half the shift we just calculated to $\delta x$, but in reality, this quantity is an unknown.

The issue is therefore how to get $x_0$. I have had two ideas so far, but both seem not very elegant. The first idea is to get an estimate for $\delta x$; the most straightforward way would be $\text{mean}(x (f(x) > 0))$, i.e. taking the average x-values for all x where the function assumes a value larger than 0. This works, but is not a "nice" solution. The second, very rough idea is as follows: The functions $f_1$ and $f_2$ are symmetric with regards to $x_0$, i.e. $f_1(x - x_0) = f_2(- (x - x0))$, so we interpolate the functions and iteratively try several values for $x_0$. The issue here is on the one hand that I want to avoid loops (but I will do them if necessary) and on the other hand that interpolating functions that include steps is somewhat problematic. I feel like there should be a way to do this with a cross-correlation, which would at least avoid the latter issue. What can I do here?

tl;dr: I want to calculate the center of rotation $x_0$ of a sampled function and its mirrored self using cross-correlation.

*: In reality, the function will have at least two steps (i.e. infinite gradient at that point) and might not be perfectly symmetric. However, the function values close to either end of the grid will be 0, and somewhere on the grid the function continuously assumes values > 0 (i.e. without "holes" where it drops to 0).

**: Technically, I rotate the function by 180° ***, but I am only interested in its projection so we can ignore the changed sign. We could also rotate the grid instead, i.e. $x2 = - \text{flip}(x1)$, and sample the unchanged function on it, which yields the same results.

***: If that helps, we measure the function not just at the angles of rotation 0° and 180°, but also 2° and 182°, 4° and 184°... 178° and 358°.

Extra information: The point $0$ on the grid $x_1$ is a first estimate of the center of rotation $x_0$, so it should be somewhere in the vicinity, probably like $\pm 3 \, \text{cm}$, maybe less.

Answer 1

The easiest seems to me to be just using a weighted sum:

$$ \hat{x}_0 = \frac{\sum f(x_k) x_k}{\sum f(x_k)} $$

My simulation of this and your original data yields:

where the red vertical line indicates the location of the orange (flipped) data and the blue vertical line indicates the location of the blue (original) data.

Then, to estimate the values, just do:

peak = max(data)
plt.plot(x, -2.5*np.power(x-x_hat,2) + peak, 'b')
plt.plot(x, data)
plt.plot(x, -2.5*np.power(x-x_hat_flipped,2) + peak, 'r')
plt.plot(x, np.flip(data))

which yields the following plot:

To get the coefficient, just do something like this:

def plot_error(k, data, x, x_hat):
    error = data-(peak+k*np.power(x-x_hat,2))
    error[data == 0] = 0
    return np.sum(np.power(error,2))

possible_ks = np.linspace(-10,10,10000)
err=np.zeros(len(possible_ks))
ix = 0
for k in possible_ks:
    err[ix] = plot_error(k,data,x,x_hat)
    ix = ix + 1

plt.plot(possible_ks, err)

min_index = np.argmin(err)
print("k is " + str(possible_ks[min_index]))

---

Code Only Below

import numpy as np
import random
import matplotlib.pyplot as plt

def my_function(x, delta_x):
  func = -2.5*np.power(x-delta_x,2)+0.016
  func[func < 0] = 0
  return func

sensor_spacing = 0.003

x = np.arange(-0.179, 0.169 + sensor_spacing, sensor_spacing)


random.seed(84134)
delta_x = 0.011
x_offset_from_zero = random.uniform(0, sensor_spacing)
data = my_function(x + x_offset_from_zero, delta_x)

x_hat = np.sum(x * data)/np.sum(data)
x_hat_flipped = np.sum(x * np.flip(data))/np.sum(data)


plt.plot(x,data)
plt.plot(x,np.flip(data))
plt.plot([x_hat,x_hat], [0,0.018], "b")
plt.plot([x_hat_flipped,x_hat_flipped], [0,0.018], "r")