Coverage for /usr/local/lib/python3.8/site-packages/spam/kalisphera/spheroid.py: 100%

1""" 

2Library of SPAM functions for generating 3D spheroids 

3Copyright (C) 2020 SPAM Contributors 

4 

5This program is free software: you can redistribute it and/or modify it 

6under the terms of the GNU General Public License as published by the Free 

7Software Foundation, either version 3 of the License, or (at your option) 

8any later version. 

9 

10This program is distributed in the hope that it will be useful, but WITHOUT 

11ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

13more details. 

14 

15You should have received a copy of the GNU General Public License along with 

16this program. If not, see <http://www.gnu.org/licenses/>. 

17""" 

18import numpy 

19 

20def makeBlurryNoisySpheroid(dims, centres, axisLengthsAC, orientations, blur=0, noise=0, flatten=True, background=0.25, foreground=0.75): 

21 """ 

22 This function creates a sheroid or series of spheroids in a 3D volume, adding noise and blur if requested. 

23 Note that this does not handle the partial volume effect like kalisphera, and so some blur or downscaling is recommended! 

24 

25 The base code is from https://sbrisard.github.io/posts/20150930-orientation_correlations_among_rice_grains-06.html 

26 

27 Parameters 

28 ---------- 

29 dims : 3-component list 

30 Dimensions of volume to create 

31 

32 centre : 1D or 2D numpy array 

33 Either a 3-component vector, or an Nx3 array of centres to draw with respect to 0,0,0 in `vol` 

34 

35 axisLengthsAC : 1D or 2D numpy array of floats 

36 Half-axis lengths for each spheroid. If c>a, a prolate (rice-like) is generated; while a>c yields an oblate (lentil-like). 

37 This is either a 2-component array for 1 particle, or a Nx2 array of values. 

38 

39 orientations : 1D or 2D numpy array of floats 

40 Orientation vectors of the axis of rotational symmetry for the spheroid. 

41 This is either a 3-component ZYX array for 1 particle, or a Nx3 array of values. 

42 

43 blur : float, optional 

44 Standard deviation of the blur kernel (in ~pixels) 

45 Default = 0 

46 

47 noise : float, optional 

48 Standard devitaion of the noise (in greylevels), reminder background=0 and sphere=1 

49 Default = 0 

50 

51 flatten : bool, optional 

52 Flatten greyvalues >1 to 1 (caused by overlaps)? 

53 Default = True 

54 

55 background : float, optional 

56 Desired mean greyvalue of the background 

57 Default = 0.25 

58 

59 foreground : float, optional 

60 Desired mean greyvalue of the background 

61 Default = 0.75 

62 

63 Returns 

64 ------- 

65 vol : the input array 

66 """ 

67 # 3D 

68 D = 3 

69 

70 vol = numpy.zeros(dims, dtype='<f4') 

71 

72 centres = numpy.array(centres) 

73 axisLengthsAC = numpy.array(axisLengthsAC) 

74 orientations = numpy.array(orientations) 

75 dims = numpy.array(dims) 

76 

77 if len(centres.shape) == 1: 

78 nParticles = 1 

79 

80 # There is just one spheroid 

81 centres = numpy.array([centres]) 

82 

83 assert len(axisLengthsAC.shape) == 1 

84 assert axisLengthsAC.shape[0] == 2 

85 axisLengthsAC = numpy.array([axisLengthsAC]) 

86 

87 assert len(orientations.shape) == 1 

88 assert orientations.shape[0] == 3 

89 assert numpy.isclose(numpy.linalg.norm(orientations), 1) 

90 orientations = numpy.array([orientations]) 

91 

92 elif len(centres.shape) == 2: 

93 nParticles = centres.shape[0] 

94 

95 # more than one spheroid 

96 assert len(axisLengthsAC.shape) == 2 

97 assert axisLengthsAC.shape[0] == nParticles 

98 assert axisLengthsAC.shape[1] == 2 

99 

100 assert len(orientations.shape) == 2 

101 assert orientations.shape[0] == nParticles 

102 assert orientations.shape[1] == 3 

103 

104 def spheroidDistanceFunction(x, invQ): 

105 """ 

106 Ordering of points: ``x[i, ...]`` is the i-th coordinate 

107 """ 

