如何对两个数组进行排序以获得平滑度

Say I have two arrays of the same size A and B. For definiteness lets assume they are two dimensional. The values in the two arrays represent some data which should smoothly depend on the position in the array. However some of the values in A have been swapped with their corresponding value in B. I would like to reverse this swapping, but am struggling to find a criterion to tell me when two values should be swapped.

假设我有两个相同大小的A和B阵列。对于明确性,我们假设它们是二维的。两个数组中的值表示一些数据,这些数据应平滑地取决于数组中的位置。但是A中的一些值已经与B中的相应值交换了。我想反转这个交换,但是我很难找到一个标准来告诉我什么时候应该交换两个值。

An example should help. Here is some python code that creates two such arrays and randomly swaps some of their elements

一个例子应该有帮助。下面是一些python代码,它创建了两个这样的数组,并随机交换了一些元素

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

### create meshgrid ###
x = np.linspace(-10,10,15);
y = np.linspace(-10,10,11);

[X,Y] = np.meshgrid(x,y);

### two sufficiently smooth functions on the meshgrid ###
A = -X**2-Y**2;
B = X**2-Y**2-100;

### plot ###
ax=plt.subplot(2,2,1)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax.set_title('A')
ax2=plt.subplot(2,2,2)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])
ax2.set_title('B')

### randomly exchange a few of the elements of A and B ###
for i in np.arange(0,15):
    for j in np.arange(0,11):
        randNumb = random.random();
        if randNumb>0.8:
            mem=A[j,i];
            A[j,i] = B[j,i];
            B[j,i] = mem;

### plot for comparison ###
ax=plt.subplot(2,2,3)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax2=plt.subplot(2,2,4)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])

plt.show()

This results in the following plot:

这导致以下情节:

The upper two plots are the original arrays A and B, the lower two are the shuffled versions. The exercise is now to reverse this process, i.e. to create the original arrays A and B from the shuffled versions.

上面两个图是原始阵列A和B,下面两个是洗牌版本。现在,练习是为了逆转这个过程,即从混洗版本中创建原始阵列A和B.

A note about what I mean by 'smooth'. Of course the algorithm will only work if the original data is actually sufficiently smooth, meaning that neighbouring data points in one array are sufficiently close in value for all point. A solution may assume that this is the case.

关于'平滑'的意思。当然,算法仅在原始数据实际上足够平滑时才起作用,这意味着一个阵列中的相邻数据点对于所有点都足够接近。解决方案可以假设是这种情况。

Also note that this exercise is super easy to do by eye in the example above. The problem I'm having in implementing this is finding a good criterion to tell me when to swap to cells.

另请注意,在上面的示例中,此练习非常容易通过眼睛进行。我在实现这个问题时遇到的问题是找到一个很好的标准来告诉我何时交换到单元格。

Note that the reverse transformation is allowed to relabel A and B of course.

注意,允许反向变换重新标记A和B.

2 个解决方案

#1

One robust method is to compare the mean of the 4 neighbor pixel values to the value actually present in each pixel. That is, for each pixel, compute the mean of the 4 neighbor pixels in both A and B, and compare these to the actual value of this pixel in both A and B. The following condition works nicely, and is really a sort of least-squares method:

一种稳健的方法是将4个相邻像素值的平均值与每个像素中实际存在的值进行比较。也就是说,对于每个像素,计算A和B中4个相邻像素的平均值,并将它们与A和B中该像素的实际值进行比较。以下条件很好地工作,并且实际上是一种最小的-squares方法:

if (  (A[i, j] - A_mean)**2 + (B[i, j] - B_mean)**2
    > (A[i, j] - B_mean)**2 + (B[i, j] - A_mean)**2
    ):
    # Do swap

Here, A_mean and B_mean are the mean of the 4 neighbor pixels.

这里,A_mean和B_mean是4个相邻像素的平均值。

Another important thing to consider is the fact that one sweep through all pixels is not necessarily enough: It might happen that after one sweep, the correction swaps has made it possible for the above condition to recognize more pixels which ought to be swapped. Thus, we have to sweep over the arrays until a "steady state" has been found.

需要考虑的另一个重要事项是,扫描所有像素并不一定是足够的:可能发生在一次扫描之后,校正交换使得上述条件可以识别应该交换的更多像素。因此,我们必须扫描阵列,直到找到“稳定状态”。

Complete code

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

### create meshgrid ###
x = np.linspace(-10,10,15);
y = np.linspace(-10,10,11);

[X,Y] = np.meshgrid(x,y);

### two sufficiently smooth functions on the meshgrid ###
A = -X**2-Y**2;
B = X**2-Y**2-100;

### plot ###
ax=plt.subplot(3,2,1)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax.set_title('A')
ax2=plt.subplot(3,2,2)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])
ax2.set_title('B')

### randomly exchange a few of the elements of A and B ###
for i in np.arange(0,15):
    for j in np.arange(0,11):
        randNumb = random.random();
        if randNumb>0.8:
            mem=A[j,i];
            A[j,i] = B[j,i];
            B[j,i] = mem;

### plot for comparison ###
ax=plt.subplot(3,2,3)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax2=plt.subplot(3,2,4)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])

