Import necessary libraries



In [1]:

    
%matplotlib inline



In [20]:

    
import matplotlib.pyplot as plt
import matplotlib
import ndreg
from ndreg import preprocessor, util, plotter
import SimpleITK as sitk



In [3]:

    
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)



In [4]:

    
def myshow(img, cmap='gray', colorbar=False):
    plt.imshow(sitk.GetArrayViewFromImage(img), cmap=cmap)
    if colorbar: plt.colorbar()
    plt.axis('off')
    plt.show()

Some metadata is required before registration



In [7]:

    
params = {
    # input image path
    'image_path': './Thy1eYFP_Control_9.tiff',
    # voxel spacing is in mm and corresponds to (x, y, z) spacing
    'image_spacing': (0.04128, 0.04128, 0.04128),
    'image_orientation': 'rpi',
    # the modality can be 'lavision' or 'colm'
    'image_modality': 'lavision',
    'atlas_spacing': (0.05, 0.05, 0.05),
    'atlas_path': './ARA_50um.tiff',
}

Load the sample data



In [8]:

    
img = util.imgRead(params['image_path'])
img.SetSpacing(params['image_spacing'])
atlas = util.imgRead(params['atlas_path'])
atlas.SetSpacing(params['atlas_spacing'])



In [9]:

    
plotter.imgShow(img, vmax=2000)



In [10]:

    
plotter.imgShow(atlas, vmax=400)

Preprocessing

This step preprocesses the input CLARITY images by resampling them to match the resolution of the atlas, bias correcting the images, and normalizing them by subtracting the mean and dividing by the standard deviation of the image intensities.



In [11]:

    
img_p = preprocessor.preprocess_brain(img, 
                                      params['atlas_spacing'], 
                                      params['image_modality'],
                                      params['image_orientation'])

Registration

We want to obtain the parameters to transform the original image to the new image. The transformation from the original image to the new image can be described as a composition of an affine transformation which can perform a combination of translation, scaling, rotation, and shear and deformable registration called LDDMM.

The output of this method is the atlas registered to the raw data



In [12]:

    
atlas_registered = ndreg.register_brain(atlas, img_p)









    



Step 0: alpha=0.05, beta=0.05, scale=0.0625
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
0.	85211.2, 89.4791, 0, 85121.7 (96.8504%), 2.750000e-03
1.	80645.1, 88.5653, 0, 80556.6 (91.6562%), 1.512500e-03
2.	74649.1, 124.356, 0, 74524.8 (84.7933%), 1.663750e-03
3.	72169.8, 157.397, 0, 72012.4 (81.9348%), 1.830125e-03
4.	70108.4, 214.212, 0, 69894.2 (79.5246%), 2.013138e-03
5.	68309.2, 280.522, 0, 68028.7 (77.4022%), 2.214451e-03
6.	67787.6, 327.251, 0, 67460.4 (76.7555%), 1.217948e-03
7.	66422.6, 376.188, 0, 66046.4 (75.1467%), 1.339743e-03
8.	66152.8, 440.439, 0, 65712.4 (74.7667%), 1.473717e-03
9.	65167.4, 511.049, 0, 64656.4 (73.5651%), 1.621089e-03
10.	64353.4, 553.485, 0, 63799.9 (72.5907%), 8.915990e-04
11.	62982.5, 600.129, 0, 62382.4 (70.9779%), 9.807589e-04
12.	62903.5, 654.72, 0, 62248.7 (70.8258%), 1.078835e-03
13.	62119.7, 719.976, 0, 61399.7 (69.8597%), 1.186718e-03
14.	61991.5, 788.65, 0, 61202.9 (69.6358%), 1.305390e-03
15.	60786.9, 829.161, 0, 59957.8 (68.2191%), 7.179645e-04
16.	60233.8, 873.171, 0, 59360.7 (67.5398%), 7.897610e-04
17.	59938.9, 923.912, 0, 59015 (67.1464%), 8.687371e-04
18.	59282, 981.598, 0, 58300.4 (66.3334%), 9.556108e-04
19.	59239.2, 997.566, 0, 58241.6 (66.2665%), 2.627930e-04
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
20.	59232.8, 998.66, 0, 58234.1 (66.258%), 1.806702e-05
21.	59229.3, 999.262, 0, 58230 (66.2533%), 9.936859e-06
22.	59228.3, 999.427, 0, 58228.9 (66.2521%), 2.732636e-06
23.	59228.2, 999.45, 0, 58228.8 (66.2519%), 3.757375e-07
24.	59228.1, 999.462, 0, 58228.7 (66.2518%), 2.066556e-07
25.	59228.1, 999.463, 0, 58228.7 (66.2518%), 1.420757e-08
E = 59278.1 (66.3086%)
Length = 27.0634
Time = 22.0763s (0.367939m)