108 y = numpy.tensordot(invQ, x, axes=([-1], [0])) 

109 numpy.multiply(x, y, y) 

110 return numpy.sum(y, axis=0) 

111 

112 # pool? 

113 for centre, axisLengthAC, orientation in zip(centres, axisLengthsAC, orientations): 

114 # Preamble 

115 p = numpy.outer(orientation, orientation) 

116 q = numpy.eye(D, dtype=numpy.float64) - p 

117 Q = axisLengthAC[1]**2*p+axisLengthAC[0]**2*q 

118 invQ = p/axisLengthAC[1]**2+q/axisLengthAC[0]**2 

119 

120 # Must implement local bounding box calculation here... 

121 bb = numpy.sqrt(numpy.diag(Q)) 

122 

123 

124 h = 1 

125 #bb = dims 

126 i_max = numpy.ceil(bb/h-0.5) 

127 bb = i_max*h 

128 shape = 2*i_max+1 

129 

130 slices = [slice(-x, x, i*1j) for (x, i) in zip(bb, shape)] 

131 x = numpy.mgrid[slices] 

132 

133 # Thresholding distance function!! 

134 spheroidLocal = spheroidDistanceFunction(x, invQ) < 1 

135 

136 spheroidLocalShape = numpy.array(spheroidLocal.shape).astype(int) 

137 spheroidLocalHalfWindow = (spheroidLocalShape - 1) / 2 

138 #GP 2021-06-23: Round to the closest int 

139 centre = numpy.floor(centre) 

140 spheroidGlobalTop = (centre-spheroidLocalHalfWindow).astype(int) 

141 spheroidGlobalBot = (centre+spheroidLocalHalfWindow+1).astype(int) 

142 

143 # How much the local BB is going over the edges of the image 

144 spheroidOffsetTop = numpy.zeros(3, dtype=int) 

145 spheroidOffsetBot = numpy.zeros(3, dtype=int) 

146 if spheroidGlobalTop[0] < 0: spheroidOffsetTop[0] = -spheroidGlobalTop[0] 

147 if spheroidGlobalTop[1] < 0: spheroidOffsetTop[1] = -spheroidGlobalTop[1] 

148 if spheroidGlobalTop[2] < 0: spheroidOffsetTop[2] = -spheroidGlobalTop[2] 

149 

150 if spheroidGlobalBot[0] > dims[0]-1: spheroidOffsetBot[0] = spheroidGlobalBot[0] - dims[0] 

151 if spheroidGlobalBot[1] > dims[1]-1: spheroidOffsetBot[1] = spheroidGlobalBot[1] - dims[1] 

152 if spheroidGlobalBot[2] > dims[2]-1: spheroidOffsetBot[2] = spheroidGlobalBot[2] - dims[2] 

153 

154 sliceGlobal = (slice(spheroidGlobalTop[0]+spheroidOffsetTop[0],spheroidGlobalBot[0]-spheroidOffsetBot[0]), 

155 slice(spheroidGlobalTop[1]+spheroidOffsetTop[1],spheroidGlobalBot[1]-spheroidOffsetBot[1]), 

156 slice(spheroidGlobalTop[2]+spheroidOffsetTop[2],spheroidGlobalBot[2]-spheroidOffsetBot[2])) 

157 sliceLocal = (slice(spheroidOffsetTop[0],spheroidLocalShape[0]-spheroidOffsetBot[0]), 

158 slice(spheroidOffsetTop[1],spheroidLocalShape[1]-spheroidOffsetBot[1]), 

159 slice(spheroidOffsetTop[2],spheroidLocalShape[2]-spheroidOffsetBot[2])) 

160 

161 vol[sliceGlobal] += spheroidLocal[sliceLocal] 

162 

163 if flatten: vol[vol>1]=1 

164 

165 if blur != 0: 

166 import scipy.ndimage 

167 vol = scipy.ndimage.gaussian_filter(vol, sigma=blur) 

168 

169 if noise != 0: 

170 vol += numpy.random.normal(size=dims, scale=noise) 

171 

172 vol *= float(foreground)-float(background) 

173 vol += float(background) 

174 

175 return vol 

176