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:

  1. select a region of interest: crop

  2. 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*')
plot multiModalRegistration

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*')
plot multiModalRegistration

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()
plot multiModalRegistration

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()
plot multiModalRegistration

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

Gallery generated by Sphinx-Gallery