Coverage for /usr/local/lib/python3.8/site-packages/spam/DIC/grid.py: 96%
85 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1"""
2Library of SPAM functions for defining a regular grid in a reproducible way.
3Copyright (C) 2020 SPAM Contributors
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.
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.
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"""
19import numpy
20import spam.deformation
21import spam.DIC
22import spam.helpers
25def makeGrid(imageSize, nodeSpacing):
26 """
27 Define a grid of correlation points.
29 Parameters
30 ----------
31 imageSize : 3-item list
32 Size of volume to spread the grid inside
34 nodeSpacing : 3-item list or int
35 Spacing between nodes
37 Returns
38 -------
39 nodePositions : nPointsx3 numpy.array
40 Array containing Z, Y, X positions of each point in the grid
41 """
43 if len(imageSize) != 3:
44 print("\tgrid.makeGrid(): imageSize doesn't have three dimensions, exiting")
45 return
47 if type(nodeSpacing) == int or type(nodeSpacing) == float:
48 nodeSpacing = [nodeSpacing] * 3
49 elif len(nodeSpacing) != 3:
50 print(
51 "\tgrid.makeGrid(): nodeSpacing is not an int or float and doesn't have three dimensions, exiting"
52 )
53 return
55 if imageSize[0] == 1:
56 twoD = True
57 else:
58 twoD = False
60 # Note: in this cheap node spacing, the first node is always at a distance of --nodeSpacing-- from the origin
61 # The following could just be done once in principle...
62 nodesMgrid = numpy.mgrid[
63 nodeSpacing[0] : imageSize[0] : nodeSpacing[0],
64 nodeSpacing[1] : imageSize[1] : nodeSpacing[1],
65 nodeSpacing[2] : imageSize[2] : nodeSpacing[2],
66 ]
68 # If twoD then overwrite nodesMgrid
69 if twoD:
70 nodesMgrid = numpy.mgrid[
71 0:1:1,
72 nodeSpacing[1] : imageSize[1] : nodeSpacing[1],
73 nodeSpacing[2] : imageSize[2] : nodeSpacing[2],
74 ]
76 nodesDim = (nodesMgrid.shape[1], nodesMgrid.shape[2], nodesMgrid.shape[3])
78 numberOfNodes = int(nodesMgrid.shape[1] * nodesMgrid.shape[2] * nodesMgrid.shape[3])
80 nodePositions = numpy.zeros((numberOfNodes, 3))
82 nodePositions[:, 0] = nodesMgrid[0].ravel()
83 nodePositions[:, 1] = nodesMgrid[1].ravel()
84 nodePositions[:, 2] = nodesMgrid[2].ravel()
86 return nodePositions, nodesDim
89def getImagettes(
90 im1,
91 nodePosition,
92 halfWindowSize,
93 Phi,
94 im2,
95 searchRange,
96 im1mask=None,
97 im2mask=None,
98 minMaskCoverage=0.0,
99 greyThreshold=[-numpy.inf, numpy.inf],
100 applyF="all",
101 twoD=False,
102):
103 """
104 This function is responsible for extracting correlation windows ("imagettes") from two larger images (im1 and im2).
105 Both spam.correlate.pixelSearch and spam.correlate.register[Multiscale] want a fixed, smaller imagette1
106 and a larger imagette 2 in which to search/interpolate.
108 Parameters
109 ----------
110 im1 : 3D numpy array
111 This is the large input reference image
113 nodePosition : 3-component numpy array of ints
114 This defines the centre of the window to extract from im1.
115 Note: for 2D Z = 0
117 halfWindowSize : 3-component numpy array of ints
118 This defines the half-size of the correlation window,
119 i.e., how many pixels to extract in Z, Y, X either side of the centre.
120 Note: for 2D Z = 0
122 Phi : 4x4 numpy array of floats
123 Phi matrix representing the movement of imagette1,
124 if not equal to `I`, imagette1 is deformed by the non-translation parts of Phi (F)
125 and the displacement is added to the search range (see below)
127 im2 : 3D numpy array
128 This is the large input deformed image
130 searchRange : 6-component numpy array of ints
131 This defines where imagette2 should be extracted with respect to imagette1's position in im1.
132 The 6 components correspond to [ Zbot Ztop Ybot Ytop Xbot Xtop ].
133 If Z, Y and X values are the same, then imagette2 will be displaced and the same size as imagette1.
134 If 'bot' is lower than 'top', imagette2 will be larger in that dimension
136 im1mask : 3D numpy array, optional
137 This needs to be same size as im1, but can be `None` if no mask is wanted.
138 This defines a mask for zones to correlate in im1, 0 means zone not to correlate
139 Default = None
141 im2mask : 3D numpy array, optional
142 This needs to be same size as im2, but can be `None` if no mask is wanted.
143 This defines a mask for zones to correlate in im2, 0 means zone not to correlate
144 Default = None
146 minMaskCoverage : float, optional
147 Threshold for imagette1 non-mask coverage, i.e. how much of imagette1 can be full of mask
148 before it is rejected with returnStatus = -5?
149 Default = 0
151 greyThreshold : two-component list of floats, optional
152 Bottom and top threshold values for mean value of imagette1 to reject it with returnStatus = -5
153 Default = no threshold
155 applyF : string, optional
156 If a non-identity Phi is passed, should the F be applied to the returned imagette1?
157 Options are: 'all', 'rigid', 'no'
158 Default = 'all'
159 Note: as of January 2021, it seems to make more sense to have this as 'all' for pixelSearch, and 'no' for local DIC
161 twoD : bool, optional
162 Are the images two-dimensional?
164 Returns
165 -------
166 Dictionary :
168 'imagette1' : 3D numpy array,
170 'imagette1mask': 3D numpy array of same size as imagette1 or None,
172 'imagette2': 3D numpy array, bigger or equal size to imagette1
174 'imagette2mask': 3D numpy array of same size as imagette2 or None,
176 'returnStatus': int,
177 Describes success in extracting imagette1 and imagette2.
178 If == 1 success, otherwise negative means failure.
180 'pixelSearchOffset': 3-component list of ints
181 Coordinates of the top of the pixelSearch range in im1, i.e., the displacement that needs to be
182 added to the raw pixelSearch output to make it a im1 -> im2 displacement
183 """
184 returnStatus = 1
185 # imagette1mask = None
186 imagette2mask = None
187 intDisplacement = numpy.round(Phi[0:3, 3]).astype(int)
189 assert (
190 len(im1.shape) == len(im2.shape) == 3
191 ), "3D images needed for im1 and im2, if you have 2D images please pad them with im[numpy.newaxis, ...]"
192 if im1mask is not None:
193 assert (
194 len(im1mask.shape) == 3
195 ), "3D image needed for im1mask, if you have 2D images please pad them with im[numpy.newaxis, ...]"
196 if im2mask is not None:
197 assert (
198 len(im2mask.shape) == 3
199 ), "3D image needed for im2mask, if you have 2D images please pad them with im[numpy.newaxis, ...]"
201 # Detect 2D images
202 # if im1.shape[0] == 1:
203 # twoD = True
204 # Impose no funny business in z if in twoD
205 # halfWindowSize[0] = 0
206 # searchRange[0:2] = 0
207 # else:
208 # twoD = False
210 PhiNoDisp = Phi.copy()
211 # PhiNoDisp[0:3,-1] -= intDisplacement
212 PhiNoDisp[0:3, -1] = numpy.zeros(3)
213 if applyF == "rigid":
214 PhiNoDisp = spam.deformation.computeRigidPhi(PhiNoDisp)
216 # If F is not the identity, create a pad to be able to apply F to imagette 1
217 if numpy.allclose(PhiNoDisp, numpy.eye(4)) or applyF == "no":
218 # 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
219 startStopIm1 = [
220 int(nodePosition[0] - halfWindowSize[0]),
221 int(nodePosition[0] + halfWindowSize[0] + 1),
222 int(nodePosition[1] - halfWindowSize[1]),
223 int(nodePosition[1] + halfWindowSize[1] + 1),
224 int(nodePosition[2] - halfWindowSize[2]),
225 int(nodePosition[2] + halfWindowSize[2] + 1),
226 ]
228 # In either case, extract imagette1, now guaranteed to be the right size
229 imagette1def = spam.helpers.slicePadded(im1, startStopIm1)
231 # Check mask
232 if im1mask is None:
233 # no mask1 --> always pas this test (e.g., labelled image)
234 maskVolumeCondition = True
235 imagette1mask = None
236 else:
237 imagette1mask = spam.helpers.slicePadded(im1mask, startStopIm1) != 0
238 maskVolumeCondition = (imagette1mask != 0).mean() >= minMaskCoverage
240 else: # This is the case that we should apply F to imagette1, which requires a pad
241 # 2020-10-06 OS and EA: Add a pad to each dimension of 25% of max(halfWindowSize) to allow space to apply F (no displacement) to imagette1
242 applyPhiPad = int(0.5 * numpy.ceil(max(halfWindowSize)))
244 if twoD:
245 applyPhiPad = (0, applyPhiPad, applyPhiPad)
246 else:
247 applyPhiPad = (applyPhiPad, applyPhiPad, applyPhiPad)
249 # 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
250 startStopIm1 = [
251 int(nodePosition[0] - halfWindowSize[0] - applyPhiPad[0]),
252 int(nodePosition[0] + halfWindowSize[0] + applyPhiPad[0] + 1),
253 int(nodePosition[1] - halfWindowSize[1] - applyPhiPad[1]),
254 int(nodePosition[1] + halfWindowSize[1] + applyPhiPad[1] + 1),
255 int(nodePosition[2] - halfWindowSize[2] - applyPhiPad[2]),
256 int(nodePosition[2] + halfWindowSize[2] + applyPhiPad[2] + 1),
257 ]
259 # In either case, extract imagette1, now guaranteed to be the right size
260 imagette1padded = spam.helpers.slicePadded(im1, startStopIm1)
262 # apply F to imagette 1 padded
263 if twoD:
264 imagette1paddedDef = spam.DIC.applyPhiPython(imagette1padded, PhiNoDisp)
265 else:
266 imagette1paddedDef = spam.DIC.applyPhi(imagette1padded, PhiNoDisp)
268 # undo padding
269 if twoD:
270 imagette1def = imagette1paddedDef[
271 :, applyPhiPad[1] : -applyPhiPad[1], applyPhiPad[2] : -applyPhiPad[2]
272 ]
274 else:
275 imagette1def = imagette1paddedDef[
276 applyPhiPad[0] : -applyPhiPad[0],
277 applyPhiPad[1] : -applyPhiPad[1],
278 applyPhiPad[2] : -applyPhiPad[2],
279 ]
281 # Check mask
282 if im1mask is None:
283 # no mask1 --> always pas this test (e.g., labelled image)
284 maskVolumeCondition = True
285 imagette1mask = None
286 else:
287 imagette1maskPadded = spam.helpers.slicePadded(im1mask, startStopIm1) != 0
289 # apply F to imagette 1 padded
290 # if twoD: imagette1maskPaddedDef = spam.DIC.applyPhiPython(imagette1maskPadded, PhiNoDisp, interpolationOrder=0)
291 # else: imagette1maskPaddedDef = spam.DIC.applyPhiPython(imagette1maskPadded, PhiNoDisp, interpolationOrder=0)
292 imagette1maskPaddedDef = spam.DIC.applyPhiPython(
293 imagette1maskPadded, PhiNoDisp, interpolationOrder=0
294 )
296 # undo padding
297 if twoD:
298 imagette1mask = imagette1maskPaddedDef[
299 :,
300 applyPhiPad[1] : -applyPhiPad[1],
301 applyPhiPad[2] : -applyPhiPad[2],
302 ]
303 else:
304 imagette1mask = imagette1maskPaddedDef[
305 applyPhiPad[0] : -applyPhiPad[0],
306 applyPhiPad[1] : -applyPhiPad[1],
307 applyPhiPad[2] : -applyPhiPad[2],
308 ]
310 maskVolumeCondition = (imagette1mask != 0).mean() >= minMaskCoverage
312 # Make sure imagette is not 0-dimensional in any dimension
313 # Check minMaskVolume
314 if numpy.all(numpy.array(imagette1def.shape) > 0) or (
315 twoD and numpy.all(numpy.array(imagette1def.shape[1:3]) > 0)
316 ):
317 # ------------ Grey threshold low --------------- and -------------- Grey threshold high -----------
318 if (
319 numpy.nanmean(imagette1def) > greyThreshold[0]
320 and numpy.nanmean(imagette1def) < greyThreshold[1]
321 ):
322 if maskVolumeCondition:
323 # Slice for image 2
324 # 2020-09-25 OS and EA: Prepare startStop array for imagette 2 to be extracted with new slicePadded
325 # Extract it...
326 startStopIm2 = [
327 int(
328 nodePosition[0]
329 - halfWindowSize[0]
330 + intDisplacement[0]
331 + searchRange[0]
332 ),
333 int(
334 nodePosition[0]
335 + halfWindowSize[0]
336 + intDisplacement[0]
337 + searchRange[1]
338 + 1
339 ),
340 int(
341 nodePosition[1]
342 - halfWindowSize[1]
343 + intDisplacement[1]
344 + searchRange[2]
345 ),
346 int(
347 nodePosition[1]
348 + halfWindowSize[1]
349 + intDisplacement[1]
350 + searchRange[3]
351 + 1
352 ),
353 int(
354 nodePosition[2]
355 - halfWindowSize[2]
356 + intDisplacement[2]
357 + searchRange[4]
358 ),
359 int(
360 nodePosition[2]
361 + halfWindowSize[2]
362 + intDisplacement[2]
363 + searchRange[5]
364 + 1
365 ),
366 ]
367 imagette2 = spam.helpers.slicePadded(im2, startStopIm2)
369 if im2mask is not None:
370 imagette2mask = spam.helpers.slicePadded(im2mask, startStopIm2)
372 # Failed minMaskVolume condition
373 else:
374 returnStatus = -5
375 imagette1def = None
376 imagette2 = None
378 # Failed greylevel condition
379 else:
380 returnStatus = -5
381 imagette1def = None
382 imagette2 = None
384 # Failed 0-dimensional imagette test
385 else:
386 returnStatus = -5
387 imagette1def = None
388 imagette2 = None
390 return {
391 "imagette1": imagette1def,
392 "imagette1mask": imagette1mask,
393 "imagette2": imagette2,
394 "imagette2mask": imagette2mask,
395 "returnStatus": returnStatus,
396 "pixelSearchOffset": searchRange[0::2] + intDisplacement,
397 }