Znajdowanie symetrycznych regionów / wzorów na obrazie

Mam zestaw zdjęć przedstawiających średnią krzywiznę ludzkiej tylnej powierzchni.

Chcę „zeskanować” obraz w poszukiwaniu punktów, które mają podobne, odbite „odpowiedniki” w innej części obrazu (najprawdopodobniej symetryczne do linii środkowej, ale niekoniecznie, ponieważ mogą występować deformacje). Niektóre techniki łączenia obrazów wykorzystują to do „automatycznego wykrywania” podobnych punktów między obrazami, ale chcę je wykryć dla obu stron tego samego obrazu.

Ostatecznym celem jest znalezienie ciągłej, najprawdopodobniej zakrzywionej linii podłużnej, która adaptacyjnie dzieli plecy na symetryczne „połówki”.

Przykładowy obraz znajduje się poniżej. Zauważ, że nie wszystkie regiony są symetryczne (konkretnie tuż nad środkiem obrazu czerwony pionowy „pasek” odchyla się w prawo). Region ten powinien otrzymać złą ocenę, czy cokolwiek innego, ale wtedy lokalna symetria byłaby zdefiniowana z symetrycznych punktów umieszczonych dalej. W każdym razie będę musiał zaadaptować dowolny algorytm do mojej domeny aplikacji, ale moim celem jest strategia korelacji / splotu / dopasowania wzorców, myślę, że musi już być coś wokół.

(EDYCJA: poniżej jest więcej zdjęć i więcej wyjaśnień)

EDYCJA: zgodnie z prośbą, dołączę bardziej typowe obrazy, dobrze wychowane i problematyczne. Ale zamiast obrazów z kolorem, są to obrazy w skali szarości, więc kolor odnosi się bezpośrednio do wielkości danych, co nie zdarzyło się w przypadku kolorowego obrazu (tylko w celu komunikacji). Mimo że szare obrazy wydają się pozbawione kontrastu w porównaniu z kolorowymi, gradienty danych istnieją i można je w razie potrzeby uzupełnić kontrastem adaptacyjnym.

1) Obraz bardzo symetrycznego obiektu:

2) Zdjęcie tego samego obiektu w innym momencie. Chociaż jest więcej „funkcji” (więcej gradientów), nie „czuje się” tak symetrycznie jak wcześniej:

3) Cienki, młody podmiot z wypukłościami (wypukłe kości, oznaczone jaśniejszymi regionami) na linii środkowej zamiast częściej występującej wklęsłej linii środkowej:

4) Młody człowiek z odchyleniem kręgosłupa potwierdzonym przez zdjęcie rentgenowskie (zauważ asymetrie):

5) Typowy „przechylony” obiekt (choć w większości symetryczny wokół zakrzywionej linii środkowej i jako taki nie jest odpowiednio „zdeformowany”):

Każda pomoc jest mile widziana!

Jak powiedziałem w komentarzach, rejestracja obrazów medycznych jest tematem z wieloma dostępnymi badaniami i nie jestem ekspertem. Z tego, co przeczytałem, podstawową powszechnie stosowaną ideą jest zdefiniowanie odwzorowania między dwoma obrazami (w twoim przypadku obrazem i jego odbiciem lustrzanym), następnie zdefiniowanie terminów energii dla gładkości i podobieństwa obrazu, jeśli zastosowane zostanie odwzorowanie, i na koniec zoptymalizuj to mapowanie przy użyciu standardowych (lub czasem specyficznych dla aplikacji) technik optymalizacji.

Zhakowałem szybki algorytm w Mathematica, aby to zademonstrować. To nie jest algorytm, którego powinieneś używać w aplikacji medycznej, a jedynie demonstracja podstawowych pomysłów.

Najpierw ładuję twój obraz, odbij go i dzielę te obrazy na małe bloki:

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]

Normalnie robilibyśmy przybliżoną sztywną rejestrację (używając np. Punktów kluczowych lub momentów obrazu), ale twój obraz jest prawie wyśrodkowany, więc pominę to.

Jeśli spojrzymy na jeden blok i jego odpowiednik odbicia lustrzanego:

{partsS[[6, 10]], partsM[[6, 10]]}

Widzimy, że są podobne, ale przesunięte. Ilość i kierunek zmiany jest tym, co próbujemy ustalić.

Aby zmierzyć podobieństwo dopasowania, mogę użyć kwadratowej odległości euklidesowej:

ListPlot3D[
  ImageData[
   ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]], 
    SquaredEuclideanDistance]]]

niestety, użycie tych danych jest tym, że optymalizacja bezpośrednio była trudniejsza niż myślałem, więc zamiast tego zastosowałem przybliżenie drugiego rzędu:

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}]

Funkcja nie jest taka sama jak faktyczna funkcja korelacji, ale jest wystarczająco blisko, aby wykonać pierwszy krok. Obliczmy to dla każdej pary bloków:

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];

To daje nam nasz pierwszy termin energetyczny do optymalizacji:

variablesX = Array[dx, Dimensions[partsS]];
variablesY = Array[dy, Dimensions[partsS]];

matchEnergyFit = 
  Total[MapThread[#1 /. {x -> #2, y -> #3} &, {distancesFit, 
     variablesX, variablesY}, 2], 3];

