Commit 25734c8f authored by Edward Andò's avatar Edward Andò
Browse files

fix for mmr example

parent a4c74b56
Pipeline #68664 passed with stages
in 28 minutes and 14 seconds
......@@ -5,7 +5,7 @@ Multimodal registration
A simple example to register two images acquired with different modalities
"""
# sphinx_gallery_thumbnail_number = 5
# sphinx_gallery_thumbnail_number = 4
import spam.DIC
import spam.deformation
import spam.datasets
......@@ -36,10 +36,12 @@ 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])))
margin = 10
cropWithMargin = (slice(int(cropRatio * xr.shape[0] + margin), int((1 - cropRatio) * xr.shape[0] - margin)),
slice(int(cropRatio * xr.shape[1] + margin), int((1 - cropRatio) * xr.shape[1] - margin)),
slice(int(cropRatio * xr.shape[2] + margin), int((1 - cropRatio) * xr.shape[2] - margin)))
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.
......@@ -57,14 +59,14 @@ ne = numpy.array(bins * (ne - neMin) / (neMax - neMin)).astype('<u1')
# where two adjacent squares represents the two images.
halfSlice = xr.shape[0] // 2
checker = spam.DIC.checkerBoard(xr[halfSlice], ne[halfSlice], n=3)
plt.figure()
plt.imshow(checker, cmap="Greys")
plt.colorbar()
#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.
PhiGuess = spam.deformation.computePhi({'t': [0.0, 0.0, 0.0], 'r': [-15.0, 0.0, 0.0]})
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']))
......@@ -80,7 +82,10 @@ print("Rotations : {:.4f}, {:.4f}, {:.4f}".format(*tmp['r']))
# The Gaussian parameters are parameters of the two ellipsoids that fit the two peaks.
print("STEP 1: Get gaussian parameters")
nPhases = 2
gaussianParameters, jointHistogram = spam.DIC.gaussianMixtureParameters(xr[cropWithMargin], neTmp[cropWithMargin], BINS=bins, NPHASES=nPhases)
gaussianParameters, jointHistogram = spam.DIC.gaussianMixtureParameters(xr[cropWithMargin],
neTmp[cropWithMargin],
BINS=bins,
NPHASES=nPhases)
plt.figure()
tmp = jointHistogram.copy()
tmp[jointHistogram <= 0] = numpy.nan
......@@ -97,7 +102,10 @@ for gp in gaussianParameters:
# Each peak corresponds to a phase (1 and 2).
# The grey background (points to far away from a peak) is ignored (0).
print("STEP 2: Create phase repartition")
phaseDiagram, actualVoxelCoverage = spam.DIC.phaseDiagram(gaussianParameters, jointHistogram, voxelCoverage=0.8, BINS=bins)
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")
......@@ -109,8 +117,15 @@ for gp in gaussianParameters:
#############################################################
# And and we use both Gaussian parameters and phase diagram as an input of the registration algorithm
print("STEP 3: Registration")
registration = spam.DIC.multimodalRegistration(xr[crop], ne[crop], phaseDiagram, gaussianParameters,
BINS=bins, PhiInit=PhiGuess, verbose=True, margin=margin)
registration = spam.DIC.multimodalRegistration(xr, ne,
phaseDiagram,
gaussianParameters,
BINS=bins,
PhiInit=PhiGuess,
verbose=True,
margin=marginPx,
maxIterations=50,
deltaPhiMin=0.005)
##########################################################
# Final transformation
......@@ -122,7 +137,7 @@ print("Rotations : {:.4f}, {:.4f}, {:.4f}".format(*registration['transformatio
###########################################################
# And check the validity of the result with a checkerboard pattern mixing the two images
checker = spam.DIC.checkerBoard(xr[xr.shape[0] // 2], neReg[neReg.shape[0] // 2], n=3)
checker = spam.DIC.checkerBoard(xr[halfSlice], neReg[halfSlice], n=3)
plt.figure()
plt.imshow(checker, cmap="Greys")
plt.colorbar()
......@@ -132,6 +147,6 @@ plt.colorbar()
# We can check that phase 1 corresponds to the mortar matrix and phase 2 to the aggregates
phaseField = registration["phaseField"]
plt.figure()
plt.imshow(phaseField[phaseField.shape[0] // 2, :, :], vmin=-0.5, vmax=nPhases + 0.5, cmap=mpl.cm.get_cmap('Set1_r', nPhases + 1))
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()
......@@ -130,7 +130,7 @@ def loadDEMtouchingGrainsAndBranchVector():
print(IOError("File %s not found" % data_file_path))
return 0
return loadtxt("DEM/spheres.txt")
def loadDEMspheresMMR():
filepath = "DEM/SpheresMMR.txt"
data_file_path = 'data/' + filepath
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment