Matlab - Alineamiento y comparación de coordenadas

Buenas, partiendo de dos imagenes DICOM adjuntadas, una de resonancia y otra de TAC del mismo maniquí, como vemos en las imágenes abajo, estoy haciendo un código en el cual se roten las imágenes para que los agujeros queden alineados, sacar las coordenadas del centroide de dichos agujeros y comparar las posiciones de cada uno de estos centroides y ver la diferencia de coordenadas de la resonancia y el TAC.
Mi primer problema es que sólo hacer la rotación de las imágenes con la función "Imrotate" probando ángulos y ver cuando queda recta, ¿alguna forma de hacerlo automáticamente que se alinee el circulo circulo central con el de arriba o el del al lado?
Otro problema que tengo es que me da las coordenadas de los centroides de los agujeros desordenadas, me gustaría ordenarlas de alguna forma, por ejemplo, de izquierda a derecha y de arriba abajo y poder comparar las coordenadas de la imagen de la resonancia con las del TAC y poder ver la diferencia entre dichas coordenadas.
Les dejo el código que tengo para que me digan sus comentarios, gracias de antemano.



clear all
close all
clc
I=dicomread('MR_T2_AX.dcm');
I=imrotate(I,8);
figure(1);imshow(I,[]);
umb=graythresh(I);
bw=im2bw(I,umb);
bw=not(bw);
s=strel('disk',1);
bw=imdilate(bw,s);
figure(2);imshow(bw,[]);
[L Ne]=bwlabel(bw);
prop=regionprops(L);
r=find([prop.Area]<100);
g=find([prop.Area]>100);

%encontrar los círculos con área pequeña
for n=1:size(r,2)
rectangle('Position',prop(r(n)).BoundingBox,'EdgeColor','r','Linewidth',1);
end

%encontrar círculos con área grande
for n=1:size(g,2)
rectangle('Position',prop(g(n)).BoundingBox,'EdgeColor','g','Linewidth',1);
end

% paso las coordenadas de pixel a mm
j=dicominfo('MR_T2_AX.dcm');
dxy=j.PixelSpacing;
for n=2:Ne;
H(n-1,:)=[prop(n).Centroid]*[dxy(1):dxy(2)];
end

% graficar las coordenadas de los agujeros
a=H(:,1)
b=H(:,2)
figure(3);plot(a,b,'o');

% hacer lo mismo con la imagen del TAC
I1=dicomread('CT.dcm');
I1=imrotate(I1,1);
figure(4);imshow(I1,[]);
s1=strel('disk',2);
j1=imdilate(I1,s);
bw1=im2bw(j1,0.02);
figure(5);imshow(bw1,[]);
[L1 Ne1]=bwlabel(bw1);
prop1=regionprops(L1);
r1=find([prop1.Area]<100);
g1=find([prop1.Area]>100);

%encontrar los círculos con área pequeña
for n=1:size(r1,2)
rectangle('Position',prop(r1(n)).BoundingBox,'EdgeColor','r','Linewidth',1);
end

%encontrar círculos con área grande
for n=1:size(g1,2)
rectangle('Position',prop(g1(n)).BoundingBox,'EdgeColor','g','Linewidth',1);
end

% paso las coordenadas de pixel a mm
j1=dicominfo('CT.dcm');
dxy1=j1.PixelSpacing;

for n=1:Ne1;
H1(n,:)=[prop1(n).Centroid]*[dxy1(1):dxy1(2)];
end
% graficar las coordenadas
for n=1:97
a1(n,:)=H1(n);
b1(n,:)=H1(n+97);
figure(6);plot(a1,b1,'o');
end