### swap back ###
N, M = A.shape
swapped = True
while swapped:
    swapped = False
    for i in range(N):
        for j in range(M):
            A_mean = np.mean([A[i - 1    , j - 1    ],
                              A[i - 1    , (j + 1)%M],
                              A[(i + 1)%N, j - 1    ],
                              A[(i + 1)%N, (j + 1)%M],
                              ])
            B_mean = np.mean([B[i - 1    , j - 1    ],
                              B[i - 1    , (j + 1)%M],
                              B[(i + 1)%N, j - 1    ],
                              B[(i + 1)%N, (j + 1)%M],
                              ])
            if (  (A[i, j] - A_mean)**2 + (B[i, j] - B_mean)**2
                > (A[i, j] - B_mean)**2 + (B[i, j] - A_mean)**2
                ):
                # Do swap
                A[i, j], B[i, j] = B[i, j], A[i, j]
                swapped = True

### plot result ###
ax=plt.subplot(3,2,5)
im1=ax.imshow(A,extent=[-10, 10, -10, 10])
ax2=plt.subplot(3,2,6)
im2=ax2.imshow(B,extent=[-10, 10, -10, 10])

plt.show()

Note that the above code consider the arrays to be periodic, in the sense that the neighbor pixels of pixels at the boundary are chosen as those on the opposite boundary (which is the case for the arrays you provided in the example). This index-wrapping happens automatically when the index becomes negative, but not when the index becomes larger than or equal to the given dimension of the array, which is why the modulo operator % is used.

请注意,上面的代码认为数组是周期性的,因为边界处像素的相邻像素被选择为相对边界上的像素(这是您在示例中提供的数组的情况)。当索引变为负数时,此索引包装会自动发生,但当索引变得大于或等于数组的给定维数时,则不会发生这种情况,这就是使用模运算符%的原因。

As a bonus trick, notice how I swap A[i, j] and B[i, j] without the need of the temporary mem variable. Also, my outer loop is over the first dimension (the one with length 11), while my inner loop is over the second dimension (the one with length 15). This is faster than doing the reverse loop order, since now each element is visited in contiguous order (the order in which they actually exist in memory).

作为一个额外的技巧,请注意我如何交换A [i,j]和B [i,j]而不需要临时mem变量。另外,我的外环在第一维(长度为11的那个)上,而我的内环在第二维(长度为15的那个)上。这比执行反向循环顺序更快,因为现在每个元素都以连续的顺序访问(它们实际存在于内存中的顺序)。

Finally, note that this method is not guaranteed to always work. It may happen that so many nearby pixels are swapped that the "correct" solution cannot be found. This will however be the case regardless of what criterion you choose to determine whether two pixels should be swapped or not.

最后,请注意,此方法不能保证始终有效。可能会发生交换这么多邻近像素的情况,无法找到“正确”的解决方案。然而,无论您选择何种标准来确定是否应交换两个像素,情况都是如此。

Edit (non-periodic arrays)

For non-periodic arrays, the boundary pixels will have fewer than 4 neighbors (3 for edge pixels, 2 for corner pixels). Something along these lines:

对于非周期性阵列,边界像素将具有少于4个邻居(边缘像素为3,角落像素为2)。这些方面的东西:

A_neighbors = []
B_neighbors = []
if i > 0 and j > 0:
    A_neighbors.append(A[i - 1, j - 1])
    B_neighbors.append(B[i - 1, j - 1])
if i > 0 and j < M - 1:
    A_neighbors.append(A[i - 1, j + 1])
    B_neighbors.append(B[i - 1, j + 1])
if i < N - 1 and j > 0:
    A_neighbors.append(A[i + 1, j - 1])
    B_neighbors.append(B[i + 1, j - 1])
if i < N - 1 and j < M - 1:
    A_neighbors.append(A[i + 1, j + 1])
    B_neighbors.append(B[i + 1, j + 1])
A_mean = np.mean(A_neighbors)
B_mean = np.mean(B_neighbors)

Note that with fewer neighbors, the method becomes less robust. You could also experiment with using the nearest 8 pixels (that is, include diagonal neighbors) rather than just 4.

请注意,使用较少的邻居,该方法变得不那么健壮。您还可以尝试使用最近的8个像素(即包括对角线邻居)而不是4个像素。

#2

Note: You might want to restate your question on MathOverflow.

注意:您可能希望在MathOverflow上重述您的问题。

First, iterate over all horizontal and vertical neighbours. (It's up to you to also consider diagonals or not.)

首先,迭代所有水平和垂直邻居。 (你也可以考虑对角线。)

Then, calculate the difference of all "neighbour" values.

然后,计算所有“邻居”值的差异。

Finally, there are two popular choices:

最后,有两个流行的选择:

• A) Sum over the abs() (absolute value) of all differences

A)对所有差异的abs()(绝对值)求和

• B) Sum over the square of all differences

B)所有差异的平方和

Your optimization goal is to minize that sum.

您的优化目标是缩小总和。

Variant A) is usually more intuitive, while variant B) is usually easier to track with optimization tools.

变体A)通常更直观,而变体B)通常更容易使用优化工具进行跟踪。

(The problem with A) is that the function is smooth, but not differentiable, while B) is smooth and differentiable, i.e. B) "behaves" better under mathematical analysis.)

(A的问题)是函数是平滑的,但不是可微分的,而B)是平滑的和可微分的,即B)在数学分析下“表现得更好”。)