Multimodal registration

A simple example to register two images acquired with different modalities

 8
 9
10
11
12
13
14
 # sphinx_gallery_thumbnail_number = 4
 import spam.DIC
 import spam.deformation
 import spam.datasets
 import numpy
 import matplotlib as mpl
 import matplotlib.pyplot as plt

Load the two images

19
20
21
22
23
24
25
26
27
 xr = spam.datasets.loadConcreteXr().astype('<f4')
 ne = spam.datasets.loadConcreteNe().astype('<f4')

 #
 #
 # plt.figure()
 # plt.imshow(xr[halfSlice, :, :])
 # plt.figure()
 # 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
36
37
38
39
40
41
42
43
44
 cropRatio = 0.1
 crop = (slice(int(cropRatio * xr.shape[0]), int((1 - cropRatio) * xr.shape[0])),
         slice(int(cropRatio * xr.shape[1]), int((1 - cropRatio) * xr.shape[1])),
         slice(int(cropRatio * xr.shape[2]), int((1 - cropRatio) * xr.shape[2])))
 cropPx = int(cropRatio * numpy.mean(xr.shape[0]))
 marginRatio = 0.1
 marginPx = int(marginRatio * numpy.mean(xr.shape[0]))
 cropWithMargin = (slice(int((cropRatio + marginRatio) * xr.shape[0]), int((1 - (cropRatio + marginRatio)) * xr.shape[0])),
                   slice(int((cropRatio + marginRatio) * xr.shape[1]), int((1 - (cropRatio + marginRatio)) * xr.shape[1])),
                   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
50
51
52
53
54
55
 bins = 128
 xrMin = xr.min()
 xrMax = xr.max()
 neMax = ne.max()
 neMin = ne.min()
 xr = numpy.array(bins * (xr - xrMin) / (xrMax - xrMin)).astype('<u1')
 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
62
63
64
65
 halfSlice = xr.shape[0] // 2
 #checker = spam.DIC.checkerBoard(xr[halfSlice], ne[halfSlice], n=3)
 #plt.figure()
 #plt.imshow(checker, cmap="Greys")
 #plt.colorbar()

We apply and initial guess transformation.

69
70
71
72
73
 PhiGuess = spam.deformation.computePhi({'t': [0.0, 0.0, 0.0], 'r': [15.0, 0.0, 0.0]})
 tmp = spam.deformation.decomposePhi(PhiGuess)
 neTmp = spam.DIC.applyPhi(ne.copy(), Phi=PhiGuess).astype('<u1')
 print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*tmp['t']))
 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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
 print("STEP 1: Get gaussian parameters")
 nPhases = 2
 gaussianParameters, jointHistogram = spam.DIC.gaussianMixtureParameters(xr[cropWithMargin],
                                                                         neTmp[cropWithMargin],
                                                                         BINS=bins,
                                                                         NPHASES=nPhases)
 plt.figure()
 tmp = jointHistogram.copy()
 tmp[jointHistogram <= 0] = numpy.nan
 tmp = numpy.log(tmp)
 plt.imshow(tmp.T, origin='lower', extent=[0.0, bins, 0.0, bins])
 plt.xlabel("x-ray grey levels")
 plt.ylabel("neutron grey levels")
 plt.colorbar()
 for gp in gaussianParameters:
     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
105
106
107
108
109
110
111
112
113
114
115
 print("STEP 2: Create phase repartition")
 phaseDiagram, actualVoxelCoverage = spam.DIC.phaseDiagram(gaussianParameters,
                                                           jointHistogram,
                                                           voxelCoverage=0.99,
                                                           BINS=bins)
 plt.figure()
 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))
 plt.xlabel("x-ray grey levels")
 plt.ylabel("neutron grey levels")
 plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
 for gp in gaussianParameters:
     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
120
121
122
123
124
125
126
127
128
 print("STEP 3: Registration")
 registration = spam.DIC.multimodalRegistration(xr, ne,
                                                phaseDiagram,
                                                gaussianParameters,
                                                BINS=bins,
                                                PhiInit=PhiGuess,
                                                verbose=True,
                                                margin=marginPx,
                                                maxIterations=50,
                                                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
135
136
 neReg = spam.DIC.applyPhi(ne, Phi=registration['Phi'])
 print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*registration['transformation']['t']))
 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
141
142
143
 checker = spam.DIC.checkerBoard(xr[halfSlice], neReg[halfSlice], n=3)
 plt.figure()
 plt.imshow(checker, cmap="Greys")
 plt.colorbar()
plot multiModalRegistration

Out:

<matplotlib.colorbar.Colorbar object at 0x7fe36ef64f50>

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
149
150
151
152
 phaseField = registration["phaseField"]
 plt.figure()
 plt.imshow(phaseField[halfSlice, :, :], vmin=-0.5, vmax=nPhases + 0.5, cmap=mpl.cm.get_cmap('Set1_r', nPhases + 1))
 plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
 plt.show()
plot multiModalRegistration

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

Gallery generated by Sphinx-Gallery