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 677cf728 authored by Edward Andò's avatar Edward Andò
Browse files

spheroid generator added to kalisphera (not exactly the right place, but OK)

parent 6d68f32c
Pipeline #61279 passed with stages
in 25 minutes and 47 seconds
__all__ = [ "kalisphera" ]
__all__ = [ "kalisphera", "spheroid"]
from .kalisphera import *
from .spheroid import *
......@@ -98,17 +98,17 @@ def makeBlurryNoisySphere(dims, centre, radius, blur=0, noise=0, flatten=True, b
Dimensions of volume to create
centre : 1D or 2D numpy array
Either a 3-component vector, or an Nx3 matrix of sphere centres to draw with respect to 0,0,0 in `vol`
Either a 3-component vector, or an Nx3 array of sphere centres to draw with respect to 0,0,0 in `vol`
radius : float or 1D numpy array
Raduis(ii) of spheres to draw in `vol`
Radius(ii) of spheres to draw in `vol`
blur : float, optional
Standard deviation of the blur kernel (in ~pixels)
Default = 0
noise : float, optional
Standard devitaion of the noise (in greylevels), reminder background=0 and sphere=1
Standard devitaion of the noise (in greylevels), noise is on the scale: background=0 and sphere=1
Default = 0
flatten : bool, optional
......
"""
Library of SPAM functions for generating 3D spheroids
Copyright (C) 2020 SPAM Contributors
This program is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option)
any later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
"""
import numpy
def makeBlurryNoisySpheroid(dims, centres, axisLengthsAC, orientations, blur=0, noise=0, flatten=True, background=0.25, foreground=0.75):
"""
This function creates a sheroid or series of spheroids in a 3D volume, adding noise and blur if requested.
Note that this does not handle the partial volume effect like kalisphera, and so some blur or downscaling is recommended!
The base code is from https://sbrisard.github.io/posts/20150930-orientation_correlations_among_rice_grains-06.html
Parameters
----------
dims : 3-component list
Dimensions of volume to create
centre : 1D or 2D numpy array
Either a 3-component vector, or an Nx3 array of centres to draw with respect to 0,0,0 in `vol`
axisLengthsAC : 1D or 2D numpy array of floats
Half-axis lengths for each spheroid. If c>a, a prolate (rice-like) is generated; while a>c yields an oblate (lentil-like).
This is either a 2-component array for 1 particle, or a Nx2 array of values.
orientations : 1D or 2D numpy array of floats
Orientation vectors of the axis of rotational symmetry for the spheroid.
This is either a 3-component ZYX array for 1 particle, or a Nx3 array of values.
blur : float, optional
Standard deviation of the blur kernel (in ~pixels)
Default = 0
noise : float, optional
Standard devitaion of the noise (in greylevels), reminder background=0 and sphere=1
Default = 0
flatten : bool, optional
Flatten greyvalues >1 to 1 (caused by overlaps)?
Default = True
background : float, optional
Desired mean greyvalue of the background
Default = 0.25
foreground : float, optional
Desired mean greyvalue of the background
Default = 0.75
Returns
-------
vol : the input array
"""
# 3D
D = 3
vol = numpy.zeros(dims, dtype='<f4')
centres = numpy.array(centres)
axisLengthsAC = numpy.array(axisLengthsAC)
orientations = numpy.array(orientations)
dims = numpy.array(dims)
if len(centres.shape) == 1:
nParticles = 1
# There is just one spheroid
centres = numpy.array([centres])
assert len(axisLengthsAC.shape) == 1
assert axisLengthsAC.shape[0] == 2
axisLengthsAC = numpy.array([axisLengthsAC])
assert len(orientations.shape) == 1
assert orientations.shape[0] == 3
assert numpy.isclose(numpy.linalg.norm(orientations), 1)
orientations = numpy.array([orientations])
elif len(centres.shape) == 2:
nParticles = centres.shape[0]
# more than one spheroid
assert len(axisLengthsAC.shape) == 2
assert axisLengthsAC.shape[0] == nParticles
assert axisLengthsAC.shape[1] == 2
assert len(orientations.shape) == 2
assert orientations.shape[0] == nParticles
assert orientations.shape[1] == 3
def spheroidDistanceFunction(x, invQ):
"""
Ordering of points: ``x[i, ...]`` is the i-th coordinate
"""
y = numpy.tensordot(invQ, x, axes=([-1], [0]))
numpy.multiply(x, y, y)
return numpy.sum(y, axis=0)
# pool?
for centre, axisLengthAC, orientation in zip(centres, axisLengthsAC, orientations):
# Preamble
p = numpy.outer(orientation, orientation)
q = numpy.eye(D, dtype=numpy.float64) - p
Q = axisLengthAC[1]**2*p+axisLengthAC[0]**2*q
invQ = p/axisLengthAC[1]**2+q/axisLengthAC[0]**2
# Must implement local bounding box calculation here...
bb = numpy.sqrt(numpy.diag(Q))
h = 1
#bb = dims
i_max = numpy.ceil(bb/h-0.5)
bb = i_max*h
shape = 2*i_max+1
slices = [slice(-x, x, i*1j) for (x, i) in zip(bb, shape)]
x = numpy.mgrid[slices]
# Thresholding distance function!!
spheroidLocal = spheroidDistanceFunction(x, invQ) < 1
spheroidLocalShape = numpy.array(spheroidLocal.shape).astype(int)
spheroidLocalHalfWindow = (spheroidLocalShape - 1) / 2
spheroidGlobalTop = (centre-spheroidLocalHalfWindow).astype(int)
spheroidGlobalBot = (centre+spheroidLocalHalfWindow+1).astype(int)
# How much the local BB is going over the edges of the image
spheroidOffsetTop = numpy.zeros(3, dtype=int)
spheroidOffsetBot = numpy.zeros(3, dtype=int)
if spheroidGlobalTop[0] < 0: spheroidOffsetTop[0] = -spheroidGlobalTop[0]
if spheroidGlobalTop[1] < 0: spheroidOffsetTop[1] = -spheroidGlobalTop[1]
if spheroidGlobalTop[2] < 0: spheroidOffsetTop[2] = -spheroidGlobalTop[2]
if spheroidGlobalBot[0] > dims[0]-1: spheroidOffsetBot[0] = spheroidGlobalBot[0] - dims[0]
if spheroidGlobalBot[1] > dims[1]-1: spheroidOffsetBot[1] = spheroidGlobalBot[1] - dims[1]
if spheroidGlobalBot[2] > dims[2]-1: spheroidOffsetBot[2] = spheroidGlobalBot[2] - dims[2]
sliceGlobal = (slice(spheroidGlobalTop[0]+spheroidOffsetTop[0],spheroidGlobalBot[0]-spheroidOffsetBot[0]),
slice(spheroidGlobalTop[1]+spheroidOffsetTop[1],spheroidGlobalBot[1]-spheroidOffsetBot[1]),
slice(spheroidGlobalTop[2]+spheroidOffsetTop[2],spheroidGlobalBot[2]-spheroidOffsetBot[2]))
sliceLocal = (slice(spheroidOffsetTop[0],spheroidLocalShape[0]-spheroidOffsetBot[0]),
slice(spheroidOffsetTop[1],spheroidLocalShape[1]-spheroidOffsetBot[1]),
slice(spheroidOffsetTop[2],spheroidLocalShape[2]-spheroidOffsetBot[2]))
vol[sliceGlobal] += spheroidLocal[sliceLocal]
if flatten: vol[vol>1]=1
if blur != 0:
import scipy.ndimage
vol = scipy.ndimage.filters.gaussian_filter(vol, sigma=blur)
if noise != 0:
vol += numpy.random.normal(size=dims, scale=noise)
vol *= float(foreground)-float(background)
vol += float(background)
return vol
......@@ -386,7 +386,7 @@ def getLabel(labelledVolume, label, boundingBoxes=None, centresOfMass=None, marg
"""
Helper function to extract labels from a labelled image/volume.
A dictionary is returned with the a subvolume around the particle.
Passing boundingBoxes and centresOfMass is highly recommended.
Passing boundingBoxes and centresOfMass is highly recommended.
Parameters
----------
......@@ -1504,7 +1504,7 @@ class Spheroid:
#imLab = imLab + labelDummy
#Update the labels
#imLab = spam.label.makeLabelsSequential(imLab)
#pbar.finish()
#print("\tlabel.fixUndersegmentation(): From "+str(len(listLabels))+" undersegmented labels, we successfully worked on "+str(successCounter))
#return imLab
......@@ -1544,7 +1544,7 @@ def convexVolume(lab, boundingBoxes=None, centresOfMass=None, volumes=None, numb
--------
convexVolume : lab.max()x1 array of floats with the convex volume.
Note
----
convexVolume can only be computed for particles with volume greater than 3 voxels. If it is not the case, it will return 0.
......@@ -1588,22 +1588,22 @@ def convexVolume(lab, boundingBoxes=None, centresOfMass=None, volumes=None, numb
except:
qResults.put( [job, 0])
numberOfJobs = numpy.max(lab)
qJobs = multiprocessing.Queue()
qResults = multiprocessing.Queue()
# print "Master: Adding jobs to queues"
for x in range(1,numberOfJobs+1):
qJobs.put(x)
for i in range(numberOfThreads):
qJobs.put("STOP")
# print "Master: Launching workers"
for i in range(numberOfThreads):
p = multiprocessing.Process(target=worker, args=(i, qJobs, qResults, ))
p.start()
if verbose:
pbar = progressbar.ProgressBar(maxval=numberOfJobs).start()
finishedThreads = 0
......@@ -1635,15 +1635,15 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
-----------
lab : 3D numpy array
Labelled image
PhiField : (multidimensional x 4 x 4 numpy array of floats)
Spatial field of Phis
returnStatus : lab.max()x1 array of ints, optional
Array with the return status for each label (usually returned by ``spam-ddic``)
If not defined (Default = None), all the labels will be moved
If returnStatus[i] == 2, the label will be moved, otherwise is omitted and erased from the final image
boundingBoxes : lab.max()x6 array of ints, optional
Bounding boxes in format returned by ``boundingBoxes``.
If not defined (Default = None), it is recomputed by running ``boundingBoxes``
......@@ -1651,28 +1651,28 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
centresOfMass : lab.max()x3 array of floats, optional
Centres of mass in format returned by ``centresOfMass``.
If not defined (Default = None), it is recomputed by running ``centresOfMass``
margin : int, optional
Margin, in pixels, to take in each label.
Default = 3
PhiCOM : bool, optional
Apply Phi to centre of mass of particle?, otherwise it will be applied in the middle of the particle\'s bounding box.
Default = True
threshold : float, optional
Threshold to keep interpolated voxels in the binary image.
Default = 0.5
labelDilate : int, optional
Number of times label should be dilated/eroded before returning it.
If ``labelDilate > 0`` a dilated label is returned, while ``labelDilate < 0`` returns an eroded label.
Default = 0
numberOfThreads : int, optional
Number of Threads for multiprocessing
Default = 1
Returns
--------
labOut : 3D numpy array
......@@ -1682,7 +1682,7 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
# Define worker for multiprocessing
def worker(worker_number, q_jobs, q_results):
while True:
job = q_jobs.get()
......@@ -1700,7 +1700,7 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
centresOfMass=centresOfMass,
extractCube=True)
# Check that the label exist
if getLabelReturn is not None:
if getLabelReturn is not None:
# Get Phi field
Phi = PhiField[label].copy()
# Phi will be split into a local part and a part of floored displacements
......@@ -1730,7 +1730,7 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
topOfSlice = numpy.array([getLabelReturn['boundingBoxCube'][0] + disp[0],
getLabelReturn['boundingBoxCube'][2] + disp[1],
getLabelReturn['boundingBoxCube'][4] + disp[2]])
botOfSlice = numpy.array([getLabelReturn['boundingBoxCube'][1] + disp[0],
getLabelReturn['boundingBoxCube'][3] + disp[1],
getLabelReturn['boundingBoxCube'][5] + disp[2]])
......@@ -1745,7 +1745,7 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
grainSlice = (slice(topOfSliceCrop[0], botOfSliceCrop[0]),
slice(topOfSliceCrop[1], botOfSliceCrop[1]),
slice(topOfSliceCrop[2], botOfSliceCrop[2]))
# Update labSubvolDefMask
labSubvolDefMaskCrop = labSubvolDefMask[topOfSliceCrop[0]-topOfSlice[0] : labSubvolDefMask.shape[0]-1 + (botOfSliceCrop[0]-botOfSlice[0]),
topOfSliceCrop[1]-topOfSlice[1] : labSubvolDefMask.shape[1]-1 + (botOfSliceCrop[1]-botOfSlice[1]),
......@@ -1770,7 +1770,7 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
labOut = numpy.zeros_like(lab, dtype=spam.label.labelType)
# Get number of labels
numberOfLabels = min(lab.max(), PhiField.shape[0]-1)
# Setting up queues
q_jobs = multiprocessing.Queue()
q_results = multiprocessing.Queue()
......@@ -1783,7 +1783,7 @@ def moveLabels(lab, PhiField, returnStatus = None, boundingBoxes=None, centresOf
else: # Add the particles
q_jobs.put(label)
numberOfLabelsToMove += 1
for i in range(numberOfThreads): q_jobs.put("STOP")
# Launching workers
for i in range(numberOfThreads):
......@@ -1819,10 +1819,10 @@ def erodeLabels(lab, erosion=1, boundingBoxes=None, centresOfMass=None, numberOf
-----------
lab : 3D numpy array
Labelled image
erosion : int, optional
Erosion level
boundingBoxes : lab.max()x6 array of ints, optional
Bounding boxes in format returned by ``boundingBoxes``.
If not defined (Default = None), it is recomputed by running ``boundingBoxes``
......@@ -1830,16 +1830,16 @@ def erodeLabels(lab, erosion=1, boundingBoxes=None, centresOfMass=None, numberOf
centresOfMass : lab.max()x3 array of floats, optional
Centres of mass in format returned by ``centresOfMass``.
If not defined (Default = None), it is recomputed by running ``centresOfMass``
numberOfThreads : int, optional
Number of Threads for multiprocessing
Default = 1
Returns
--------
erodeImage : 3D numpy array
New labelled image with the eroded labels.
Note
----
The function makes use of spam.label.moveLabels() to generate the eroded image.
......@@ -1872,7 +1872,7 @@ def convexFillHoles(lab, boundingBoxes=None, centresOfMass=None):
-----------
lab : 3D numpy array
Labelled image
boundingBoxes : lab.max()x6 array of ints, optional
Bounding boxes in format returned by ``boundingBoxes``.
If not defined (Default = None), it is recomputed by running ``boundingBoxes``
......@@ -1885,7 +1885,7 @@ def convexFillHoles(lab, boundingBoxes=None, centresOfMass=None):
--------
labOut : 3D numpy array
New labelled image.
Note
----
The function works nicely for convex particles. For non-convex particles, it will alter the shape.
......@@ -1947,11 +1947,11 @@ def convexFillHoles(lab, boundingBoxes=None, centresOfMass=None):
# Update the progressbar
widgets[0] = progressbar.FormatLabel("{}/{} ".format(i, numberOfLabels))
pbar.update(i)
return labOut
def getNeighbours(lab, listOfLabels, method = 'getLabel', parameter = None, centresOfMass = None, boundingBoxes = None):
"""
This function computes the neighbours for a list of labels.
......@@ -1959,21 +1959,21 @@ def getNeighbours(lab, listOfLabels, method = 'getLabel', parameter = None, cent
-----------
lab : 3D numpy array
Labelled image
listOfLabels : list of ints
List of labels to which the neighbours will be computed
method : string
Method to compute the neighbours.
'getLabel' : The neighbours are the labels inside the subset obtained through spam.getLabel()
'mesh' : The neighbours are computed using a tetrahedral connectivity matrix
Default = 'getLabel'
parameter : int
Parameter controlling each method.
For 'getLabel', it correspond to the size of the subset. Default = 3
For 'mesh', it correspond to the size of the alpha shape used for carving the mesh. Default = 5*meanDiameter.
boundingBoxes : lab.max()x6 array of ints, optional
Bounding boxes in format returned by ``boundingBoxes``.
If not defined (Default = None), it is recomputed by running ``boundingBoxes``
......@@ -2020,7 +2020,7 @@ def getNeighbours(lab, listOfLabels, method = 'getLabel', parameter = None, cent
neighboursLabel = neighboursLabel[~numpy.in1d(neighboursLabel, 0)]
# Add the neighbours to the list
neighbours.append(neighboursLabel)
elif method == 'mesh':
# Compute parameter if needed
if parameter == None:
......@@ -2036,5 +2036,5 @@ def getNeighbours(lab, listOfLabels, method = 'getLabel', parameter = None, cent
neighbours.append(neighboursLabel)
else:
print('spam.label.getNeighbours(): Wrong method, aborting')
return neighbours
......@@ -79,6 +79,50 @@ class TestFunctionLabel(unittest.TestCase):
self.assertAlmostEqual(numpy.sum(Box)/theoreticalVolume, 1.0, places=2)
def test_spheroid(self):
DIMS = (50,50,50)
# start with a noisy blurry sphere and check volume
sphereRadius = 20
sphereTheoreticalVolume = (4/3)*numpy.pi*sphereRadius**3
sphere = spam.kalisphera.makeBlurryNoisySpheroid(DIMS, numpy.array(DIMS)/2, [sphereRadius,sphereRadius], [1,0,0], blur=0.5, noise=0.01, background=0.25, foreground=0.75)
sphereVolume = numpy.sum(sphere > 0.5)
self.assertTrue(numpy.isclose(sphereVolume/sphereTheoreticalVolume, 1, atol=0.2))
# Crash outside the edges, centre at 0
sphereEighth = spam.kalisphera.makeBlurryNoisySpheroid(DIMS, [0,0,0], [sphereRadius,sphereRadius], [1,0,0], blur=0, noise=0, background=0.0, foreground=1.0)
self.assertTrue(numpy.isclose(sphereEighth.sum()/sphereTheoreticalVolume, 0.125, atol=0.05))
# Crash outside the edges, centre at opposite edge
sphereEighth = spam.kalisphera.makeBlurryNoisySpheroid(DIMS, DIMS, [sphereRadius,sphereRadius], [1,0,0], blur=0, noise=0, background=0.0, foreground=1.0)
self.assertTrue(numpy.isclose(sphereEighth.sum()/sphereTheoreticalVolume, 0.125, atol=0.05))
# Let's make two objects -- use the two corner cases
sphereTwoEighths = spam.kalisphera.makeBlurryNoisySpheroid(DIMS,
[[0,0,0], DIMS],
[[sphereRadius,sphereRadius],[sphereRadius,sphereRadius]],
[[1,0,0], [1,0,0]],
blur=0,
noise=0,
background=0.0,
foreground=1.0)
self.assertTrue(numpy.isclose(sphereTwoEighths.sum()/sphereTheoreticalVolume, 0.25, atol=0.05))
# Let's actually make a spheroid... Lentils first (always) -- means A>C
lentil = spam.kalisphera.makeBlurryNoisySpheroid(DIMS,
numpy.array(DIMS)/2,
[20,10],
[1,0,0], blur=0, noise=0, background=0.0, foreground=1.0)
lentilTheoreticalVolume = (4/3)*numpy.pi*(20*20*10)
self.assertTrue(numpy.isclose(lentil.sum()/lentilTheoreticalVolume, 1, atol=0.1))
# check orientation, max eigenvector for a lentil
orientation = spam.label.momentOfInertia(lentil > 0.5)[1][1, 0:3]
self.assertTrue(numpy.allclose(numpy.abs(orientation), [1,0,0], atol=0.1))
ellipseAxes = spam.label.ellipseAxes(lentil > 0.5)[1, 0:3]
self.assertTrue(numpy.allclose(ellipseAxes , [20,20,10], atol=0.1))
if __name__ == '__main__':
unittest.main()
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