Coverage for /usr/local/lib/python3.10/site-packages/spam/DIC/pixelSearch.py: 91%
97 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-24 17:26 +0000
1"""
2Library of SPAM image correlation functions.
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 multiprocessing
21import numpy
22import progressbar
23import spam.label
24from spambind.DIC.DICToolkit import pixelSearch as pixelSearchCPP
27def _errorCalc(im1, im2):
28 return numpy.nansum(numpy.square(numpy.subtract(im1, im2))) / numpy.nansum(im1)
31def _pixelSearch(imagette1, imagette2, imagette1mask=None, imagette2mask=None, returnError=False):
32 """
33 This function performs a pixel-by-pixel search in 3D for a small reference imagette1 within a larger imagette2.
35 The normalised correlation coefficient (NCC) is computed for EVERY combination of the displacements of imagette1 defined within imagette2.
36 What is returned is the highest NCC value obtained and the position where it was obtained with respect to the origin in imagette2.
38 Values of NCC > 0.99 can generally be trusted.
40 Parameters
41 ----------
42 imagette1 : 3D numpy array of floats
43 Imagette 1 is the smaller reference image
45 imagette2 : 3D numpy array of floats
46 Imagette 2 is the bigger image inside which to search
48 imagette1mask : 3D numpy array of bools, optional
49 A mask for imagette1 to define which pixels to include in the correlation
50 (True = Correlate these pixels, False = Skip),
51 Default = no mask
53 imagette2mask : 3D numpy array of bools, optional
54 A mask for imagette2 to define which pixels to include in the correlation
55 (True = Correlate these pixels, False = Skip),
56 Default = no mask
58 returnError : bool, optional
59 Return a normalised sum of squared differences error compatible with
60 register() family of functions? If yes, it's the last returned variable
61 Default = False
63 Returns
64 -------
65 p : 3-component vector
66 Z, Y, X position with respect to origin in imagette2 of imagette1 to get best NCC
68 cc : float
69 Normalised correlation coefficient, ~0.5 is random, >0.99 is very good correlation.
70 """
71 # Note
72 # ----
73 # It important to remember that the C code runs A BIT faster in its current incarnation when it has
74 # a cut-out im2 to deal with (this is related to processor optimistaions).
75 # Cutting out imagette2 to just fit around the search range might save a bit of time
76 assert numpy.all(imagette2.shape >= imagette1.shape), "spam.DIC.pixelSearch(): imagette2 should be bigger or equal to imagette1 in all dimensions"
78 if imagette1mask is not None:
79 assert imagette1.shape == imagette1mask.shape, "spam.DIC.pixelSearch: imagette1mask ({}) should have same size as imagette1 ({})".format(imagette1mask.shape, imagette1.shape)
80 imagette1 = imagette1.astype("<f4")
81 imagette1[imagette1mask == 0] = numpy.nan
83 if imagette2mask is not None:
84 assert imagette2.shape == imagette2mask.shape, "spam.DIC.pixelSearch: imagette2mask ({}) should have same size as imagette2 ({})".format(imagette2mask.shape, imagette2.shape)
85 imagette2 = imagette2.astype("<f4")
86 imagette2[imagette2mask == 0] = numpy.nan
88 # Run the actual pixel search
89 returns = numpy.zeros(4, dtype="<f4")
90 pixelSearchCPP(imagette1.astype("<f4"), imagette2.astype("<f4"), returns)
92 if returnError:
93 error = _errorCalc(
94 imagette1,
95 imagette2[
96 int(returns[0]) : int(returns[0]) + imagette1.shape[0],
97 int(returns[1]) : int(returns[1]) + imagette1.shape[1],
98 int(returns[2]) : int(returns[2]) + imagette1.shape[2],
99 ],
100 )
101 return numpy.array(returns[0:3]), returns[3], error
103 else:
104 return numpy.array(returns[0:3]), returns[3]
107def pixelSearchDiscrete(
108 lab1, im1, im2, searchRange, PhiField=None, labelsToCorrelate=None, boundingBoxes=None, centresOfMass=None, applyF="all", labelDilate=1, volThreshold=100, numProc=multiprocessing.cpu_count()
109):
110 """
111 Discrete pixel search over all labels in a labelled volume, given a displacement search range.
112 This is the python function called by spam-pixelSearch when given a labelled image.
114 Parameters
115 ----------
116 lab1 : 2D or 3D numpy array of ints
117 Array representing L labelled objects to search for in the reference configuration max(lab1) = L
119 im1 : 2D or 3D numpy array of greylevels
120 Array representing greylevels in the reference configuration, same size as lab1
122 im2 : 2D or 3D numpy array of greylevels
123 Array representing greylevels in the deformed configuration, same size as lab1
125 searchRange : 1D numpy array of signed ints
126 Array defining search range in [low Z, high Z, low Y, high Y, low X, high X]
128 PhiField : L+1 x 4 x 4 numpy array of floats, optional
129 Optional initial guess for Phi for each particle, default guess for each particle = identity matrix
131 labelsToCorrelate : 1D numpy array of ints, optional
132 Array containing which labels to correlate, default=all of them
134 centresOfMass : L+1 x 3 numpy array of floats, optional
135 Label centres of mass for particles as coming out of spam.label.boundingBoxes, default = it's recalculated on lab1
137 boundingBoxes : L+1 x 3 numpy array of ints, optional
138 Bounding boxes for particles as coming out of spam.label.boundingBoxes, default = it's recalculated on lab1
140 applyF : string, optional
141 Apply the F part of Phi guess? Accepted values are:\n\t"all": apply all of F' + '\n\t"rigid": apply rigid part (mostly rotation) \n\t"no": don\'t apply it "all" is default
143 labelDilate : int, optional
144 Number of times to dilate labels. Default = 1
146 volThreshold : int, optional
147 Volume threshold below which labels are ignored. Default = 100
149 numProc : int, optional
150 Number of processes to use for the calculation, default = multiprocessing.cpu_count()
152 Returns
153 -------
154 Dictionary including the following keys for all labels, but only filled in correlated labels:
155 - PhiField : L+1 x 4 x 4 numpy array of floats: Phi for each label
156 - pixelSearchCC: L+1 numpy array of floats: Correlation Coefficient for each label bounded [0, 1]
157 - error: L+1 numpy array of floats: SSQD per pixel for each label
158 - returnStatus: L+1 numpy array of ints: returnStatus for each label
159 - deltaPhiNorm: L+1 numpy array of floats: 1 for each label TODO decide if this is worth doing here
160 - iterations: L+1 numpy array of ints: 1 for each label TODO decide if this is worth doing here
161 """
163 nLabels = numpy.max(lab1)
165 if PhiField is None:
166 PhiField = numpy.zeros((nLabels + 1, 4, 4), dtype="<f4")
167 for label in range(nLabels + 1):
168 PhiField[label] = numpy.eye(4)
169 else:
170 assert PhiField.shape[0] == nLabels + 1
171 assert PhiField.shape[1] == 4
172 assert PhiField.shape[2] == 4
174 if labelsToCorrelate is None:
175 labelsToCorrelate = numpy.arange(1, nLabels + 1)
176 else:
177 assert labelsToCorrelate.dtype.kind in numpy.typecodes["AllInteger"]
178 assert numpy.min(labelsToCorrelate) >= 1
179 assert numpy.max(labelsToCorrelate) <= nLabels
181 if boundingBoxes is None:
182 boundingBoxes = spam.label.boundingBoxes(lab1)
184 if centresOfMass is None:
185 centresOfMass = spam.label.centresOfMass(lab1)
187 assert applyF in ["all", "no", "rigid"]
189 # Create pixelSearchCC vector
190 pixelSearchCC = numpy.zeros((nLabels + 1), dtype=float)
191 # Error compatible with register()
192 error = numpy.zeros((nLabels + 1), dtype=float)
193 returnStatus = numpy.ones((nLabels + 1), dtype=int)
194 # deltaPhiNorm = numpy.ones((nLabels+1), dtype=int)
195 iterations = numpy.ones((nLabels + 1), dtype=int)
197 numberOfNodes = len(labelsToCorrelate)
199 # CHECK THIS
200 # firstNode = 1
201 # finishedNodes = 1
202 returnStatus[0] = 0
204 global _multiprocessingPixelSearchOneNodeDiscrete
206 def _multiprocessingPixelSearchOneNodeDiscrete(nodeNumber):
207 """
208 Function to be called by multiprocessing parallelisation for pixel search in one position.
209 This function will call getImagettes, or the equivalent for labels and perform the pixel search
211 Parameters
212 ----------
213 nodeNumber : int
214 node number to work on
216 Returns
217 -------
218 List with:
219 - nodeNumber (needed to write result in right place)
220 - displacement vector
221 - NCC value
222 - error value
223 - return Status
224 """
226 imagetteReturns = spam.label.getImagettesLabelled(
227 lab1,
228 nodeNumber,
229 PhiField[nodeNumber].copy(),
230 im1,
231 im2,
232 searchRange.copy(),
233 boundingBoxes,
234 centresOfMass,
235 margin=labelDilate,
236 labelDilate=labelDilate,
237 applyF=applyF,
238 volumeThreshold=volThreshold,
239 )
240 imagetteReturns["imagette2mask"] = None
242 # If getImagettes was successful (size check and mask coverage check)
243 if imagetteReturns["returnStatus"] == 1:
244 PSreturns = _pixelSearch(
245 imagetteReturns["imagette1"], imagetteReturns["imagette2"], imagette1mask=imagetteReturns["imagette1mask"], imagette2mask=imagetteReturns["imagette2mask"], returnError=True
246 )
247 pixelSearchOffset = imagetteReturns["pixelSearchOffset"]
249 return (nodeNumber, PSreturns[0] + pixelSearchOffset, PSreturns[1], PSreturns[2], imagetteReturns["returnStatus"])
251 # Failed to extract imagettes or something
252 else:
253 return (nodeNumber, numpy.array([numpy.nan] * 3), 0.0, numpy.inf, imagetteReturns["returnStatus"])
255 print("\n\tStarting Pixel Search Discrete (with {} process{})".format(numProc, "es" if numProc > 1 else ""))
257 widgets = [
258 progressbar.FormatLabel(""),
259 " ",
260 progressbar.Bar(),
261 " ",
262 progressbar.AdaptiveETA(),
263 ]
264 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes)
265 pbar.start()
266 finishedNodes = 0
268 with multiprocessing.Pool(processes=numProc) as pool:
269 for returns in pool.imap_unordered(_multiprocessingPixelSearchOneNodeDiscrete, labelsToCorrelate):
270 finishedNodes += 1
272 # Update progres bar if point is not skipped
273 if returns[4] > 0:
274 widgets[0] = progressbar.FormatLabel(" CC={:0>7.5f} ".format(returns[2]))
275 pbar.update(finishedNodes)
277 PhiField[returns[0], 0:3, -1] = returns[1]
278 # Create pixelSearchCC vector
279 pixelSearchCC[returns[0]] = returns[2]
280 error[returns[0]] = returns[3]
281 returnStatus[returns[0]] = returns[4]
282 pool.close()
283 pool.join()
285 pbar.finish()
287 return {
288 "PhiField": PhiField,
289 "pixelSearchCC": pixelSearchCC,
290 "error": error,
291 "returnStatus": returnStatus,
292 # "deltaPhiNorm" : deltaPhiNorm,
293 "iterations": iterations,
294 }
297def pixelSearchLocal(
298 im1,
299 im2,
300 hws,
301 searchRange,
302 nodePositions,
303 PhiField=None,
304 # twoD=False,
305 im1mask=None,
306 im2mask=None,
307 applyF="all",
308 maskCoverage=0.5,
309 greyLowThresh=-numpy.inf,
310 greyHighThresh=numpy.inf,
311 numProc=multiprocessing.cpu_count(),
312):
313 """
314 Performs a series of pixel searches at given coordinates extracting a parallelepiped axis aligned box.
316 Parameters
317 ----------
318 im1 : 2D or 3D numpy array of greylevels
319 Array representing greylevels in the reference configuration
321 im2 : 2D or 3D numpy array of greylevels
322 Array representing greylevels in the reference configuration, same size as im1
324 hws : a 3-element list of ints
325 Z, Y, X extents of the "Half-window" size for extraction of the correlation window around the pixel of interest defined in nodePositions
327 searchRange : 1D numpy array of signed ints
328 Array defining search range in [low Z, high Z, low Y, high Y, low X, high X]
330 nodePositions : N x 3 numpy array of ints
331 Z, Y, X positions of points to correlate, around which to extract the HWS box
333 PhiField : N x 4 x 4 numpy array of floats, optional
334 Optional initial guess for Phi for each correlation point.
335 Default = 4x4 identity matrix
337 im1mask : 3D numpy array of bools, optional
338 Binary array, same size as im1, which is True in for pixels in im1 that should be included in the correlation
339 Default = all the pixels in im1
341 im2mask : 3D numpy array of bools, optional
342 Binary array, same size as im1, which is True in for pixels in im2 that should be included in the correlation
343 Default = all the pixels in im2
345 applyF : string, optional
346 Apply the F part of Phi guess? Accepted values are:\n\t"all": apply all of F' + '\n\t"rigid": apply rigid part (mostly rotation) \n\t"no": don\'t apply it "all" is default
348 maskCoverage : float, optional
349 For a given correlation window, what fraction of pixels should be within the mask (if passed) to consider the correlation window for correlation,
350 otherwise it will be skipped with return status = -5 as per spam.DIC.grid.getImagettes().
351 Default = 0.5
353 greyLowThresh : float, optional
354 Mean grey value in im1 above which a correlation window will be correlated, below this value it will be skipped.
355 Default = -numpy.inf
357 greyHighThresh : float, optional
358 Mean grey value in im1 below which a correlation window will be correlated, above this value it will be skipped.
359 Default = numpy.inf
361 numProc : int, optional
362 Number of processes to use for the calculation, default = multiprocessing.cpu_count()
364 Returns
365 -------
366 Dictionary including the following keys for all window, but only filled in correlated window:
367 - PhiField : L+1 x 4 x 4 numpy array of floats: Phi for each window
368 - pixelSearchCC: L+1 numpy array of floats: Correlation Coefficient for each window bounded [0, 1]
369 - error: L+1 numpy array of floats: SSQD per pixel for each window
370 - returnStatus: L+1 numpy array of ints: returnStatus for each window
371 - deltaPhiNorm: L+1 numpy array of floats: 1 for each window TODO decide if this is worth doing here
372 - iterations: L+1 numpy array of ints: 1 for each window TODO decide if this is worth doing here
373 """
374 numberOfNodes = nodePositions.shape[0]
375 if len(im1.shape) == 2:
376 twoD = True
377 else:
378 twoD = False
380 # Create pixelSearchCC vector
381 pixelSearchCC = numpy.zeros((numberOfNodes), dtype=float)
382 # Error compatible with register()
383 error = numpy.zeros((numberOfNodes), dtype=float)
384 returnStatus = numpy.ones((numberOfNodes), dtype=int)
385 deltaPhiNorm = numpy.ones((numberOfNodes), dtype=int)
387 iterations = numpy.ones((numberOfNodes), dtype=int)
389 firstNode = 0
390 finishedNodes = 0
392 global _multiprocessingPixelSearchOneNodeLocal
394 def _multiprocessingPixelSearchOneNodeLocal(nodeNumber):
395 """
396 Function to be called by multiprocessing parallelisation for pixel search in one position.
397 This function will call getImagettes, or the equivalent for labels and perform the pixel search
399 Parameters
400 ----------
401 nodeNumber : int
402 node number to work on
404 Returns
405 -------
406 List with:
407 - nodeNumber (needed to write result in right place)
408 - displacement vector
409 - NCC value
410 - error value
411 - return Status
412 """
413 imagetteReturns = spam.DIC.getImagettes(
414 im1,
415 nodePositions[nodeNumber],
416 hws,
417 PhiField[nodeNumber].copy(),
418 im2,
419 searchRange.copy(),
420 im1mask=im1mask,
421 im2mask=im2mask,
422 minMaskCoverage=maskCoverage,
423 greyThreshold=[greyLowThresh, greyHighThresh],
424 applyF=applyF,
425 twoD=twoD,
426 )
428 # If getImagettes was successful (size check and mask coverage check)
429 if imagetteReturns["returnStatus"] == 1:
430 PSreturns = _pixelSearch(
431 imagetteReturns["imagette1"], imagetteReturns["imagette2"], imagette1mask=imagetteReturns["imagette1mask"], imagette2mask=imagetteReturns["imagette2mask"], returnError=True
432 )
433 pixelSearchOffset = imagetteReturns["pixelSearchOffset"]
435 return (nodeNumber, PSreturns[0] + pixelSearchOffset, PSreturns[1], PSreturns[2], imagetteReturns["returnStatus"])
437 # Failed to extract imagettes or something
438 else:
439 return (nodeNumber, numpy.array([numpy.nan] * 3), 0.0, numpy.inf, imagetteReturns["returnStatus"])
441 print("\n\tStarting Pixel Search Local (with {} process{})".format(numProc, "es" if numProc > 1 else ""))
443 widgets = [
444 progressbar.FormatLabel(""),
445 " ",
446 progressbar.Bar(),
447 " ",
448 progressbar.AdaptiveETA(),
449 ]
450 pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes)
451 pbar.start()
452 # finishedNodes = 0
454 with multiprocessing.Pool(processes=numProc) as pool:
455 for returns in pool.imap_unordered(_multiprocessingPixelSearchOneNodeLocal, range(firstNode, numberOfNodes)):
456 finishedNodes += 1
458 # Update progres bar if point is not skipped
459 if returns[4] > 0:
460 widgets[0] = progressbar.FormatLabel(" CC={:0>7.5f} ".format(returns[2]))
461 pbar.update(finishedNodes)
463 PhiField[returns[0], 0:3, -1] = returns[1]
464 # Create pixelSearchCC vector
465 pixelSearchCC[returns[0]] = returns[2]
466 error[returns[0]] = returns[3]
467 returnStatus[returns[0]] = returns[4]
468 pool.close()
469 pool.join()
471 pbar.finish()
473 return PhiField, pixelSearchCC, error, returnStatus, deltaPhiNorm, iterations