<jupyter_start><jupyter_text>This exercise is based on the MATLAB implementation of LC2. The goal here is to become familiar with registration ofCT-US images. First we try LC2 on one sample image and compare it to SSD.Then we evaluate LC2 and SSD on registration of a pair of CT-US images.<jupyter_code>import imageio
import matplotlib.pyplot as plt
import numpy as np
from skimage import color
from skimage.transform import resize
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load the phantom image and resize<jupyter_code>img = imageio.imread("../data/SheppLogan_Phantom.png")
img = color.rgb2gray(img)
img = resize(img, (200, 200))<jupyter_output><empty_output><jupyter_text>Add noise to the phantom image to make the moving image<jupyter_code>img1 = img + np.random.normal(0, 0.15, img.shape)
plt.imshow(img1)<jupyter_output><empty_output><jupyter_text>Stack the gradient of the second image<jupyter_code>from skimage import filters

edges = filters.sobel(img)
img2 = np.stack((img, edges), axis=2)<jupyter_output><empty_output><jupyter_text>First, We define the global LC2 metric. The goal is to find the coefficients that best simulates the first image(later the US image) given the second image (later the CT image). The simulation is assumed a linear function ofintensity and gradient of the second image.<jupyter_code>def LC2Similarity(img1, img2):
    similarity, weight, measure = -1, 0, 0

    num = img1.shape[0] * img2.shape[1]
    ids = img1 > 0
    v1 =  np.var(img1[ids])

    if v1 > 1e-12:
        if np.sum(ids) > num / 2:
            img2r = np.concatenate((np.reshape(img2, (num, img2.shape[2])), np.ones((num, 1))), axis=1)
            pimg2r = np.linalg.pinv(img2r[ids.flatten(), :])
            parameter = np.matmul(pimg2r, img1[ids])
            similarity = 1 - (np.var(img1[ids] - (np.matmul(img2r[ids.flatten(), :], parameter))) / v1)
            weight = np.sqrt(v1)
            measure = weight * similarity
    if similarity == -1:
        similarity = 0

    return similarity, measure, weight<jupyter_output><empty_output><jupyter_text>Test against rotation<jupyter_code># constants:
ang_range = 25
ang_bins = 31
tra_range = 11
tra_step_size = 1

from scipy.ndimage import rotate

angles = np.linspace(-ang_range, ang_range, ang_bins)
sims_lc2 = []
sims_ssd = []
for angle in angles:
    img1_r = rotate(img1, angle, reshape=False)
    similarity, _, _ = LC2Similarity(img1_r, img2)
    sims_lc2.append(similarity)
    sims_ssd.append(np.sum((img1_r - img) ** 2))

sims_ssd = 1 - sims_ssd / np.max(sims_ssd)
plt.plot(angles, sims_lc2)
plt.plot(angles, sims_ssd)
plt.xlabel("angle")
plt.ylabel("similarity")<jupyter_output><empty_output><jupyter_text>Test against translation<jupyter_code>from scipy.ndimage import affine_transform

affine = np.eye(3)
sims = []

grid_x, grid_y = np.mgrid[-tra_range:tra_range:tra_step_size, -tra_range:tra_range:tra_step_size]

for x, y in zip(grid_x.flatten(), grid_y.flatten()):
    affine[0, 2] = x
    affine[1, 2] = y
    img1_t = affine_transform(img1, affine)
    similarity, _, _ = LC2Similarity(img1_t, img2)
    sims.append(similarity)

sims = np.reshape(sims, grid_x.shape)
h = plt.contourf(grid_x, grid_y, sims)<jupyter_output><empty_output><jupyter_text>As can be seen from the graphs, our global LC2 reaches its maximum value on the correct alignment. For the phantom image,SSD also reaches its extremum on the correct alignment.Now, we move to a pair of US-CT images. First we load the two image and pre-calculate the gradient of the second image.<jupyter_code>img1 = imageio.imread("../data/Ultrasound_slice.png")
img1 = color.rgb2gray(img1)
img1 = resize(img1, (200, 200))
plt.imshow(img1)
img = imageio.imread("../data/CT_slice.png")
img = color.rgb2gray(img)
img = resize(img, (200, 200))
plt.imshow(img)

from skimage import filters

edges = filters.sobel(img)
img2 = np.stack((img, edges), axis=2)<jupyter_output><ipython-input-9-a8da66ead2ab>:2: FutureWarning: The behavior of rgb2gray will change in scikit-image 0.19. Currently, rgb2gray allows 2D grayscale image to be passed as inputs and leaves them unmodified as outputs. Starting from version 0.19, 2D arrays will be treated as 1D images with 3 channels.
  img = color.rgb2gray(img)<jupyter_text>Now test against rotation and translation as before and compare it to SSD.<jupyter_code>from scipy.ndimage import rotate