variablesX/Yzawiera przesunięcia dla każdego bloku i matchEnergyFitprzybliża kwadratową różnicę euklidesową między obrazem oryginalnym a obrazem lustrzanym z zastosowanymi przesunięciami.

Sama optymalizacja tej energii dałaby słabe wyniki (gdyby w ogóle się zbiegła). Chcemy również, aby przesunięcia były gładkie, a podobieństwo bloku nie mówi nic o przesunięciu (np. Wzdłuż linii prostej lub na białym tle).

Dlatego ustanowiliśmy drugi termin energetyczny dla gładkości:

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]];

Na szczęście w Mathematica jest wbudowana ograniczona optymalizacja:

allVariables = Flatten[{variablesX, variablesY}];
constraints = -blockSize/3. < # < blockSize/3. & /@ allVariables;
initialValues = {#, 0} & /@ allVariables;
solution = 
  FindMinimum[{matchEnergyFit + 0.1 smoothnessEnergy, constraints}, 
   initialValues];

Spójrzmy na wynik:

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.1Czynnikiem zanim smoothnessEnergyto masa względna energia gładkość dostaje w stosunku do terminu obraz mecz energii. Oto wyniki dla różnych wag:

Możliwe ulepszenia:

• Jak powiedziałem, najpierw wykonaj sztywną rejestrację. Na białym tle prosta rejestracja oparta na momentach obrazu powinna działać poprawnie.
• To tylko jeden krok. Możesz użyć przesunięć znalezionych w jednym kroku i poprawić je w drugim kroku, być może z mniejszym oknem wyszukiwania lub mniejszymi rozmiarami bloku
• Czytałem artykuły, w których robią to w ogóle bez bloków, ale optymalizuję przesunięcie na piksel.
• Wypróbuj różne funkcje wygładzania

Interesujące pytanie. Po pierwsze, być może poszukujesz podejść opartych na detektorze kluczowych punktów i dopasowaniu. Obejmuje to SIFT (Scale-Invariant Feature Transform), SURF, ORB itp. Lub nawet prostsze podejście oparte wyłącznie na operatorze Harrisa (csce.uark.edu/~jgauch/library/Features/Harris.1988.pdf ). Z twojego postu nie jest jasne, czego próbowałeś, więc przepraszam, jeśli jestem tu naiwny.

Mówiąc to, pozwólcie, że zastosuję prostsze podejście z Morfologią Matematyczną (MM) dla zabawy :) Zdjęcia do wizualizacji wszystkich kroków są na końcu.

Wziąłem twój przykładowy obraz i przekonwertowałem go do przestrzeni kolorów L a b * za pomocą ImageMagick i użyłem tylko pasma L *:

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

0.png odpowiada pasmowi L *. Teraz jestem pewien, że masz rzeczywiste dane obrazu, ale mam do czynienia z artefaktami kompresji jpg, a co nie. Aby częściowo poradzić sobie z tym problemem, wykonałem otwarcie morfologiczne, a następnie zamknięcie morfologiczne płaską tarczą o promieniu 5. Jest to podstawowy sposób na zmniejszenie szumu przy MM, a biorąc pod uwagę promień dysku, niewielka część obrazu jest zmieniana. Następnie mój pomysł został oparty na tym pojedynczym obrazie, który ma duże szanse na niepowodzenie w innych przypadkach. Twój obszar zainteresowania jest wizualnie wyróżniony przez to, że jest ciemniejszy („gorętszy” na kolorowym obrazie), więc przypuszczałem, że statystyczny binarizator może działać dobrze. Użyłem podejścia Otsu, które jest automatyczne.

W tym momencie można wyraźnie wizualizować centralny obszar zainteresowania. Problem polega na tym, że w moim podejściu chciałem, aby był to zamknięty komponent, ale tak nie jest. Zaczynam od odrzucenia każdego podłączonego komponentu, który jest mniejszy niż największy (nie licząc tła jako jednego z nich). Ma to większą szansę na działanie w innych przypadkach, jeśli wynik binaryzacji był dobry. W twoim przykładowym obrazie jeden element jest podłączony do tła, więc nie jest odrzucany, ale nie powoduje problemów.

Jeśli nadal mnie śledzisz, nie udało nam się znaleźć faktycznego centralnego regionu zainteresowania. Oto moje zdanie na ten temat. Bez względu na to, jak osoba jest zakrzywiona (właściwie widzę pewne problematyczne przypadki), region przypomina pionową linię. W tym celu upraszczam bieżący obraz, wykonując otwarcie morfologiczne z pionową linią o długości 100. Ta długość jest czysto dowolna, jeśli nie masz problemów ze skalowaniem, nie jest to trudne do ustalenia. Teraz ponownie odrzucamy komponenty, ale na tym etapie byłem nieco bardziej ostrożny. Użyłem otwierania według obszaru z uzupełnieniem obrazu, aby odrzucić to, co uważałem za małe regiony, można to zrobić w bardziej kontrolowany sposób, wykonując coś w formie analizy granulometrii (również z MM).

Z grubsza mamy teraz trzy części: lewą część obrazu, środkową część i prawą część obrazu. Oczekuje się, że środkowa część będzie mniejszym składnikiem trójki, więc jest uzyskiwana w sposób trywialny.

Oto końcowy wynik, prawy dolny obraz jest tylko nałożonym obrazem po lewej stronie z oryginalnym. Poszczególne liczby nie są wyrównane, przepraszam za pośpiech.