正如我在评论中所说，医学图像配准是一个有大量研究可用的主题，我不是专家。根据我的阅读，常用的基本思想是定义两个图像之间的映射（在您的情况下是图像及其镜像），然后如果应用映射，则定义平滑度和图像相似性的能量项，最后使用标准（或有时特定于应用程序的）优化技术优化此映射。
我在 Mathematica 中编写了一个快速算法来演示这一点。这不是您应该在医疗应用中使用的算法，只是基本思想的演示。
首先，我加载您的图像，对其进行镜像并将这些图像分割成小块：
src = ColorConvert[Import["http://i.stack.imgur.com/jf709.jpg"],
"Grayscale"];
mirror = ImageReflect[src, Left -> Right];
blockSize = 30;
partsS = ImagePartition[src, {blockSize, blockSize}];
partsM = ImagePartition[mirror, {blockSize, blockSize}];
GraphicsGrid[partsS]
通常，我们会进行近似刚性配准（使用例如关键点或图像时刻），但您的图像几乎居中，所以我会跳过这个。
如果我们看一个块，它是镜像对应的：
{partsS[[6, 10]], partsM[[6, 10]]}
我们可以看到它们很相似，但发生了变化。转变的数量和方向是我们试图找出的。
为了量化匹配相似度，我可以使用平方欧几里得距离：
ListPlot3D[
ImageData[
ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]],
SquaredEuclideanDistance]]]
可悲的是，直接使用这些数据进行优化比我想象的要难，所以我使用了二阶近似值：
fitTerms = {1, x, x^2, y, y^2, x*y};
fit = Fit[
Flatten[MapIndexed[{#2[[1]] - blockSize/2, #2[[2]] -
blockSize/2, #1} &,
ImageData[
ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]],
SquaredEuclideanDistance]], {2}], 1], fitTerms, {x, y}];
Plot3D[fit, {x, -25, 25}, {y, -25, 25}]
该函数与实际的相关函数不同，但对于第一步来说已经足够接近了。让我们为每对块计算这个：
distancesFit = MapThread[
Function[{part, template},
Fit[Flatten[
MapIndexed[{#2[[2]] - blockSize/2, #2[[1]] - blockSize/2, #1} &,
ImageData[
ImageCorrelate[part, template,
SquaredEuclideanDistance]], {2}], 1],
fitTerms, {x, y}]], {partsM, partsS}, 2];
这为我们提供了优化的第一个能量项：
variablesX = Array[dx, Dimensions[partsS]];
variablesY = Array[dy, Dimensions[partsS]];
matchEnergyFit =
Total[MapThread[#1 /. {x -> #2, y -> #3} &, {distancesFit,
variablesX, variablesY}, 2], 3];
variablesX/Y包含每个块的偏移量，并在matchEnergyFit应用了偏移量的情况下近似原始图像和镜像图像之间的平方欧几里得差。
单独优化这种能量会产生糟糕的结果（如果它完全收敛的话）。我们还希望偏移是平滑的，其中块相似性对偏移没有任何影响（例如，沿着直线或在白色背景中）。
所以我们为平滑度设置了第二个能量项：
smoothnessEnergy = Total[Flatten[
{
Table[
variablesX[[i, j - 1]] - 2 variablesX[[i, j]] +
variablesX[[i, j + 1]], {i, 1, Length[partsS]}, {j, 2,
Length[partsS[[1]]] - 1}],
Table[
variablesX[[i - 1, j]] - 2 variablesX[[i, j]] +
variablesX[[i + 1, j]], {i, 2, Length[partsS] - 1}, {j, 1,
Length[partsS[[1]]]}],
Table[
variablesY[[i, j - 1]] - 2 variablesY[[i, j]] +
variablesY[[i, j + 1]], {i, 1, Length[partsS]}, {j, 2,
Length[partsS[[1]]] - 1}],
Table[
variablesY[[i - 1, j]] - 2 variablesY[[i, j]] +
variablesY[[i + 1, j]], {i, 2, Length[partsS] - 1}, {j, 1,
Length[partsS[[1]]]}]
}^2]];
幸运的是，Mathematica 内置了约束优化：
allVariables = Flatten[{variablesX, variablesY}];
constraints = -blockSize/3. < # < blockSize/3. & /@ allVariables;
initialValues = {#, 0} & /@ allVariables;
solution =
FindMinimum[{matchEnergyFit + 0.1 smoothnessEnergy, constraints},
initialValues];
让我们看看结果：
grid = Table[{(j - 0.5)*blockSize - dx[i, j], (i - 0.5)*blockSize -
dy[i, j]}, {i, Length[partsS]}, {j, Length[partsS[[1]]]}] /.
solution[[2]];
Show[src, Graphics[
{Red,
Line /@ grid,
Line /@ Transpose[grid]
}]]
之前的0.1因素smoothnessEnergy是平滑能量相对于图像匹配能量项的相对权重。这些是不同权重的结果：
可能的改进：
- 就像我说的，首先执行严格的注册。在白色背景下，简单的基于图像时刻的配准应该可以正常工作。
- 这只是一步。您可以使用在一步中找到的偏移量并在第二步中改进它们，可能使用更小的搜索窗口或更小的块大小
- 我读过文章，他们完全没有块地这样做，但优化了每像素的偏移量。
- 尝试不同的平滑函数
