Note
Go to the end to download the full example code
Multimodal registration
A simple example to register two images acquired with different modalities
8 import matplotlib as mpl
9 import matplotlib.pyplot as plt
10 import numpy
11 import spam.datasets
12 import spam.deformation
13
14 # sphinx_gallery_thumbnail_number = 4
15 import spam.DIC
16 import spam.helpers
Load the two images
21 xr = spam.datasets.loadConcreteXr().astype("<f4")
22 ne = spam.datasets.loadConcreteNe().astype("<f4")
23
24 #
25 #
26 # plt.figure()
27 # plt.imshow(xr[halfSlice, :, :])
28 # plt.figure()
29 # 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
37 cropRatio = 0.1
38 crop = (
39 slice(int(cropRatio * xr.shape[0]), int((1 - cropRatio) * xr.shape[0])),
40 slice(int(cropRatio * xr.shape[1]), int((1 - cropRatio) * xr.shape[1])),
41 slice(int(cropRatio * xr.shape[2]), int((1 - cropRatio) * xr.shape[2])),
42 )
43 cropPx = int(cropRatio * numpy.mean(xr.shape[0]))
44 marginRatio = 0.1
45 marginPx = int(marginRatio * numpy.mean(xr.shape[0]))
46 cropWithMargin = (
47 slice(
48 int((cropRatio + marginRatio) * xr.shape[0]),
49 int((1 - (cropRatio + marginRatio)) * xr.shape[0]),
50 ),
51 slice(
52 int((cropRatio + marginRatio) * xr.shape[1]),
53 int((1 - (cropRatio + marginRatio)) * xr.shape[1]),
54 ),
55 slice(
56 int((cropRatio + marginRatio) * xr.shape[2]),
57 int((1 - (cropRatio + marginRatio)) * xr.shape[2]),
58 ),
59 )
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
).
64 bins = 128
65 xrMin = xr.min()
66 xrMax = xr.max()
67 neMax = ne.max()
68 neMin = ne.min()
69 xr = numpy.array(bins * (xr - xrMin) / (xrMax - xrMin)).astype("<u1")
70 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.
76 halfSlice = xr.shape[0] // 2
77 # checker = spam.helpers.checkerBoard(xr[halfSlice], ne[halfSlice], n=3)
78 # plt.figure()
79 # plt.imshow(checker, cmap="Greys")
80 # plt.colorbar()
We apply and initial guess transformation.
84 PhiGuess = spam.deformation.computePhi({"t": [0.0, 0.0, 0.0], "r": [15.0, 0.0, 0.0]})
85 tmp = spam.deformation.decomposePhi(PhiGuess)
86 neTmp = spam.DIC.applyPhi(ne.copy(), Phi=PhiGuess).astype("<u1")
87 print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*tmp["t"]))
88 print("Rotations : {:.4f}, {:.4f}, {:.4f}".format(*tmp["r"]))
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.
98 print("STEP 1: Get gaussian parameters")
99 nPhases = 2
100 gaussianParameters, jointHistogram = spam.DIC.gaussianMixtureParameters(xr[cropWithMargin], neTmp[cropWithMargin], BINS=bins, NPHASES=nPhases)
101 plt.figure()
102 tmp = jointHistogram.copy()
103 tmp[jointHistogram <= 0] = numpy.nan
104 tmp = numpy.log(tmp)
105 plt.imshow(tmp.T, origin="lower", extent=[0.0, bins, 0.0, bins])
106 plt.xlabel("x-ray grey levels")
107 plt.ylabel("neutron grey levels")
108 plt.colorbar()
109 for gp in gaussianParameters:
110 plt.plot(gp[1], gp[2], "b*")

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
/usr/local/lib/python3.8/site-packages/spam/DIC/multimodal.py:675: RuntimeWarning: overflow encountered in cast
field[nx, ny] = float(GLOBALphi) * numpy.exp(-computeLambda(a, b, c, x, GLOBALmu1, y, GLOBALmu2))
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).
116 print("STEP 2: Create phase repartition")
117 phaseDiagram, actualVoxelCoverage = spam.DIC.phaseDiagram(gaussianParameters, jointHistogram, voxelCoverage=0.99, BINS=bins)
118 plt.figure()
119 plt.imshow(
120 phaseDiagram.T,
121 origin="lower",
122 extent=[0.0, bins, 0.0, bins],
123 vmin=-0.5,
124 vmax=nPhases + 0.5,
125 cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),
126 )
127 plt.xlabel("x-ray grey levels")
128 plt.ylabel("neutron grey levels")
129 plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
130 for gp in gaussianParameters:
131 plt.plot(gp[1], gp[2], "b*")

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.
/builds/ttk/spam/examples/DIC/plot_multiModalRegistration.py:125: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),
And and we use both Gaussian parameters and phase diagram as an input of the registration algorithm
135 print("STEP 3: Registration")
136 registration = spam.DIC.multimodalRegistration(
137 xr,
138 ne,
139 phaseDiagram,
140 gaussianParameters,
141 BINS=bins,
142 PhiInit=PhiGuess,
143 verbose=True,
144 margin=marginPx,
145 maxIterations=50,
146 deltaPhiMin=0.005,
147 )
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
153 neReg = spam.DIC.applyPhi(ne, Phi=registration["Phi"])
154 print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*registration["transformation"]["t"]))
155 print("Rotations : {:.4f}, {:.4f}, {:.4f}".format(*registration["transformation"]["r"]))
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
159 checker = spam.helpers.checkerBoard(xr[halfSlice], neReg[halfSlice], n=3)
160 plt.figure()
161 plt.imshow(checker, cmap="Greys")
162 plt.colorbar()

<matplotlib.colorbar.Colorbar object at 0x7fb1c0ed31c0>
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
167 phaseField = registration["phaseField"]
168 plt.figure()
169 plt.imshow(
170 phaseField[halfSlice, :, :],
171 vmin=-0.5,
172 vmax=nPhases + 0.5,
173 cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),
174 )
175 plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
176 plt.show()

/builds/ttk/spam/examples/DIC/plot_multiModalRegistration.py:173: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),
Total running time of the script: ( 0 minutes 30.585 seconds)