Note
Go to the end to download the full example code.
Expectation of global descriptors in 3D#
This example shows how to compute theoretical expectations of total volume (L3), surface area (L2), mean curvature measure (L1) and Euler characteristic (L0) of excursion sets as a function of the excursion’s threshold. Excrusions are define here as the subdomain of the domain where the Random Field is defined where the Random Field values are above a certain threshold.
Here the two measures are seen as function of the threshold value.
Monte Carlo results are confronted to the theory.
import spam.excursions
import spam.measurements
import matplotlib.pyplot as plt
import numpy
Compute the four theoretical expected measures#
The four measures of the excursion are computed and ploted for every thresholds
# spatial dimension
spatialDimension = 3
# the measure number 3
totalVolume = spam.excursions.expectedMesures(thresholds, 3, spatialDimension, std=std, lc=correlationLength, a=length)
# the measure number 2
totalSurface = 2.0 * spam.excursions.expectedMesures(thresholds, 2, spatialDimension, std=std, lc=correlationLength, a=length)
# the measure number 1
caliperDiameter = 0.5 * spam.excursions.expectedMesures(thresholds, 1, spatialDimension, std=std, lc=correlationLength, a=length)
# the measure number 0
eulerCharac = spam.excursions.expectedMesures(thresholds, 0, spatialDimension, std=std, lc=correlationLength, a=length)
plt.figure()
plt.xlabel("Threshold")
plt.title("Total volume")
plt.plot(thresholds, totalVolume, 'r')
plt.figure()
plt.xlabel("Threshold")
plt.title("Total surface area")
plt.plot(thresholds, totalSurface, 'r')
plt.figure()
plt.xlabel("Threshold")
plt.title("Mean caliper diametre")
plt.plot(thresholds, caliperDiameter, 'r')
plt.figure()
plt.xlabel("Threshold")
plt.title("Euler characteristic")
plt.plot(thresholds, eulerCharac, 'r')
[<matplotlib.lines.Line2D object at 0x7f11a7d4c310>]
Generate 5 realisations of the correlated Random Field#
In order to compare the theoretical values to Monte Carlo results, first 5 realisations of a correlated Random Field are generated.
# number of realisations
nRea = 5
nNodes = 100
# define the covariance
covarianceParameters = {'var': variance, 'len_scale': correlationLength}
# generate realisations
realisations = spam.excursions.simulateRandomField(
lengths=length,
nNodes=nNodes,
covarianceParameters=covarianceParameters,
dim=spatialDimension,
nRea=nRea)
Generating Field_0000... 32.23 seconds
Generating Field_0001... 32.26 seconds
Generating Field_0002... 32.27 seconds
Generating Field_0003... 32.25 seconds
Generating Field_0004... 32.24 seconds
We now show, on the left, a slice a the first realisation and on the right the corresponding excrusion for 3 thresholds values: -2, 0 and 2 which corresponds to three different topologies.
plt.figure()
plt.title("Middle slice of the first realisation")
plt.imshow(realisations[0][int(nNodes / 2), :, :])
plt.figure()
plt.title("Middle slice of the excursion with a threshold of -2")
plt.imshow(realisations[0][int(nNodes / 2), :, :] > -2)
plt.figure()
plt.title("Middle slice of the excursion with a threshold of 0")
plt.imshow(realisations[0][int(nNodes / 2), :, :] > 0)
plt.figure()
plt.title("Middle slice of the excursion with a threshold of 2")
plt.imshow(realisations[0][int(nNodes / 2), :, :] > 2)
<matplotlib.image.AxesImage object at 0x7f11a8fe2200>
Compute the four averaged measures (actually three)#
For every thresholds, three average measures (L0, L2 and L3) over all the realisations are compute and compared to the theoretical values.
# save average for every thresholds
totalVolumeMC = numpy.zeros_like(thresholdsMC)
totalSurfaceMC = numpy.zeros_like(thresholdsMC)
caliperDiameterMC = numpy.zeros_like(thresholdsMC)
eulerCharacMC = numpy.zeros_like(thresholdsMC)
# compute the aspect ratio
ar = length / float(nNodes)
# loop over the thresholds
for i, t in enumerate(thresholdsMC):
print(f"Threshold: {i + 1}/{len(thresholdsMC)}: {t:.2f}")
# loop over the realisations
for r in realisations:
totalVolumeMC[i] += spam.measurements.volume((r > t), aspectRatio=(ar, ar, ar)) / float(nRea)
totalSurfaceMC[i] += spam.measurements.surfaceArea(r, level=t, aspectRatio=(ar, ar, ar)) / float(nRea)
# commented because requires too much ressources
# caliperDiameterMC[i] += spam.measurements.totalCurvature(r, level=t, aspectRatio=(ar, ar, ar))/(float(nRea)*2.0*numpy.pi)
eulerCharacMC[i] += spam.measurements.eulerCharacteristic((r > t)) / float(nRea)
print(f"\t volume = {totalVolumeMC[i]:.2f} \t (err = {abs(totalVolume[i] - totalVolumeMC[i]) / abs(totalVolume[i]):.2f})")
print(f"\t surface = {totalSurfaceMC[i]:.2f} \t (err = {abs(totalSurface[i] - totalSurfaceMC[i]) / abs(totalSurface[i]):.2f})")
print(f"\t caliper = {caliperDiameterMC[i]:.2f} \t (err = {abs(caliperDiameter[i] - caliperDiameterMC[i]) / abs(caliperDiameter[i]):.2f})")
print(f"\t euler = {eulerCharacMC[i]:.2f} \t (err = {abs(eulerCharac[i] - eulerCharacMC[i]) / abs(eulerCharac[i]):.2f})")
Threshold: 1/20: -4.00
volume = 997.41 (err = 0.00)
surface = 627.83 (err = 0.04)
caliper = 0.00 (err = 1.00)
euler = 65.60 (err = 4.16)
Threshold: 2/20: -3.58
volume = 993.91 (err = 0.01)
surface = 662.83 (err = 0.10)
caliper = 0.00 (err = 1.00)
euler = 104.80 (err = 6.50)
Threshold: 3/20: -3.16
volume = 986.65 (err = 0.01)
surface = 727.70 (err = 0.20)
caliper = 0.00 (err = 1.00)
euler = 156.80 (err = 9.22)
Threshold: 4/20: -2.74
volume = 972.94 (err = 0.03)
surface = 838.73 (err = 0.39)
caliper = 0.00 (err = 1.00)
euler = 191.60 (err = 10.39)
Threshold: 5/20: -2.32
volume = 948.63 (err = 0.05)
surface = 1011.31 (err = 0.67)
caliper = 0.00 (err = 1.00)
euler = 185.40 (err = 9.05)
Threshold: 6/20: -1.89
volume = 909.15 (err = 0.09)
surface = 1244.58 (err = 1.05)
caliper = 0.00 (err = 1.00)
euler = 100.00 (err = 3.95)
Threshold: 7/20: -1.47
volume = 850.73 (err = 0.15)
surface = 1517.77 (err = 1.50)
caliper = 0.00 (err = 1.00)
euler = -71.80 (err = 4.25)
Threshold: 8/20: -1.05
volume = 772.27 (err = 0.23)
surface = 1781.81 (err = 1.93)
caliper = 0.00 (err = 1.00)
euler = -285.60 (err = 12.83)
Threshold: 9/20: -0.63
volume = 675.00 (err = 0.32)
surface = 1985.43 (err = 2.26)
caliper = 0.00 (err = 1.00)
euler = -473.00 (err = 18.96)
Threshold: 10/20: -0.21
volume = 563.77 (err = 0.44)
surface = 2072.43 (err = 2.40)
caliper = 0.00 (err = 1.00)
euler = -568.80 (err = 20.81)
Threshold: 11/20: 0.21
volume = 446.95 (err = 0.55)
surface = 2009.10 (err = 2.29)
caliper = 0.00 (err = 1.00)
euler = -529.40 (err = 17.94)
Threshold: 12/20: 0.63
volume = 334.53 (err = 0.67)
surface = 1795.64 (err = 1.93)
caliper = 0.00 (err = 1.00)
euler = -348.80 (err = 11.27)
Threshold: 13/20: 1.05
volume = 235.62 (err = 0.76)
surface = 1474.70 (err = 1.40)
caliper = 0.00 (err = 1.00)
euler = -80.40 (err = 3.18)
Threshold: 14/20: 1.47
volume = 155.97 (err = 0.84)
surface = 1119.30 (err = 0.82)
caliper = 0.00 (err = 1.00)
euler = 130.00 (err = 2.25)
Threshold: 15/20: 1.89
volume = 96.74 (err = 0.90)
surface = 780.59 (err = 0.26)
caliper = 0.00 (err = 1.00)
euler = 278.40 (err = 5.44)
Threshold: 16/20: 2.32
volume = 56.42 (err = 0.94)
surface = 504.87 (err = 0.18)
caliper = 0.00 (err = 1.00)
euler = 315.20 (err = 5.74)
Threshold: 17/20: 2.74
volume = 30.79 (err = 0.97)
surface = 302.61 (err = 0.51)
caliper = 0.00 (err = 1.00)
euler = 287.00 (err = 4.69)
Threshold: 18/20: 3.16
volume = 15.63 (err = 0.98)
surface = 166.48 (err = 0.73)
caliper = 0.00 (err = 1.00)
euler = 229.40 (err = 3.22)
Threshold: 19/20: 3.58
volume = 7.43 (err = 0.99)
surface = 84.76 (err = 0.86)
caliper = 0.00 (err = 1.00)
euler = 158.20 (err = 1.71)
Threshold: 20/20: 4.00
volume = 3.32 (err = 1.00)
surface = 40.42 (err = 0.94)
caliper = 0.00 (err = 1.00)
euler = 93.20 (err = 0.48)
We can now plot the theory and Monte Carlo measures over the thresholds
# plot volume
plt.figure()
plt.xlabel("Threshold")
plt.title("Total volume")
plt.plot(thresholds, totalVolume, 'r', label='Theory')
plt.plot(thresholdsMC, totalVolumeMC, '*b', label='Monte Carlo')
plt.legend()
# plot surface
plt.figure()
plt.xlabel("Threshold")
plt.title("Total surface area")
plt.plot(thresholds, totalSurface, 'r', label='Theory')
plt.plot(thresholdsMC, totalSurfaceMC, '*b', label='Monte Carlo')
plt.legend()
# plot caliper diametre
plt.figure()
plt.plot(thresholds, caliperDiameter, 'r')
plt.plot(thresholdsMC, caliperDiameterMC, '*b', label='Monte Carlo')
plt.xlabel("Threshold")
plt.title("Mean caliper diametre")
plt.legend()
# plot Euler characteristic
plt.figure()
plt.xlabel("Threshold")
plt.title("Euler characteristic")
plt.plot(thresholds, eulerCharac, 'r', label='Theory')
plt.plot(thresholdsMC, eulerCharacMC, '*b', label='Monte Carlo')
plt.legend()
plt.show()
Total running time of the script: (2 minutes 55.435 seconds)