As I've said in the comments, medical image registration is a topic with lots of research available, and I'm not an expert. From what I've read, the basic idea commonly used is to define a mapping between two images (in your case an image and its mirror image), then define energy terms for smoothness and for image similarity if the mapping is applied, and finally optimize this mapping using standard (or sometimes application-specific) optimization techniques.
I've hacked together a quick algorithm in Mathematica to demonstrate this. This is not an algorithm you should use in a medical application, only a demonstration of the basic ideas.
First, I load your image, mirror it and split these images into little blocks:
src = ColorConvert[Import["https://i.sstatic.net/jf709.jpg"],
"Grayscale"];
mirror = ImageReflect[src, Left -> Right];
blockSize = 30;
partsS = ImagePartition[src, {blockSize, blockSize}];
partsM = ImagePartition[mirror, {blockSize, blockSize}];
GraphicsGrid[partsS]
Normally, we would do an approximate rigid registration (using e.g. keypoints or image moments), but your image is almost centered, so I'll skip this.
If we look at one block and it's mirror-image counterpart:
{partsS[[6, 10]], partsM[[6, 10]]}
We can see that they're similar, but shifted. The amount and direction of shift is what we're trying to find out.
To quantify the match similarity, I can use the squared euclidean distance:
ListPlot3D[
ImageData[
ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]],
SquaredEuclideanDistance]]]
sadly, using this data is the optimization directly was harder than I thought, so I used a 2nd order approximation instead:
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}]
The function is not the same as the actual correlation function, but it's close enough for a first step. Let's calculate this for every pair of blocks:
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];
This gives us our first energy term for the optimization:
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 contains the offsets for each block, and matchEnergyFit approximates the squared euclidean difference between the original image and the mirrored image with the offsets applied.
Optimizing this energy alone would give poor results (if it converged at all). We also want the offsets to be smooth, where the block similarity tells as nothing about the offset (e.g. along a straight line or in the white background).
So we set up a second energy term for smoothness:
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]];
Fortunately, constrained optimization is built-in in Mathematica:
allVariables = Flatten[{variablesX, variablesY}];
constraints = -blockSize/3. < # < blockSize/3. & /@ allVariables;
initialValues = {#, 0} & /@ allVariables;
solution =
FindMinimum[{matchEnergyFit + 0.1 smoothnessEnergy, constraints},
initialValues];
Let's look at the result:
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]
}]]
The 0.1 factor before smoothnessEnergy is the relative weight the smoothness energy gets in relation to the image match energy term. These are results for different weights:
Possible improvements:
- Like I said, perform a rigid registration first. With a white background, simple image moments-based registration should work fine.
- This is only one step. You can use the offsets you found in one step and improve them in a second step, maybe with a smaller search window or smaller block sizes
- I've read articles where they do this without blocks at all, but optimize an offset per pixel.
- Try different smoothness functions
