# Finding symmetric regions/patterns in image

I have a set of images representing the mean curvature of a human back surface.

What I want to do is to "scan" the image for points which have similar, reflected "counterparts" in some other part of the image (most likely symmetrical to the midline, but not necessarily since there might be deformities). Some image-stitching techniques use this to "auto-detect" similar points between images, but I want to detect them for both sides of the same image.

The ultimate goal is to find a continuous, most probably curved, longitudinal line that adaptively divides the back in symmetric "halves".

A sample image is placed below. Notice that not all the regions are symmetrical (specifically, just above the center of the image, the red vertical "strip" deviates to the right). That region should receive a bad score, or whatever, but then the local symmetry would be defined from symmetrical points placed farther away. In any case, I'll have to adapt any algorythm to my application domain, but what I'm after is som correlation/convolution/pattern-matching strategy, there must already be something around, I think.

(EDIT: there are more images below, and some more explanation)

EDIT: as requested, I'll include more typical images, either well-behaved and problematic. But instead of colormapped images, they are grayscale ones, so that color relates directly to data magnitude, which didn't happened with the colored image (provided just for communication). Although the gray images seem to lack contrast compared to the colored one, the data gradients are there and can be brought up with some adaptive contrast if desired.

1) Image of a very symmetrical subject:

2) Image of the same subject at a different moment. Although there are more "features" (more gradients), it doesn't "feel" so symmetric as before:

3) A thin young subject, with convexities (bony protrusions, denoted by lighter regions) at the midline instead of the more common concave midline:

4) A young person with spinal deviation confirmed by X-Ray (note the asymmetries):

5) The typical "tilted" subject (although mostly symmetrical around the curved midline, and as such not properly "deformed"):

Any help is much welcome!

-
Why not just use the spine as the divider? –  Jim Clay Nov 28 '12 at 20:31
@JimClay: I suspect the spine is the part being measured, relative to the actual axis of symmetry of the rest of the image –  endolith Nov 28 '12 at 21:46
"Some image-stitching techniques use this to "auto-detect" similar points between images" Make a flipped copy of the image and then use one of those. :) –  endolith Nov 28 '12 at 21:48
Couldn't you simply mirror the image along the Y axis and use a registration algorithm? Because there's already lots of research on flexible/nonparametric registration algorithms you could build upon. –  nikie Nov 28 '12 at 22:16
JimClay, the spine is what I want to find, I don't know where it is; Endolith, my question involves people telling me the names of some of those algorithms, I didn't find any yet. And Nikie, that's the whole point, but I don't KNOW any of those algorithms, that's why I'm asking in the question in the first place :o) –  heltonbiker Nov 28 '12 at 22:55

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["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]


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
-
answer too long to read just for fun, but the final image is pretty indicative: it looks amazing :D –  penelope Nov 29 '12 at 10:38
This answer was very enlightening. I'll need some time to swallow it, but most probably the non-rigid registration technique is what I'll need to use. Fortunately you provided some conceptual details, so in the worst case I can figure out some similar approach. In the meantime, I'll update the question with more images. Thanks for now! –  heltonbiker Nov 29 '12 at 12:34

Interesting question. First, maybe you are after approaches based on interest keypoint detector and matching. This would include SIFT (Scale-Invariant Feature Transform), SURF, ORB, etc ... or even a simpler approach based solely on the Harris operator (csce.uark.edu/~jgauch/library/Features/Harris.1988.pdf). It is not clear from your post what you have tried, so I'm sorry if I'm being naïve here.

Said that, let me take a simpler approach with Mathematical Morphology (MM) just for fun :) Images for visualization of all of the steps are at the end.

I took your sample image and converted it to the L*a*b* colorspace using ImageMagick and used only the L* band:

convert x.jpg -colorspace Lab -separate %d.png


0.png corresponds to the L* band. Now, I'm sure you have the actual image data, but I'm dealing with jpg compression artifacts and what not. To partially handle this issue, I performed a morphological opening followed by a morphological closing with a flat disk of radius 5. This is a basic way to reduce noise with MM, and given the disk radius not much of the image is changed. Next my idea was based on this single image, which has a great chance of failing for other cases. Your region of interest is visually distinguished by being darker ("hotter" in your color image), so I supposed a statistically-based binarizer could perform well. I used the Otsu's approach, which is an automatic one.

At this point, it is possible to clearly visualize the central region of interest. The problem is that, in my approach, I wanted it to be a closed component but it is not. I start by discarding every connected component that is smaller than the largest one (not counting background as one of them). This has greater chance to work in other cases if the binarization result was a good one. In your example image, there is one component connected to the background, so it is not discarded but it does not cause problems.

If you are still following me, we have yet to find the actual supposed central region of interest. Here is my take on it. No matter how curved the person is (actually I can see certain problematic cases), the region resembles a vertical line. To that end, I simplify the current image by performing a morphological opening with a vertical line of length 100. This length is purely arbitrary, if you do not have scaling issues then this is not a hard value to determine. Now we again discard components, but I was a bit more careful at this step. I used opening by area with the image's complement to discard what I considered small regions, this could be done in a more controlled way by performing something in the form of granulometrics analysis (from MM too).

We roughly have three pieces now: the left part of the image, the central part, and the right part of the image. The central part is expected to be the smaller component of the three, so it is trivially obtained.

Here is the final result, the bottom right image is just the superposed image to its left with the original one. The individual figures are not all aligned, sorry for the hurry.

-
Thank you very very much for your kind interest, but your approach should consider certain properties of my data (not a complaint, just a detailing): 1) The actual data is a 2D array of floats, colormapped with a divergent red-yellow-green colormap in Python's matplotlib. I don't think working with the color data would be conceptually correct, the images are shown for communication purposes only; 2) The actual data relates to surface curvature (convex vs. concave), with red parts being concave and green ones being convex. Symmetrical axis doesn't necessarily falls in a concave region. –  heltonbiker Nov 29 '12 at 12:24
I'll very soon add some more images (and replace this one) with grayscale ones, so that the images themselves might be used for testing, eliminating the danger of dynamic range distortion due to color. –  heltonbiker Nov 29 '12 at 12:26
The data is not available yet, unfortunately. The grayscale images are at best an approximation of it. –  mmgp Nov 29 '12 at 13:19
I believe the approximation most probably is enough, but I don't mind providing the actual data. I can post some public DropBox download links, just don't know in which file format. –  heltonbiker Nov 29 '12 at 13:27