Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

Commit 595c3380 authored by Alessandro's avatar Alessandro
Browse files

First day

parent be865c3f
Pipeline #50970 failed with stages
in 12 minutes and 47 seconds
......@@ -78,9 +78,10 @@ print("Rotations : {:.4f}, {:.4f}, {:.4f}".format(*tmp['r']))
# 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.
print("STEP 1: Get gaussian parameters")
print("\nSTEP 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, verbose=True, GRAPHS = True, INTERACTIVE = True)
exit()
plt.figure()
tmp = jointHistogram.copy()
tmp[jointHistogram <= 0] = numpy.nan
......@@ -96,7 +97,7 @@ for gp in gaussianParameters:
# 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).
print("STEP 2: Create phase repartition")
print("\nSTEP 2: Create phase repartition")
phaseDiagram, actualVoxelCoverage = spam.DIC.phaseDiagram(gaussianParameters, jointHistogram, voxelCoverage=0.8, 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))
......@@ -108,7 +109,7 @@ 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")
print("\nSTEP 3: Registration")
registration = spam.DIC.multimodalRegistration(xr[crop], ne[crop], phaseDiagram, gaussianParameters,
BINS=bins, PhiInit=PhiGuess, verbose=True, margin=margin)
......
......@@ -219,8 +219,8 @@ for i, (bin, nPhases) in enumerate(zip(bins, args.PHASES)):
# except:
gaussianParameters, jointHistogram = spam.DIC.gaussianMixtureParameters(im1[cropWithMargin],
im2def[cropWithMargin],
BINS=args.JOINT_HISTO_BINS,
NPHASES=nPhases,
bins=args.JOINT_HISTO_BINS,
nPhases=nPhases,
im1threshold=args.IM1_THRESHOLD, im2threshold=args.IM2_THRESHOLD,
distanceMaxima=distanceMaxima,
fitDistance=args.FIT_DISTANCE,
......
......@@ -88,8 +88,9 @@ def multimodalRegistration(im1, im2, phaseDiagram, gaussianParameters, im1mask=N
Order of the greylevel interpolation for applying F to im1 when correlating. Reccommended value is 3, but you can get away with 1 for faster calculations. Default = 3
verbose : bool, optional
Get to know what the function is really thinking, recommended for debugging only. Default = False
Get to know what the function is really thinking, recommended for debugging only.
Default = False
GRAPHS : bool, optional
Pop up a window showing the image differences (im1-im2) as im1 is progressively deformed.
......@@ -488,49 +489,71 @@ def multimodalRegistration(im1, im2, phaseDiagram, gaussianParameters, im1mask=N
# helper functions for correlation with different modalities
################################################################
def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2threshold=0, distanceMaxima=None, fitDistance=None, GRAPHS=False, INTERACTIVE=False, sliceAxis=0, rootPath=".", suffix=""):
def gaussianMixtureParameters(im1, im2, bins=64, nPhases=2, im1threshold=0, im2threshold=0, distanceMaxima=None, fitDistance=None, GRAPHS=False, INTERACTIVE=False, verbose=False, sliceAxis=0, rootPath=".", suffix=""):
"""
This function, given two different modality images, builds the joint histogram in BINS bins,
and fits NPHASES peaks with bivariate Gaussian distributions.
This function, given two different modality images, builds the joint histogram in bins×bins,
and fits nPhases peaks with bivariate Gaussian distributions.
Parameters
----------
im1 : 3D numpy array of floats
One image, values should be rescaled to 0:BIN-1
One image, values should be rescaled to 0:bins-1
im2 : 3D numpy array of floats
Another image, values should be rescaled to 0:BIN-1
Another image, values should be rescaled to 0:bins-1
BINS : int, optional
bins : int, optional
Number of bins for the joint histogram, recommend 2^x
Default = 64
NPHASES : int, optional
Number of phases to automatically fit
nPhases : int, optional
Number of phases to automatically fit with a 2D-Gaussian function
Default = 2
im1threshold : float, optional
im1threshold : float, optional
Greyvalue threshold above which the bivariate histogram fits the 2D-Gaussian functions for im1
Default = 0
im2threshold : float, optional
Greyvalue threshold above which the bivariate histogram fits the 2D-Gaussian functions for im2
Default = 0
distanceMaxima : float, optional
Minimum number of points (rescaled greyvalues) separating peaks in a region of 2 * distanceMaxima + 1 (i.e. peaks are separated by at least distanceMaxima)
Default = bins/25.0
fitDistance : float, optional
Radius of the area within which the 2D-Gaussian is fitted
Default = bins/2.0
GRAPHS : bool, optional
Interactively plots graphs
Default = False
INTERACTIVE : bool, optional
Activate interactive mode
Default = False
verbose : bool, optional
Get to know what the function is really thinking, recommended for debugging only.
Default = False
sliceAxis=0
Axis of the image in which the plot evolves if INTERACTIVE=True
Default = 0 (i.e. the Z-axis is shown)
rootPath="."
Path in which output plots are saved if INTERACTIVE=False
Defaults to working directory
suffix=""
Suffix of the output plots if INTERACTIVE=False
Defaults to file name, e.g. "GaussianMixture_jointHistogram-i-j.png"
Returns
-------
gaussianParameters : list of lists
List of NPHASES components, containing the parameters of the bivariate Gaussian fit for each phase:
List of nPhases components, containing the parameters of the bivariate Gaussian fit for each phase:
[0] = GLOBALphi -- Normlised "Height" of the fit
[1] = GLOBALmu1 -- Mean in F of bivairate Gaussian
[2] = GLOBALmu2 -- Mean in G of bivairate Gaussian
......@@ -539,7 +562,7 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
[5] = c --
p : BINSxBINS 2D numpy array of floats
p : binsxbins 2D numpy array of floats
The normalised joint histogram, the sum of which will be =1 if all of your image values are 0:BIN-1
"""
......@@ -586,18 +609,20 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
im2max = im2.max()
# f,g discretisation
X = numpy.linspace(0, BINS-1, BINS)
Y = numpy.linspace(0, BINS-1, BINS)
X = numpy.linspace(0, bins-1, bins)
Y = numpy.linspace(0, bins-1, bins)
print("\tim1 from {:.2f} to {:.2f}".format(im1min, im1max))
print("\tim2 from {:.2f} to {:.2f}".format(im2min, im2max))
if verbose:
print("\tim1 from {:.2f} to {:.2f}".format(im1min, im1max))
print("\tim2 from {:.2f} to {:.2f}".format(im2min, im2max))
# Calculate probability distribution p(f,g)
p, _, _ = numpy.histogram2d(im1.ravel(), im2.ravel(), bins=BINS, range=[[0.0, BINS], [0.0, BINS]], normed=False, weights=None)
p, _, _ = numpy.histogram2d(im1.ravel(), im2.ravel(), bins=bins, range=[[0.0, bins], [0.0, bins]], normed=False, weights=None)
p /= float(len(im1.ravel()))
p_sum = p.sum()
print("\tp normalisation: {:.2f}".format(p_sum))
if verbose:
print("\tp normalisation: {:.2f}".format(p_sum))
# write joint histogram in a dat file for tikz
# if DATA:
......@@ -614,52 +639,62 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
# f.write(string)
# Get disordered maxima of the joint histogram
print("\tFind maxima")
if verbose:
print("\tFind maxima")
if distanceMaxima is None:
distanceMaxima = int(BINS/25.0)
print("\t\tMin distance between maxima: {:.0f}".format(distanceMaxima))
distanceMaxima = int(bins/25.0)
if verbose:
print("\t\tMin distance between maxima: {:.0f}".format(distanceMaxima))
maxima = skimage.feature.peak_local_max(p, min_distance=int(distanceMaxima))
nMax = len(maxima)
if(nMax < NPHASES):
print("\t\t{} phases asked but only {} maxima...".format(NPHASES, nMax))
NPHASES = min(NPHASES, nMax)
if(nMax < nPhases):
print("\n\nspam.DIC.gaussianMixtureParameters():{} phases asked but only {} maxima...".format(nPhases, nMax))
nPhases = min(nPhases, nMax)
# print "\t\t Number of maxima: {:2}".format(nMax)
if nMax == 0: # no cover
print("In gaussian fit: no maxium found... Stoping...")
print("\n\nspam.DIC.gaussianMixtureParameters(): no maxium found... Exiting...")
exit()
# Organise maxima into muF, muG, p(f,g) the sort at take the maximum
maxTmp = numpy.zeros((nMax, 3))
for i, (f, g) in enumerate(maxima):
#print(i,f,g)
if f > im1threshold or g > im2threshold:
maxTmp[i] = [f, g, p[f, g]]
# print("\t\t max {:2} --- f: {:.2f} g: {:.2f} hits: {}".format(i+1,*maxTmp[i]))
#print("\t\t max {:2} --- f: {:.2f} g: {:.2f} hits: {}".format(i+1,*maxTmp[i]))
nMax = 0
percentage = 0.1
while nMax < NPHASES:
while nMax < nPhases:
maxSorted = [m for m in maxTmp[maxTmp[:, 2].argsort()[::-1]] if float(m[2]) > percentage*float(p_sum)]
nMax = len(maxSorted)
print("\t\t{:02d} maxima found over the {} needed with {:.2e} times of the total count".format(nMax, NPHASES, percentage))
if verbose:
print("\t\t{:02d} maxima found over the {} needed with {:.2e} times of the total count".format(nMax, nPhases, percentage))
percentage /= 10.0
for i, (GLOBALmu1, GLOBALmu2, GLOBALphi) in enumerate(maxSorted):
print("\t\tMaximum {}:\t mu1={:.2f}\t mu2={:.2f}\t Phi={:.2e}".format(i+1, GLOBALmu1, GLOBALmu2, GLOBALphi))
print("")
if verbose:
for i, (GLOBALmu1, GLOBALmu2, GLOBALphi) in enumerate(maxSorted):
print("\t\tMaximum {}:\t mu1={:.2f}\t mu2={:.2f}\t Phi={:.2e}".format(i+1, GLOBALmu1, GLOBALmu2, GLOBALphi))
print("")
# output of the function: list of fitting gaussian parameters
# size NPHASES x 6, the 6 parameters being GLOBALlogPhi, GLOBALmu1, GLOBALmu2, a, b, c
gaussianParameters = numpy.zeros((NPHASES, 6)).astype('<f4')
# size nPhases x 6, the 6 parameters being GLOBALlogPhi, GLOBALmu1, GLOBALmu2, a, b, c
gaussianParameters = numpy.zeros((nPhases, 6)).astype('<f4')
# loop over phases
maxEllipsoid = numpy.zeros_like(p)
for iPhase in range(NPHASES):
for iPhase in range(nPhases):
GLOBALmu1, GLOBALmu2, GLOBALphi = maxSorted[iPhase]
print("\tPhase {:2}:\t\t mu1={:.2f}\t mu2={:.2f}\t Phi={:.2e}".format(iPhase+1, GLOBALmu1, GLOBALmu2, GLOBALphi))
if verbose:
print("\tPhase {:2}:\t\t mu1={:.2f}\t mu2={:.2f}\t Phi={:.2e}".format(iPhase+1, GLOBALmu1, GLOBALmu2, GLOBALphi))
if fitDistance is None:
fitDistance = BINS/2.0
fitDistance = bins/2.0
# fit the probability distribution p(f,g) with an gaussian ellipsoids
pFit = p.copy()
......@@ -674,12 +709,16 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
pFit[nf, ng] = 0.0
(a, b, c), _ = curve_fit(gaussian2Delliptical, (X, Y), pFit.ravel(), p0=(1, 1, 1))
print("\t\tFit:\t\t a={:.2f}\t b={:.2f}\t c={:.2f}\t Hessian: {:.2f}".format(a, b, c, a*c-b**2))
if verbose:
print("\t\tFit:\t\t a={:.2f}\t b={:.2f}\t c={:.2f}\t Hessian: {:.2f}".format(a, b, c, a*c-b**2))
while a*c-b**2 < 0:
print("\t\t\tWarning: Hessian < 0")
print("\t\t\t The gaussian doesn't have a local maximum.")
if verbose:
print("\t\t\tWarning: Hessian < 0")
print("\t\t\t The gaussian doesn't have a local maximum.")
fitDistance /= 2.0
print("\t\t\t Try with fit distance = {} ".format(fitDistance))
if verbose:
print("\t\t\t Try with fit distance = {} ".format(fitDistance))
for nf in range(pFit.shape[0]):
for ng in range(pFit.shape[1]):
......@@ -691,11 +730,13 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
pFit[nf, ng] = 0.0
(a, b, c), _ = curve_fit(gaussian2Delliptical, (X, Y), pFit.ravel(), p0=(1, 1, 1))
print("\t\tFit:\t\t a={:.2f}\t b={:.2f}\t c={:.2f}\t Hessian: {:.2f}".format(a, b, c, a*c-b**2))
if verbose:
print("\t\tFit:\t\t a={:.2f}\t b={:.2f}\t c={:.2f}\t Hessian: {:.2f}".format(a, b, c, a*c-b**2))
# compute covariance function
cov = scipy.linalg.inv(numpy.array([[a, b], [b, c]]))
print("\t\tCov:\t\t 1,1={:.4f}\t 1,2={:.4f}\t 2,2={:.4f}".format(cov[0, 0], cov[1, 0], cov[1, 1]))
if verbose:
print("\t\tCov:\t\t 1,1={:.4f}\t 1,2={:.4f}\t 2,2={:.4f}".format(cov[0, 0], cov[1, 0], cov[1, 1]))
gaussianParameters[iPhase, 0] = GLOBALphi
gaussianParameters[iPhase, 1] = GLOBALmu1
......@@ -715,22 +756,22 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
plt.title("im1 (from {:.2f} to {:.2f})".format(im1.min(), im1.max()))
plt.axis('off')
if sliceAxis == 0:
plt.imshow(im1[im1.shape[0]//2, :, :], vmin=0.0, vmax=BINS)
plt.imshow(im1[im1.shape[0]//2, :, :], vmin=0.0, vmax=bins)
elif sliceAxis == 1:
plt.imshow(im1[:, im1.shape[1]//2, :], vmin=0.0, vmax=BINS)
plt.imshow(im1[:, im1.shape[1]//2, :], vmin=0.0, vmax=bins)
elif sliceAxis == 2:
plt.imshow(im1[:, :, im1.shape[2]//2], vmin=0.0, vmax=BINS)
plt.imshow(im1[:, :, im1.shape[2]//2], vmin=0.0, vmax=bins)
plt.colorbar()
plt.subplot(222)
plt.title("im2 (from {:.2f} to {:.2f})".format(im2.min(), im2.max()))
plt.axis('off')
if sliceAxis == 0:
plt.imshow(im2[im2.shape[0]//2, :, :], vmin=0.0, vmax=BINS)
plt.imshow(im2[im2.shape[0]//2, :, :], vmin=0.0, vmax=bins)
elif sliceAxis == 1:
plt.imshow(im2[:, im2.shape[1]//2, :], vmin=0.0, vmax=BINS)
plt.imshow(im2[:, im2.shape[1]//2, :], vmin=0.0, vmax=bins)
elif sliceAxis == 2:
plt.imshow(im2[:, :, im2.shape[2]//2], vmin=0.0, vmax=BINS)
plt.imshow(im2[:, :, im2.shape[2]//2], vmin=0.0, vmax=bins)
plt.colorbar()
plt.subplot(223)
......@@ -738,7 +779,7 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
tmp[p <= 0] = numpy.nan
tmp = numpy.log(tmp)
plt.title("Log Probability log(p(im1,im2))")
plt.imshow(tmp.T, origin='low', extent=[0.0, BINS, 0.0, BINS])
plt.imshow(tmp.T, origin='low', extent=[0.0, bins, 0.0, bins])
for gp in maxSorted:
plt.plot(gp[0], gp[1], 'b*')
plt.plot(GLOBALmu1, GLOBALmu2, 'r*')
......@@ -754,7 +795,7 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
tmp[currentEllipsoid <= 0] = numpy.nan
tmp = numpy.log(tmp)
plt.title("Gaussian ellipsoid")
plt.imshow(tmp.T, origin='low', extent=[0.0, BINS, 0.0, BINS])
plt.imshow(tmp.T, origin='low', extent=[0.0, bins, 0.0, bins])
plt.plot(GLOBALmu1, GLOBALmu2, 'r*')
plt.xlabel("im1")
plt.ylabel("im2")
......@@ -797,8 +838,8 @@ def gaussianMixtureParameters(im1, im2, BINS=64, NPHASES=2, im1threshold=0, im2t
# plt.clf()
# tmp = maxEllipsoid.copy()
# plt.title("Maximum gaussian ellispoid")
# X = numpy.linspace(0, 1, BINS)
# Y = numpy.linspace(0, 1, BINS)
# X = numpy.linspace(0, 1, bins)
# Y = numpy.linspace(0, 1, bins)
# CS = plt.contour(X, Y, tmp.T, levels=levels)
# for gp in maxSorted:
# plt.plot(gp[0], gp[1], 'b*')
......
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