Step 1: alpha=0.05, beta=0.05, scale=0.125
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
0.	77539.7, 16.346, 0, 77523.3 (95.9591%), 2.750000e-03
1.	76357, 22.3631, 0, 76334.6 (94.4876%), 3.025000e-03
2.	74659.4, 34.5098, 0, 74624.9 (92.3714%), 1.663750e-03
3.	74222.6, 45.9297, 0, 74176.7 (91.8166%), 1.830125e-03
4.	73400.4, 52.5644, 0, 73347.8 (90.7906%), 1.006569e-03
5.	73020.9, 61.1694, 0, 72959.7 (90.3102%), 1.107226e-03
6.	72753.1, 71.4347, 0, 72681.7 (89.9661%), 1.217948e-03
7.	72439.8, 83.4901, 0, 72356.3 (89.5633%), 1.339743e-03
8.	72114.6, 97.6144, 0, 72017 (89.1433%), 1.473717e-03
9.	71760, 114.15, 0, 71645.9 (88.6839%), 1.621089e-03
10.	71379.3, 133.575, 0, 71245.7 (88.1886%), 1.783198e-03
11.	71005.9, 156.187, 0, 70849.7 (87.6984%), 1.961518e-03
12.	70627.5, 182.414, 0, 70445.1 (87.1975%), 2.157670e-03
13.	70432.8, 212.341, 0, 70220.5 (86.9195%), 2.373436e-03
14.	70134.7, 229.795, 0, 69904.9 (86.5289%), 1.305390e-03
15.	69990.1, 247.866, 0, 69742.2 (86.3275%), 1.435929e-03
16.	69882, 270.555, 0, 69611.5 (86.1657%), 1.579522e-03
17.	69499.4, 281.099, 0, 69218.3 (85.6791%), 8.687371e-04
18.	69403.9, 294.599, 0, 69109.3 (85.5441%), 9.556108e-04
19.	69272.1, 309.469, 0, 68962.6 (85.3625%), 1.051172e-03
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
20.	69148.9, 326.254, 0, 68822.7 (85.1893%), 1.156289e-03
21.	68965.6, 345.109, 0, 68620.5 (84.9391%), 1.271918e-03
22.	68808.1, 366.449, 0, 68441.6 (84.7177%), 1.399110e-03
23.	68620.6, 390.511, 0, 68230.1 (84.4558%), 1.539021e-03
24.	68444.9, 417.819, 0, 68027.1 (84.2045%), 1.692923e-03
25.	68270.5, 448.477, 0, 67822 (83.9507%), 1.862215e-03
26.	68070.6, 483.877, 0, 67586.8 (83.6595%), 2.048437e-03
27.	67975.6, 502.786, 0, 67472.8 (83.5184%), 1.126640e-03
28.	67814.3, 524.935, 0, 67289.4 (83.2914%), 1.239304e-03
29.	67765.5, 548.727, 0, 67216.8 (83.2016%), 1.363235e-03
30.	67698, 576.78, 0, 67121.2 (83.0832%), 1.499558e-03
31.	67456.3, 590.771, 0, 66865.5 (82.7667%), 8.247569e-04
32.	67357.5, 607.778, 0, 66749.7 (82.6234%), 9.072326e-04
33.	67300.9, 626.197, 0, 66674.7 (82.5306%), 9.979559e-04
34.	67209, 647.007, 0, 66562 (82.3911%), 1.097751e-03
35.	67112.8, 669.815, 0, 66443 (82.2437%), 1.207527e-03
36.	67030.7, 695.584, 0, 66335.2 (82.1102%), 1.328279e-03
37.	66964.7, 723.733, 0, 66240.9 (81.9936%), 1.461107e-03
38.	66933, 756.304, 0, 66176.7 (81.9141%), 1.607218e-03
39.	66805, 772.968, 0, 66032 (81.735%), 8.839698e-04
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
40.	66692.4, 792.942, 0, 65899.4 (81.5709%), 9.723668e-04
41.	66670.2, 814.083, 0, 65856.1 (81.5173%), 1.069604e-03
42.	66572.2, 838.601, 0, 65733.6 (81.3657%), 1.176564e-03
43.	66570, 864.628, 0, 65705.3 (81.3306%), 1.294220e-03
44.	66513, 895.415, 0, 65617.6 (81.222%), 1.423642e-03
45.	66395.6, 910.823, 0, 65484.8 (81.0576%), 7.830032e-04
46.	66332.4, 929.375, 0, 65403 (80.9564%), 8.613036e-04
47.	66296.9, 949.192, 0, 65347.7 (80.888%), 9.474339e-04
48.	66228, 971.662, 0, 65256.3 (80.7749%), 1.042177e-03
49.	66181.8, 996.155, 0, 65185.7 (80.6874%), 1.146395e-03
E = 66181.8 (80.6874%)
Length = 24.9747
Time = 39.1043s (0.651739m)