angles = np.linspace(-ang_range, ang_range, ang_bins)
sims_lc2 = []
sims_ssd = []
for angle in angles:
    img1_r = rotate(img1, angle, reshape=False)
    similarity, _, _ = LC2Similarity(img1_r, img2)
    sims_lc2.append(similarity)
    sims_ssd.append(np.sum((img1_r - img) ** 2))

sims_ssd = 1 - sims_ssd / np.max(sims_ssd)
plt.plot(angles, sims_lc2)
plt.plot(angles, sims_ssd)
plt.xlabel("angle")
plt.ylabel("similarity")
from scipy.ndimage import affine_transform

affine = np.eye(3)
sims = []

grid_x, grid_y = np.mgrid[-tra_range:tra_range:tra_step_size, -tra_range:tra_range:tra_step_size]

for x, y in zip(grid_x.flatten(), grid_y.flatten()):
    affine[0, 2] = x
    affine[1, 2] = y
    img1_t = affine_transform(img1, affine)
    similarity, _, _ = LC2Similarity(img1_t, img2)
    sims.append(similarity)

sims = np.reshape(sims, grid_x.shape)
h = plt.contourf(grid_x, grid_y, sims)<jupyter_output><empty_output><jupyter_text>While SSD is failing, our global LC2 is showing aFor the next part your task is to write the code for patch-wise LC2 using the one we already created.<jupyter_code>def LC2SimilarityPatch(img1, img2, patchsize):
    max_x = img1.shape[0]
    max_y = img1.shape[1]
    min_size = ((2 * patchsize + 1)^2) / 2
    measure = np.zeros(img1.shape)
    weights = np.zeros(img1.shape)

    for y in np.arange(0, max_y, step = patchsize):
        for x in np.arange(0, max_x, step = patchsize):
            # TODO: calculate the metrics for each patch
            x_min = np.max((0, x - patchsize))
            x_max = np.min((max_x - 1, x + patchsize))
            y_min = np.max((1, y - patchsize))
            y_max = np.min((max_y - 1, y + patchsize))
            patch1 = img1[x_min:x_max, y_min:y_max]
            patch2 = img2[x_min:x_max, y_min:y_max, :]
            
            if np.size(patch1) > min_size:
                _, measure[x, y], weights[x, y] = LC2Similarity(patch1, patch2)

    similarity = np.nansum(measure) / np.nansum(weights) # TODO: average over patches
    return similarity<jupyter_output><empty_output><jupyter_text>Test against rotation and translation. Check different values for patchsize.<jupyter_code>LC2SimilarityPatch(img1, img2, patchsize=5)
angles = np.linspace(-ang_range, ang_range, ang_bins)
sims_lc2 = []
sims_ssd = []
for angle in angles:
    img1_r = rotate(img1, angle, reshape=False)
    similarity = LC2SimilarityPatch(img1_r, img2, patchsize=5)
    sims_lc2.append(similarity)
    sims_ssd.append(np.sum((img1_r - img) ** 2))

sims_ssd = 1 - sims_ssd / np.max(sims_ssd)
plt.plot(angles, sims_lc2)
plt.plot(angles, sims_ssd)
plt.xlabel("angle")
plt.ylabel("similarity")
affine = np.eye(3)
sims = []

grid_x, grid_y = np.mgrid[-tra_range:tra_range:tra_step_size, -tra_range:tra_range:tra_step_size]

for x, y in zip(grid_x.flatten(), grid_y.flatten()):
    affine[0, 2] = x
    affine[1, 2] = y
    img1_t = affine_transform(img1, affine)
    similarity = LC2SimilarityPatch(img1_t, img2, patchsize=5)
    sims.append(similarity)

sims = np.reshape(sims, grid_x.shape)
h = plt.contourf(grid_x, grid_y, sims)<jupyter_output>D:\Anaconda\lib\site-packages\numpy\core\fromnumeric.py:3702: RuntimeWarning: Degrees of freedom <= 0 for slice
  return _methods._var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
D:\Anaconda\lib\site-packages\numpy\core\_methods.py:221: RuntimeWarning: invalid value encountered in true_divide
  arrmean = um.true_divide(arrmean, div, out=arrmean, casting='unsafe',
D:\Anaconda\lib\site-packages\numpy\core\_methods.py:253: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
{'fraction_non_alphanumeric': 0.11863136863136862, 'fraction_numerical': 0.0265984015984016, 'mean_word_length': 4.342895263509006, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '304214', 'license': 'MIT', 'n_tokens_mistral': 3147, 'n_tokens_neox': 2975, 'n_words': 715, 'path': 'IN2021 CAMP Medical Image Registration Exercise/intensity-based/US-CT registration/US-CT-registration.ipynb', 'repo_name': 'CornellLenard/TUM-Computer-Aided-Medical-Procedures', 'size': '8008'}