Совпадение угловых точек в КТ и МР позвоночника
Я пытаюсь сопоставить угловые точки Харриса, обнаруженные при КТ и МРТ. Однако, когда я пытаюсь сопоставить эти точки, используя нормализованную взаимную корреляцию и NMI, но это показывает неожиданный результат. Я не могу определить резонанс. Нужно ли предварительно обрабатывать входные изображения, используя какую-то технику, поскольку изображения являются мультимодальными?
import numpy as np
import cv2
from matplotlib import pyplot as plt
from skimage import exposure
import scipy.ndimage as ndimage
import time
from skimage.feature import corner_harris, corner_subpix, corner_peaks
from sklearn.metrics.cluster import normalized_mutual_info_score
start_time = time.time()
EPS = np.finfo(float).eps
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
plt.close('all')
def extract_keypoints(keypoints):
i = 0
temp_array=np.zeros((500,2),dtype=np.int)
for point in keypoints:
temp=point.pt
temp_array[i][0]=temp[0]
temp_array[i][1]=temp[1]
i=i+1
return temp_array
def correlation_coefficient(patch1, patch2):
product = np.mean((patch1 - patch1.mean()) * (patch2 - patch2.mean()))
stds = patch1.std() * patch2.std()
if stds == 0:
return 0
else:
product /= stds
return product
def get_descriptors(image,filtered_coords,wid=5):
""" For each point return pixel values around the point
using a neighbourhood of width 2*wid+1. (Assume points are
extracted with min_distance > wid). """
desc = []
for coords in filtered_coords:
patch = image[coords[0]-wid:coords[0]+wid+1,
coords[1]-wid:coords[1]+wid+1].flatten()
desc.append(patch)
return desc
def match(desc1,desc2,threshold=0.1):
""" For each corner point descriptor in the first image,
select its match to second image using
normalized cross correlation. """
n = len(desc1[0])
# pair-wise distances
d = -np.ones((len(desc1),len(desc2)))
for i in range(len(desc1)):
for j in range(len(desc2)):
d1 = (desc1[i] - np.mean(desc1[i])) / np.std(desc1[i])
d2 = (desc2[j] - np.mean(desc2[j])) / np.std(desc2[j])
ncc_value = sum(d1 * d2) / (n-1)
if ncc_value > threshold:
d[i,j] = ncc_value
ndx = np.argsort(-d)
matchscores = ndx[:,0]
return matchscores
def match_twosided(desc1,desc2,threshold=0.1):
""" Two-sided symmetric version of match(). """
matches_12 = match(desc1,desc2,threshold)
matches_21 = match(desc2,desc1,threshold)
ndx_12 = np.where(matches_12 >= 0)[0]
# remove matches that are not symmetric
for n in ndx_12:
if matches_21[matches_12[n]] != n:
matches_12[n] = -1
return matches_12
def match_MI(desc1,desc2,threshold=0.001):
""" For each corner point descriptor in the first image,
select its match to second image using
normalized cross correlation. """
# pair-wise distances
d = -np.ones((len(desc1),len(desc2)))
for i in range(len(desc1)):
for j in range(len(desc2)):
MI_value = normalized_mutual_info_score(desc1[i],desc2[j])
if MI_value > threshold:
d[i,j] = MI_value
ndx = np.argsort(-d)
matchscores = ndx[:,0]
return matchscores
def match_twosided_MI(desc1,desc2,threshold=0.001):
""" Two-sided symmetric version of match(). """
matches_12 = match_MI(desc1,desc2,threshold)
matches_21 = match_MI(desc2,desc1,threshold)
ndx_12 = np.where(matches_12 >= 0)[0]
# remove matches that are not symmetric
for n in ndx_12:
if matches_21[matches_12[n]] != n:
matches_12[n] = -1
return matches_12
def appendimages(im1,im2):
""" Return a new image that appends the two images side-by-side. """
# select the image with the fewest rows and fill in enough empty rows
rows1 = im1.shape[0]
rows2 = im2.shape[0]
if rows1 < rows2:
im1 = np.concatenate((im1,np.zeros((rows2-rows1,im1.shape[1]))),axis=0)
elif rows1 > rows2:
im2 = np.concatenate((im2,np.zeros((rows1-rows2,im2.shape[1]))),axis=0)
# if none of these cases they are equal, no filling needed.
return np.concatenate((im1,im2), axis=1)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
""" Show a figure with lines joining the accepted matches
input: im1,im2 (images as arrays), locs1,locs2 (feature locations),
matchscores (as output from 'match()'),
show_below (if images should be shown below matches). """
im3 = appendimages(im1,im2)
if show_below:
im3 = np.vstack((im3,im3))
plt.show(im3)
cols1 = im1.shape[1]
for i,m in enumerate(matchscores):
if m>0:
plt.plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
plt.axis('off')
CT= cv2.imread('F:\MyTransforms_HP\DataSet1\CT.jpg')
MR= cv2.imread('F:\MyTransforms_HP\DataSet1\MR.jpg')
img1= cv2.imread('F:\MyTransforms_HP\DataSet1\CT.jpg',0)
img2= cv2.imread('F:\MyTransforms_HP\DataSet1\MR.jpg',0)
#img1=cv2.resize(img,(10,10),0,0)
CT_size=img1.shape
MR_size=img2.shape
#img11 = exposure.adjust_gamma(img1, 0.6)
p2, p98 = np.percentile(img2, (2, 60))
img21 = exposure.rescale_intensity(img2, in_range=(p2, p98))
coords1 = corner_peaks(corner_harris(img1), min_distance=10)
coords2 = corner_peaks(corner_harris(img21), min_distance=10)
NoCT_Corners=len(coords1)
NoMR_Corners=len(coords2)
print("Harris Corners detetcted--- %s seconds ---" % (time.time() - start_time))
print("No. of Corners in CT : %s" %NoCT_Corners)
print("No. of Corners in MR : %s" %NoMR_Corners)
des1=get_descriptors(img1,coords1)
des2=get_descriptors(img2,coords2)
matches_NCC=match_twosided(des1,des2)
des2=get_descriptors(img21,coords2)
matches_MI=match_twosided_MI(des1,des2)
#plot_matches(img1, img2, coords1, coords2, matches_12)
FinalMatch_NCC=np.zeros(shape=(36,2),dtype=np.int64)
k1=0
for i in range (0,NoCT_Corners):
if (matches_NCC[i]!=-1):
FinalMatch_NCC[k1][0]=i
FinalMatch_NCC[k1][1]=matches_NCC[i]
k1=k1+1
FinalMatch_MI=np.zeros(shape=(36,2),dtype=np.int64)
k2=0
for i in range (0,NoCT_Corners):
if (matches_MI[i]!=-1):
FinalMatch_MI[k2][0]=i
FinalMatch_MI[k2][1]=matches_MI[i]
k2=k2+1
im3 = appendimages(img1,img2)
im4 = appendimages(img1,img2)
#Draw lines
count = 0
for i in range (0,k1):
x1=coords1[FinalMatch_NCC[i][0]][1]
y1=coords1[FinalMatch_NCC[i][0]][0]
x2=coords2[FinalMatch_NCC[i][1]][1]+512
y2=coords2[FinalMatch_NCC[i][1]][0]
cv2.line(im3,(x1,y1),(x2,y2),(255,0,0),2)
cv2.circle(im3,(x1,y1),3,(255,255,255), -1)
cv2.circle(im3,(x2,y2),3,(255,255,255), -1)
#plt.close('all')
for i in range (0,k2):
x1=coords1[FinalMatch_MI[i][0]][1]
y1=coords1[FinalMatch_MI[i][0]][0]
x2=coords2[FinalMatch_MI[i][1]][1]+512
y2=coords2[FinalMatch_MI[i][1]][0]
cv2.line(im4,(x1,y1),(x2,y2),(255,0,0),1)
cv2.circle(im4,(x1,y1),3,(255,255,255), -1)
cv2.circle(im4,(x2,y2),3,(255,255,255), -1)
fig1, ax = plt.subplots(1,2)
ax[0].imshow(CT, interpolation='nearest', cmap=plt.cm.gray)
ax[0].plot(coords1[:, 1], coords1[:, 0], '.r', markersize=3)
ax[0].axis((0, 512,512, 0))
ax[1].imshow(MR, interpolation='nearest', cmap=plt.cm.gray)
ax[1].plot(coords2[:, 1], coords2[:, 0], '.r', markersize=3)
ax[1].axis((0, 512,512, 0))
plt.show()
plt.figure(2);
plt.imshow(im3);
plt.show();
plt.figure(3);
plt.imshow(im4);
plt.show();
[Harris Corners][1]