Step 2: alpha=0.05, beta=0.05, scale=0.25
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
0.	63799, 5.6476, 0, 63793.3 (98.4367%), 2.750000e-03
1.	62206.1, 4.5564, 0, 62201.6 (95.9806%), 1.512500e-03
2.	61478.6, 5.46345, 0, 61473.1 (94.8565%), 8.318750e-04
3.	61291.3, 7.06033, 0, 61284.3 (94.5651%), 9.150625e-04
4.	61115.4, 9.05287, 0, 61106.4 (94.2906%), 1.006569e-03
5.	60958.4, 11.4717, 0, 60946.9 (94.0445%), 1.107226e-03
6.	60799.5, 14.414, 0, 60785 (93.7948%), 1.217948e-03
7.	60642.5, 17.9339, 0, 60624.5 (93.5471%), 1.339743e-03
8.	60447, 22.1425, 0, 60424.8 (93.2389%), 1.473717e-03
9.	60247.7, 27.1286, 0, 60220.6 (92.9238%), 1.621089e-03
10.	60041.4, 33.0758, 0, 60008.4 (92.5963%), 1.783198e-03
11.	59861, 40.0686, 0, 59820.9 (92.3071%), 1.961518e-03
12.	59758, 48.4517, 0, 59709.5 (92.1352%), 2.157670e-03
13.	59538.9, 57.982, 0, 59481 (91.7825%), 2.373436e-03
14.	59529.8, 63.9057, 0, 59465.9 (91.7592%), 1.305390e-03
15.	59365.3, 69.9745, 0, 59295.3 (91.4961%), 1.435929e-03
16.	59255.2, 73.5901, 0, 59181.6 (91.3206%), 7.897610e-04
17.	59165.2, 77.4935, 0, 59087.7 (91.1756%), 8.687371e-04
18.	59103.9, 81.9245, 0, 59022 (91.0742%), 9.556108e-04
19.	59029.3, 86.8711, 0, 58942.4 (90.9516%), 1.051172e-03
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
20.	58961.1, 92.4463, 0, 58868.6 (90.8377%), 1.156289e-03
21.	58874.6, 98.6993, 0, 58775.9 (90.6946%), 1.271918e-03
22.	58815.3, 105.741, 0, 58709.5 (90.5921%), 1.399110e-03
23.	58708.7, 113.631, 0, 58595.1 (90.4156%), 1.539021e-03
24.	58658.2, 122.551, 0, 58535.7 (90.3239%), 1.692923e-03
25.	58510.7, 132.524, 0, 58378.2 (90.0809%), 1.862215e-03
26.	58494.6, 138.146, 0, 58356.5 (90.0473%), 1.024218e-03
27.	58419.9, 144.337, 0, 58275.5 (89.9224%), 1.126640e-03
28.	58391.1, 151.282, 0, 58239.8 (89.8673%), 1.239304e-03
29.	58312.6, 158.982, 0, 58153.6 (89.7343%), 1.363235e-03
30.	58282.7, 167.616, 0, 58115.1 (89.6749%), 1.499558e-03
31.	58175.4, 177.185, 0, 57998.2 (89.4946%), 1.649514e-03
32.	58159.2, 182.521, 0, 57976.7 (89.4614%), 9.072326e-04
33.	58110.2, 188.41, 0, 57921.8 (89.3766%), 9.979559e-04
34.	58071.2, 194.961, 0, 57876.2 (89.3063%), 1.097751e-03
35.	58022, 202.232, 0, 57819.8 (89.2192%), 1.207527e-03
36.	57987.3, 210.327, 0, 57777 (89.1532%), 1.328279e-03
37.	57913.4, 219.329, 0, 57694.1 (89.0253%), 1.461107e-03
38.	57877.8, 229.403, 0, 57648.4 (88.9547%), 1.607218e-03
39.	57795.1, 240.586, 0, 57554.5 (88.8098%), 1.767940e-03
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
40.	57771.4, 246.83, 0, 57524.6 (88.7638%), 9.723668e-04
41.	57734.1, 253.717, 0, 57480.3 (88.6954%), 1.069604e-03
42.	57699, 261.382, 0, 57437.6 (88.6295%), 1.176564e-03
43.	57655.2, 269.884, 0, 57385.3 (88.5487%), 1.294220e-03
44.	57632, 279.354, 0, 57352.7 (88.4985%), 1.423642e-03
45.	57557.7, 289.835, 0, 57267.9 (88.3676%), 1.566006e-03
46.	57549.1, 301.593, 0, 57247.5 (88.3361%), 1.722607e-03
47.	57499.3, 314.573, 0, 57184.7 (88.2393%), 1.894868e-03
48.	57449.1, 321.822, 0, 57127.3 (88.1506%), 1.042177e-03
49.	57414.7, 329.764, 0, 57084.9 (88.0853%), 1.146395e-03
E = 57414.7 (88.0853%)
Length = 13.2401
Time = 87.4781s (1.45797m)

