Note
Click here to download the full example code
Multimodal registration¶
A simple example to register two images acquired with different modalities
8 # sphinx_gallery_thumbnail_number = 4
9 import spam.DIC
10 import spam.deformation
11 import spam.datasets
12 import numpy
13 import matplotlib as mpl
14 import matplotlib.pyplot as plt
Load the two images¶
19 xr = spam.datasets.loadConcreteXr().astype('<f4')
20 ne = spam.datasets.loadConcreteNe().astype('<f4')
21
22 #
23 #
24 # plt.figure()
25 # plt.imshow(xr[halfSlice, :, :])
26 # plt.figure()
27 # plt.imshow(ne[halfSlice, :, :])
Preparation of the images¶
- We start to crop the images to:
select a region of interest:
crop
keep some margin to feed the transformation:
margin
35 cropRatio = 0.1
36 crop = (slice(int(cropRatio * xr.shape[0]), int((1 - cropRatio) * xr.shape[0])),
37 slice(int(cropRatio * xr.shape[1]), int((1 - cropRatio) * xr.shape[1])),
38 slice(int(cropRatio * xr.shape[2]), int((1 - cropRatio) * xr.shape[2])))
39 cropPx = int(cropRatio * numpy.mean(xr.shape[0]))
40 marginRatio = 0.1
41 marginPx = int(marginRatio * numpy.mean(xr.shape[0]))
42 cropWithMargin = (slice(int((cropRatio + marginRatio) * xr.shape[0]), int((1 - (cropRatio + marginRatio)) * xr.shape[0])),
43 slice(int((cropRatio + marginRatio) * xr.shape[1]), int((1 - (cropRatio + marginRatio)) * xr.shape[1])),
44 slice(int((cropRatio + marginRatio) * xr.shape[2]), int((1 - (cropRatio + marginRatio)) * xr.shape[2])))
We rescale the two images between 0 and the number of bins int the joint histogram.
For sake of efficiency they have to be saved as 8bits integers (<u1
).
49 bins = 128
50 xrMin = xr.min()
51 xrMax = xr.max()
52 neMax = ne.max()
53 neMin = ne.min()
54 xr = numpy.array(bins * (xr - xrMin) / (xrMax - xrMin)).astype('<u1')
55 ne = numpy.array(bins * (ne - neMin) / (neMax - neMin)).astype('<u1')
To see the spatial incoherency we can plot them both on the same image with a checkerbord pattern, where two adjacent squares represents the two images.
61 halfSlice = xr.shape[0] // 2
62 #checker = spam.DIC.checkerBoard(xr[halfSlice], ne[halfSlice], n=3)
63 #plt.figure()
64 #plt.imshow(checker, cmap="Greys")
65 #plt.colorbar()
We apply and initial guess transformation.
69 PhiGuess = spam.deformation.computePhi({'t': [0.0, 0.0, 0.0], 'r': [15.0, 0.0, 0.0]})
70 tmp = spam.deformation.decomposePhi(PhiGuess)
71 neTmp = spam.DIC.applyPhi(ne.copy(), Phi=PhiGuess).astype('<u1')
72 print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*tmp['t']))
73 print("Rotations : {:.4f}, {:.4f}, {:.4f}".format(*tmp['r']))
Out:
Translations: 0.0000, 0.0000, 0.0000
Rotations : 15.0000, 0.0000, 0.0000
Registration algorithm¶
We first get the parameters of the fitted gaussians and the joint histogram. x axis corresponds to the first image grey levels. y axis corresponds to the second image grey levels. The Gaussian parameters are parameters of the two ellipsoids that fit the two peaks.
83 print("STEP 1: Get gaussian parameters")
84 nPhases = 2
85 gaussianParameters, jointHistogram = spam.DIC.gaussianMixtureParameters(xr[cropWithMargin],
86 neTmp[cropWithMargin],
87 BINS=bins,
88 NPHASES=nPhases)
89 plt.figure()
90 tmp = jointHistogram.copy()
91 tmp[jointHistogram <= 0] = numpy.nan
92 tmp = numpy.log(tmp)
93 plt.imshow(tmp.T, origin='lower', extent=[0.0, bins, 0.0, bins])
94 plt.xlabel("x-ray grey levels")
95 plt.ylabel("neutron grey levels")
96 plt.colorbar()
97 for gp in gaussianParameters:
98 plt.plot(gp[1], gp[2], 'b*')

