Coverage for /usr/local/lib/python3.8/site-packages/spam/measurements/covariance.py: 100%
63 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1# Library of SPAM functions for measuring the covariance
2# Copyright (C) 2020 SPAM Contributors
3#
4# This program is free software: you can redistribute it and/or modify it
5# under the terms of the GNU General Public License as published by the Free
6# Software Foundation, either version 3 of the License, or (at your option)
7# any later version.
8#
9# This program is distributed in the hope that it will be useful, but WITHOUT
10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
12# more details.
13#
14# You should have received a copy of the GNU General Public License along with
15# this program. If not, see <http://www.gnu.org/licenses/>.
18"""
19Measure correlation of a given field with two different approach.
21View examples in the gallery.
23How to import
24-------------
25>>> import spam.measurements
26>>> spam.measurements.covarianceAlongAxis(im, d)
27>>> spam.measurements.covarianceSubPixel(im)
30"""
32import numpy
33import spam.helpers
36def covarianceAlongAxis(im, d, mask=None, axis=[0, 1, 2]):
37 """
38 Compute the covariance function of 2D or 3D images along specific axis.
40 This function computes the covariance of a field (grey values or binary) with or without a mask.
41 The covariance is computed at integer distances (in voxels). No interpolation is made but specific direction can be asked.
43 Parameters
44 ----------
45 im : array, string
46 The image as a name of an input file (grey or binary) or directly an array
48 d : array
49 The list of distances (in voxels) considered to compute the covariance. It can't be bigger than the size of image and it has to be integers.
51 mask : array, optional
52 The name of the mask file (binary, 1 for phase of interest) or directly an array
54 axis : array, default=[0,1,2]
55 The list of axis in wich the direction is computed.
57 Returns
58 -------
59 c : array
60 A 2D array with list of values of the covariance at the different distances for each axis.
62 Examples
63 --------
65 >>> import numpy
66 >>> import tifffile
67 >>> im = tifffile.imread( "snow.tif" )
68 This image of size 100x100x100 is a field of grey values
69 >>> d = numpy.arange( 0, )
70 array([ 0, 1, 2, 3, 4])
71 >>> c = spam.measurements.covarianceAlongAxis( im, d )
72 array([[86833030., 74757580., 53643410., 35916468.], # axis 0
73 . [86833030., 76282920., 56910720., 39792388.], # axis 1
74 . [86833030., 76410500., 57040076., 39938920.]], dtype=float32) # axis 2
76 """
78 # convert scalar into array of size 1
79 if isinstance(d, int):
80 d = numpy.array([d])
81 if type(d[0]) == float:
82 print(
83 "spam.measurements.covariance.covarianceAlongAxis: d={}. Should be a list of integers.".format(
84 type(d[0])
85 )
86 )
87 print("exit function.")
88 return -1
89 if any([max(d) >= im.shape[i] for i in axis]):
90 print(
91 "spam.measurements.covariance.covarianceAlongAxis: max(d)={}. Should be smaller than the image.".format(
92 max(d)
93 )
94 )
95 print("exit function.")
96 return -1
98 # Step 0: apply mask
99 if mask is not None:
100 im = numpy.multiply(im, mask)
102 # Step 1: Calculate expectation and variance
103 if mask is not None:
104 E = (
105 numpy.mean(im, dtype=numpy.float32)
106 * numpy.size(mask)
107 / float(numpy.sum(mask))
108 )
109 # V = ((numpy.mean(numpy.multiply(im, im), dtype=numpy.float64)*numpy.size(mask)/float(numpy.sum(mask)))-E*E)
110 else:
111 E = numpy.mean(im, dtype=numpy.float32)
112 # V = numpy.var(im, dtype=numpy.float32)
114 # Step 2: Compute covariance c(d)
115 axis = [axis] if isinstance(axis, int) else axis
116 c = numpy.zeros((len(axis), len(d)), dtype=numpy.float32)
117 for j, a in enumerate(axis):
118 for (i, x) in enumerate(d):
119 if mask is not None:
120 # Step 2.1: Take the effectif part and compute the pairs of number
121 mask_eff = numpy.logical_and(
122 mask, spam.helpers.singleShift(mask, x, a, sub=False)
123 )
124 n = numpy.sum(mask_eff)
125 # Step 2.2: multiply the image
126 im_multi = numpy.multiply(im, spam.helpers.singleShift(im, x, a, sub=0))
127 im_multi = numpy.multiply(im_multi, mask_eff, dtype=numpy.float32)
128 else:
129 # Step 2.1: Compute the pairs of numbers
130 n = (im.shape[a] - x) * numpy.prod(
131 [s for i, s in enumerate(im.shape) if i != a]
132 ) #
133 # Step 2.2: Multiply the image
134 im_multi = numpy.multiply(
135 im, spam.helpers.singleShift(im, x, a, sub=0), dtype=numpy.float32
136 )
137 # n_multi = numpy.sum(im_multi)
138 # Step 2.3: Compute covariance
139 c[j, i] = float(numpy.sum(im_multi)) / float(n) - E * E if n > 0 else 0.0
141 return c
144def covarianceSubPixel(im, distance=25, step=10, normType=None, n_threads=1):
145 """Compute the covariance function of 2D or 3D images.
147 The correlation function is computed on the zero mean image.
148 The code behind -- for pre-allocation reasons -- works on a number of unique distances,
149 which in 3D is closely unders=estimated by the square of the distance asked for.
150 To change this a pre-allocated precision (i.e., 0.2 px) and distance could be used in the C-code.
151 A sub pixel interpolation of the grey values in made.
154 Parameters
155 ----------
156 im : array
157 The image. It is automatically converted into a 32 bit float ``<f4``.
159 distance : int, optional
160 The approximate distance in pixels over which to make this measurement.
162 step : int, optional
163 The step parameter, which is how many pixels to jump when moving the pixel of interest
164 when calculating the correaltion function (this is the best way to save time).
165 The smaller the step the longer the computational time and the better the measure of each value of the covariance.
167 normType : string, optional
168 Select the of normalisation for the covariance output.
169 ``normType="None"``: no normalisation,
170 ``normType="variance"``: divide the covariance by the variance of the image,
171 ``normType="first"``: divide the covariance by its value at 0.
173 n_threads : int, optional
174 Number of threads for the c++ ``#pragma`` command using ``openmp``.
176 Returns
177 -------
178 d : array
179 The list of distances in pixel where the correlation is computed.
181 c : array
182 The list of calues of the covariance at the different distances.
184 Examples
185 --------
187 >>> import numpy
188 >>> import tifffile
189 >>> im = tifffile.imread( "snow.tif" )
190 This image of size 100x100x100 is a field of grey values
191 >>> d = numpy.arange( 0, 4 )
192 array([ 0, 1, 2, 3, 4])
193 >>> c = spam.measurements.covarianceSubPixel( im, distance=2 )
194 (array([0. , 1. , 1.41421354, 1.73205078]), # distances
195 array( [89405248. , 77995006. , 69534446. , 62701681.])) # covariance
197 Warning
198 -------
199 The multithread version is bugged... openmp gives random values.
200 """
201 import spam.measurements.measurementsToolkit as mtk
203 # Make sure the image is in the right format for us:
204 im = im.astype("<f4")
206 # If 2D image pad with false 1st dimension
207 if len(im.shape) == 2:
208 # Make them degraded 3D images
209 im = im.reshape((1,) + im.shape)
211 # allocate memory
212 n_unique_distances = distance**2
213 output = numpy.zeros((n_unique_distances, 2), dtype="<f8")
215 # call cpp function
216 mtk.computeCorrelationFunction(im - im.mean(), output, step, n_threads)
218 if normType is not None:
219 if normType == "variance":
220 # normalise covariance function to a correlation function by dividing by the variance of the image.
221 imVar = im.std() ** 2
222 output[:, 1] = output[:, 1] / imVar
223 elif normType == "first":
224 # normalise covariance function by its first value. Basically it forces it to be 1 at 0.
225 output[:, 1] = output[:, 1] / output[0, 1]
227 d = output[:, 0]
228 c = output[:, 1]
230 return (d, c)
233def betaCovariance(d, lc, b):
234 """
235 Sample correlation function: Beta correlation function
237 $$C(d) = exp(-(d/l_c)^b) $$
239 Parameters
240 ----------
241 d : array
242 The list of distances where the function is evaluated.
244 lc : float
245 The correlation length.
247 b : float
248 The first parameter of the function. If ``b=2`` this is the gaussian function
250 Returns
251 -------
252 c : array
253 The list of values of the function corresponding to the given distances.
255 """
257 c = numpy.exp(-numpy.float_power(d / lc, b))
259 return c
262def gaussianCovariance(d, lc):
263 """
264 Sample correlation function: Gaussian correlation function
267 $$C(d) = exp(-(d/l_c)^2) $$
269 Parameters
270 ----------
271 d : array
272 The list of distances where the function is evaluated.
274 lc : float
275 The correlation length.
277 Returns
278 -------
279 array :
280 The list of values of the function corresponding to the given distances.
282 """
284 return numpy.exp(-numpy.float_power(d / lc, 2))
287def fitCovariance(d, c, functionType="gaussian"):
288 """For autofitting, a helper function...
290 Parameters
291 ----------
292 d : array
293 The list of distances of the covariance.
295 c: array
296 The list of value of covariance corresponding to the distances.
298 functionType : ``{"gaussian","beta"}``
299 The type of covariance function
301 Returns
302 -------
303 popt : array
304 The optimisation parameters
306 """
308 import scipy.optimize
310 d = list(d)
311 c = list(c)
312 if functionType == "gaussian":
313 [popt, pvar] = scipy.optimize.curve_fit(gaussianCovariance, d, c)
314 elif functionType == "beta":
315 [popt, pvar] = scipy.optimize.curve_fit(betaCovariance, d, c)
317 return popt