Step 3: alpha=0.05, beta=0.05, scale=0.5
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
0.	60797.4, 5.05434, 0, 60792.3 (96.2588%), 5.500000e-03
1.	60261.5, 4.03493, 0, 60257.5 (95.4119%), 3.781250e-04
2.	60261.4, 4.03483, 0, 60257.4 (95.4118%), 2.030945e-07
E = 60261.4 (95.4118%)
Length = 2.38056
Time = 61.4595s (1.02432m)

Step 4: alpha=0.05, beta=0.05, scale=1.0
	E, E_velocity, E_rate, E_image (E_image %), LearningRate
0.	55346.6, 1.89644, 0, 55344.7 (94.8336%), 1.375000e-03
1.	54978.9, 0.866762, 0, 54978 (94.2052%), 3.781250e-04
E = 54978.9 (94.2052%)
Length = 2.29652
Time = 304.056s (5.06761m)

Visualize registered image

The two images below should match if the registration worked successfully!



In [13]:

    
plotter.imgShow(atlas_registered)



In [19]:

    
plotter.imgShow(plotter.imgChecker(atlas_registered, img_p), vmax=2)

Quantitative evaluation

Here, we print out the Mean Squared Error between both the atlas and the observed data. As we can see, this metric decreases from the unprocessed data (first cell below this one) to the final atlas registered to our data (3rd cell below this one)



In [15]:

    
ndreg.imgMSE(sitk.Normalize(atlas), sitk.Normalize(img))









    Out[15]:





1.9355552516414318



In [16]:

    
ndreg.imgMSE(sitk.Normalize(atlas), sitk.Normalize(img_p))









    Out[16]:





0.7295863343639051



In [17]:

    
ndreg.imgMSE(sitk.Normalize(atlas_registered), sitk.Normalize(img_p))









    Out[17]:





0.416928609017568