Out:
STEP 1: Get gaussian parameters
im1 from 29.00 to 128.00
im2 from 23.00 to 122.00
p normalisation: 1.00
Find maxima
Min distance between maxima: 5
00 maxima found over the 2 needed with 1.00e-01 times of the total count
00 maxima found over the 2 needed with 1.00e-02 times of the total count
02 maxima found over the 2 needed with 1.00e-03 times of the total count
Maximum 1: mu1=108.00 mu2=86.00 Phi=8.63e-03
Maximum 2: mu1=109.00 mu2=67.00 Phi=3.59e-03
Phase 1: mu1=108.00 mu2=86.00 Phi=8.63e-03
Fit: a=0.18 b=0.02 c=0.04 Hessian: 0.01
Cov: 1,1=5.8691 1,2=-3.3674 2,2=30.3347
Phase 2: mu1=109.00 mu2=67.00 Phi=3.59e-03
Fit: a=0.15 b=-0.00 c=0.00 Hessian: 0.00
Cov: 1,1=6.6323 1,2=0.8708 2,2=323.2468
Then we create the phase diagram based on the joint histogram. Each peak corresponds to a phase (1 and 2). The grey background (points to far away from a peak) is ignored (0).
104 print("STEP 2: Create phase repartition")
105 phaseDiagram, actualVoxelCoverage = spam.DIC.phaseDiagram(gaussianParameters,
106 jointHistogram,
107 voxelCoverage=0.99,
108 BINS=bins)
109 plt.figure()
110 plt.imshow(phaseDiagram.T, origin='lower', extent=[0.0, bins, 0.0, bins], vmin=-0.5, vmax=nPhases + 0.5, cmap=mpl.cm.get_cmap('Set1_r', nPhases + 1))
111 plt.xlabel("x-ray grey levels")
112 plt.ylabel("neutron grey levels")
113 plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
114 for gp in gaussianParameters:
115 plt.plot(gp[1], gp[2], 'b*')

Out:
STEP 2: Create phase repartition
9-sigma: voxel coverage = 99.15 (> 99.00%)
8-sigma: voxel coverage = 98.99 (< 99.00%) -> Returning this phase diagram.
And and we use both Gaussian parameters and phase diagram as an input of the registration algorithm
119 print("STEP 3: Registration")
120 registration = spam.DIC.multimodalRegistration(xr, ne,
121 phaseDiagram,
122 gaussianParameters,
123 BINS=bins,
124 PhiInit=PhiGuess,
125 verbose=True,
126 margin=marginPx,
127 maxIterations=50,
128 deltaPhiMin=0.005)
Out:
STEP 3: Registration
Initial state LL = 13699.66
Iteration Number 000 LL = 13699.66 dPhi = 0.0000 Tr = 0.000, 0.000, 0.000 Ro = 15.000, 0.000, 0.000 Zo = 1.000, 1.000, 1.000
Iteration Number 001 LL = 13531.89 dPhi = 0.0585 Tr = -0.041, 0.003, 0.041 Ro = 15.137, 0.020, -0.010 Zo = 1.000, 0.999, 0.998
Iteration Number 002 LL = 13384.98 dPhi = 0.0550 Tr = -0.077, 0.007, 0.082 Ro = 15.282, 0.039, -0.021 Zo = 0.999, 0.997, 0.996
Iteration Number 003 LL = 13279.62 dPhi = 0.0533 Tr = -0.115, 0.011, 0.120 Ro = 15.433, 0.059, -0.034 Zo = 0.998, 0.995, 0.994
Iteration Number 004 LL = 13139.97 dPhi = 0.0531 Tr = -0.154, 0.015, 0.154 Ro = 15.594, 0.082, -0.044 Zo = 0.998, 0.994, 0.993
Iteration Number 005 LL = 13046.46 dPhi = 0.0538 Tr = -0.196, 0.019, 0.187 Ro = 15.764, 0.106, -0.058 Zo = 0.997, 0.993, 0.991
Iteration Number 006 LL = 12939.69 dPhi = 0.0549 Tr = -0.242, 0.022, 0.217 Ro = 15.938, 0.129, -0.068 Zo = 0.996, 0.992, 0.990
Iteration Number 007 LL = 12861.05 dPhi = 0.0573 Tr = -0.292, 0.026, 0.244 Ro = 16.118, 0.151, -0.078 Zo = 0.996, 0.991, 0.989
Iteration Number 008 LL = 12733.84 dPhi = 0.0576 Tr = -0.342, 0.030, 0.271 Ro = 16.307, 0.170, -0.090 Zo = 0.995, 0.991, 0.988
Iteration Number 009 LL = 12623.39 dPhi = 0.0587 Tr = -0.394, 0.034, 0.296 Ro = 16.499, 0.191, -0.105 Zo = 0.994, 0.990, 0.988
Iteration Number 010 LL = 12548.71 dPhi = 0.0623 Tr = -0.450, 0.043, 0.322 Ro = 16.698, 0.211, -0.114 Zo = 0.994, 0.990, 0.987
Iteration Number 011 LL = 12469.26 dPhi = 0.0634 Tr = -0.507, 0.053, 0.347 Ro = 16.901, 0.231, -0.125 Zo = 0.993, 0.990, 0.987
Iteration Number 012 LL = 12403.03 dPhi = 0.0638 Tr = -0.565, 0.064, 0.369 Ro = 17.103, 0.254, -0.133 Zo = 0.993, 0.990, 0.987
Iteration Number 013 LL = 12341.03 dPhi = 0.0642 Tr = -0.623, 0.075, 0.391 Ro = 17.307, 0.279, -0.141 Zo = 0.992, 0.990, 0.987
Iteration Number 014 LL = 12287.28 dPhi = 0.0613 Tr = -0.679, 0.087, 0.412 Ro = 17.506, 0.302, -0.145 Zo = 0.992, 0.990, 0.987
Iteration Number 015 LL = 12228.45 dPhi = 0.0599 Tr = -0.733, 0.100, 0.434 Ro = 17.708, 0.325, -0.147 Zo = 0.992, 0.990, 0.987
Iteration Number 016 LL = 12196.08 dPhi = 0.0552 Tr = -0.782, 0.113, 0.454 Ro = 17.904, 0.345, -0.143 Zo = 0.992, 0.990, 0.987
Iteration Number 017 LL = 12155.41 dPhi = 0.0530 Tr = -0.828, 0.127, 0.474 Ro = 18.094, 0.364, -0.140 Zo = 0.992, 0.990, 0.987
Iteration Number 018 LL = 12102.12 dPhi = 0.0485 Tr = -0.871, 0.141, 0.492 Ro = 18.273, 0.383, -0.133 Zo = 0.992, 0.991, 0.987
Iteration Number 019 LL = 12051.18 dPhi = 0.0449 Tr = -0.909, 0.156, 0.507 Ro = 18.437, 0.402, -0.125 Zo = 0.992, 0.991, 0.988
Iteration Number 020 LL = 12021.00 dPhi = 0.0424 Tr = -0.946, 0.171, 0.521 Ro = 18.586, 0.422, -0.116 Zo = 0.992, 0.991, 0.988
Iteration Number 021 LL = 12014.26 dPhi = 0.0391 Tr = -0.980, 0.184, 0.535 Ro = 18.720, 0.441, -0.108 Zo = 0.992, 0.992, 0.988
Iteration Number 022 LL = 12006.01 dPhi = 0.0350 Tr = -1.009, 0.197, 0.548 Ro = 18.841, 0.461, -0.103 Zo = 0.992, 0.992, 0.988
Iteration Number 023 LL = 11982.48 dPhi = 0.0312 Tr = -1.035, 0.208, 0.560 Ro = 18.950, 0.477, -0.096 Zo = 0.992, 0.992, 0.988
Iteration Number 024 LL = 11965.45 dPhi = 0.0272 Tr = -1.058, 0.217, 0.572 Ro = 19.047, 0.489, -0.091 Zo = 0.992, 0.993, 0.989
Iteration Number 025 LL = 11948.00 dPhi = 0.0243 Tr = -1.079, 0.224, 0.581 Ro = 19.134, 0.500, -0.088 Zo = 0.992, 0.993, 0.989
Iteration Number 026 LL = 11934.43 dPhi = 0.0221 Tr = -1.097, 0.231, 0.590 Ro = 19.210, 0.510, -0.084 Zo = 0.992, 0.993, 0.989
Iteration Number 027 LL = 11924.29 dPhi = 0.0202 Tr = -1.115, 0.238, 0.597 Ro = 19.277, 0.519, -0.081 Zo = 0.992, 0.994, 0.989
Iteration Number 028 LL = 11917.08 dPhi = 0.0175 Tr = -1.130, 0.243, 0.603 Ro = 19.336, 0.527, -0.077 Zo = 0.992, 0.994, 0.989
Iteration Number 029 LL = 11912.87 dPhi = 0.0146 Tr = -1.144, 0.246, 0.608 Ro = 19.386, 0.534, -0.074 Zo = 0.992, 0.995, 0.989
Iteration Number 030 LL = 11902.09 dPhi = 0.0130 Tr = -1.155, 0.249, 0.612 Ro = 19.428, 0.540, -0.073 Zo = 0.992, 0.995, 0.989
Iteration Number 031 LL = 11891.75 dPhi = 0.0117 Tr = -1.166, 0.251, 0.616 Ro = 19.467, 0.546, -0.072 Zo = 0.992, 0.995, 0.989
Iteration Number 032 LL = 11882.82 dPhi = 0.0096 Tr = -1.175, 0.253, 0.618 Ro = 19.501, 0.551, -0.072 Zo = 0.992, 0.995, 0.990
Iteration Number 033 LL = 11869.45 dPhi = 0.0081 Tr = -1.183, 0.254, 0.619 Ro = 19.528, 0.555, -0.072 Zo = 0.992, 0.996, 0.990
Iteration Number 034 LL = 11873.17 dPhi = 0.0064 Tr = -1.189, 0.255, 0.619 Ro = 19.552, 0.558, -0.073 Zo = 0.992, 0.996, 0.990
Iteration Number 035 LL = 11883.93 dPhi = 0.0053 Tr = -1.194, 0.256, 0.620 Ro = 19.573, 0.560, -0.074 Zo = 0.992, 0.996, 0.990
Iteration Number 036 LL = 11891.06 dPhi = 0.0045 Tr = -1.199, 0.256, 0.620 Ro = 19.590, 0.563, -0.075 Zo = 0.992, 0.996, 0.990
-> Converged
Final transformation¶
We can now apply the final transformation
134 neReg = spam.DIC.applyPhi(ne, Phi=registration['Phi'])
135 print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*registration['transformation']['t']))
136 print("Rotations : {:.4f}, {:.4f}, {:.4f}".format(*registration['transformation']['r']))
Out:
Translations: -1.1988, 0.2564, 0.6200
Rotations : 19.5900, 0.5626, -0.0752
And check the validity of the result with a checkerboard pattern mixing the two images
140 checker = spam.DIC.checkerBoard(xr[halfSlice], neReg[halfSlice], n=3)
141 plt.figure()
142 plt.imshow(checker, cmap="Greys")
143 plt.colorbar()

Out:
<matplotlib.colorbar.Colorbar object at 0x7fb7f29c96d0>
From the phase diagram a segemntation can also directly be obtained We can check that phase 1 corresponds to the mortar matrix and phase 2 to the aggregates
148 phaseField = registration["phaseField"]
149 plt.figure()
150 plt.imshow(phaseField[halfSlice, :, :], vmin=-0.5, vmax=nPhases + 0.5, cmap=mpl.cm.get_cmap('Set1_r', nPhases + 1))
151 plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
152 plt.show()

Total running time of the script: ( 0 minutes 43.110 seconds)