spam-ldic 18.9 KB
Newer Older
1
#!/usr/bin/env python
2

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
"""
This python script performs Local Digital Image Correlation using SPAM functions
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/>.
"""

21
import os
22
23
os.environ['OPENBLAS_NUM_THREADS'] = '1'

24
25
26
import spam.DIC
import spam.deformation
import spam.helpers
27
import argparse
28
29
30

import numpy
import scipy.ndimage
31
32
33
34
import tifffile
import multiprocessing
import progressbar

35
36
37
38
39





Edward Andò's avatar
Edward Andò committed
40
# Define argument parser object
Edward Andò's avatar
Edward Andò committed
41
42
43
44
45
46
parser = argparse.ArgumentParser(description="spam-ldic "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
                                             "This script performs Local Digital Image Correlation script between a series of at least two 3D greyscale images"+\
                                             "with independent measurement points spread on a regular grid (with -ns spacing in pixels between points). "+\
                                             "Around each point a cubic subvolume of +-hws (Half-window size) is extracted and correlated"+\
                                             "\nSee for more details: https://ttk.gricad-pages.univ-grenoble-alpes.fr/spam/tutorial-02b-DIC-practice.html",
                                 formatter_class=argparse.RawTextHelpFormatter)
47

Edward Andò's avatar
Edward Andò committed
48
# Parse arguments with external helper function
49
args = spam.helpers.optionsParser.ldicParser(parser)
50

Emmanuel Roubin's avatar
Emmanuel Roubin committed
51
if len(args.inFiles) < 2:
52
    print("\nldic: Did not receive enough input images... you need (at least) two to tango...")
Edward Andò's avatar
Edward Andò committed
53
    exit()
54

55
print("spam-ldic -- Current Settings:")
56
57
58
59
60
61
62
63
64
65
66
67
argsDict = vars(args)
for key in sorted(argsDict):
    print("\t{}: {}".format(key, argsDict[key]))

# Load reference image
im1 = tifffile.imread(args.inFiles[0].name)

# Detect unpadded 2D image first:
if len(im1.shape) == 2:
    im1 = im1[numpy.newaxis, ...]
if im1.shape[0] == 1:
    twoD = True
Edward Andò's avatar
Edward Andò committed
68
else:
69
    twoD = False
70

71
72
73
74
75
76
if args.MASK1:
    im1mask = tifffile.imread(args.MASK1.name) != 0
    if len(im1mask.shape) == 2:
        im1mask = im1mask[numpy.newaxis, ...]
else:
    im1mask = None
77

78

79
80
81
82
83
84
85
86
87
88
89
90
91
92
### Interpolation settings
if args.INTERPOLATION_ORDER == 1:
    interpolator = 'C'
else:
    interpolator = 'python'
# Override interpolator for python in 2D
if twoD:
    interpolator = 'python'

margin = [-args.MARGIN[0], args.MARGIN[0],
          -args.MARGIN[1], args.MARGIN[1],
          -args.MARGIN[2], args.MARGIN[2]]


93
firstCorrelation = True
94

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
# Bad to redifine this for every loop, so it's defined here, to be called by the pool
def correlateOneNode(nodeNumber):
    """
    This function does a correlation at one point and returns:

    Returns
    -------
        List of:
        - nodeNumber
        - Phi
        - returnStatus
        - error
        - iterations
        - deltaPhiNorm
    """
    PhiInit = PhiField[nodeNumber]

112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    if numpy.isfinite(PhiInit).sum() == 16:
        imagetteReturns = spam.DIC.getImagettes(im1, nodePositions[nodeNumber], args.HWS, PhiInit.copy(), im2, margin, im1mask=im1mask, minMaskCoverage=args.MASK_COVERAGE, greyThreshold=[args.GREY_LOW_THRESH, args.GREY_HIGH_THRESH], applyF='no', twoD=twoD)

        if imagetteReturns['returnStatus'] == 1:
            # compute displacement that will be taken by the getImagettes
            initialDisplacement = numpy.round(PhiInit[0:3, 3]).astype(int)
            PhiInit[0:3,-1] -= initialDisplacement

            #print(imagetteReturns['imagette1'].shape)
            #print(imagetteReturns['imagette2'].shape)
            #print()

            registerReturns = spam.DIC.register(imagetteReturns['imagette1'],
                                                imagetteReturns['imagette2'],
                                                im1mask=imagetteReturns['imagette1mask'],
                                                PhiInit=PhiInit, # minus initial displacement above, which is in the search range and thus taken into account in imagette2
                                                margin=1,  # see top of this file for compensation
                                                maxIterations=args.MAX_ITERATIONS,
                                                deltaPhiMin=args.MIN_DELTA_PHI,
                                                updateGradient=args.UPDATE_GRADIENT,
                                                interpolationOrder=args.INTERPOLATION_ORDER,
                                                verbose=False,
                                                imShowProgress=False)
            goodPhi = registerReturns['Phi']
            goodPhi[0:3,-1] += initialDisplacement
            return nodeNumber, goodPhi, registerReturns['returnStatus'], registerReturns['error'], registerReturns['iterations'], registerReturns['deltaPhiNorm']

        else:
            badPhi = numpy.eye(4)
            badPhi[0:3, 3] = numpy.nan
            return nodeNumber, badPhi, imagetteReturns['returnStatus'], numpy.inf, 0, numpy.inf
143
    else:
144
        ### Phi has nans or infs
145
146
        badPhi = numpy.eye(4)
        badPhi[0:3, 3] = numpy.nan
147
        return nodeNumber, badPhi, -7, numpy.inf, 0, numpy.inf
148

149

150
151
152
153
154
155
156
157
158
# Loop over input images
for im2number in range(1, len(args.inFiles)):
    # Variables to track last correlation in order to ask MPI workers to hang up
    if im2number == len(args.inFiles)-1: lastCorrelation = True
    else: lastCorrelation = False

    # decide on number, in input files list, of the reference image
    if args.SERIES_INCREMENTAL:
        im1number = im2number-1
159
    else:
160
        im1number = 0
161

162
163
164
    # Output file name prefix
    if args.PREFIX is None or len(args.inFiles) > 2:
        args.PREFIX = os.path.splitext(os.path.basename(args.inFiles[im1number].name))[0]+"-"+os.path.splitext(os.path.basename(args.inFiles[im2number].name))[0]
165

166
167
168
    # If not first correlation and we're interested in loading previous Ffile:
    if not firstCorrelation and args.SERIES_PHIFILE:
        args.PHIFILE = previousPhiFile
169

170
    print("\nCorrelating:", args.PREFIX)
171

172
173
174
    im2 = tifffile.imread(args.inFiles[im2number].name)
    if len(im2.shape) == 2:
        im2 = im2[numpy.newaxis, ...]
175

176
177
178
179
    assert(im1.shape == im2.shape), "\nim1 and im2 must have the same size! Exiting."
    if args.MASK1:
        assert(im1.shape == im1mask.shape), "\nim1 and im1mask must have the same size! Exiting."

180
181
182
183
    # Three cases to handle:
    #   1. phi file is reg   -> define nodes and apply reg
    #   2. phi file is field -> take everything and check NS if passed
    #   3. no phi file       -> define nodes
184
    if args.PHIFILE is not None:
185
        PhiFromFile = spam.helpers.readCorrelationTSV(args.PHIFILE.name, fieldBinRatio=args.PHIFILE_BIN_RATIO)
186
187
188
        if PhiFromFile is None:
            print(f"\tFailed to read your TSV file passed with -pf {args.PHIFILE.name}")
            exit()
189

190
191
        # If the read Phi-file has only one line -- it's a single point registration!
        if PhiFromFile['fieldCoords'].shape[0] == 1:
192
193
            PhiInit = PhiFromFile['PhiField'][0]
            print(f"\tI read a registration from a file in binning {args.PHIFILE_BIN_RATIO}")
194

195
            decomposedPhiInit = spam.deformation.decomposePhi(PhiInit)
196
            print("\tTranslations (px)")
197
            print("\t\t", decomposedPhiInit['t'])
198
            print("\tRotations (deg)")
199
200
201
            print("\t\t", decomposedPhiInit['r'])
            print("\tZoom")
            print("\t\t", decomposedPhiInit['z'])
202
203
            del decomposedPhiInit

204
205
206
207
208
209
            # Create nodes
            if args.NS is None:
                print(f"spam-ldic: You passed a registration file {args.PHIFILE.name}, I need -ns to be defined")
                exit()
            nodePositions, nodesDim = spam.DIC.makeGrid(im1.shape, args.NS)
            numberOfNodes = nodePositions.shape[0]
210
211
212
213
214

            # We have a registration to apply to all points.
            # This is done in 2 steps:
            #   1. by copying the registration's little F to the Fs of all points
            #   2. by calling the decomposePhi function to compute the translation of each point
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
            # Now that we know how many points we want to correlate, initalise PhiField
            PhiField = numpy.zeros((numberOfNodes, 4, 4))
            for node in range(numberOfNodes):
                if args.APPLY_F == "all":
                    PhiField[node] = PhiInit.copy()
                elif args.APPLY_F == "rigid":
                    PhiField[node] = spam.deformation.computeRigidPhi(PhiInit.copy())
                else:
                    PhiField[node] = numpy.eye(4)
                # trasnlation application
                PhiField[node][0:3, -1] = spam.deformation.decomposePhi(PhiInit.copy(), PhiCentre=PhiFromFile["fieldCoords"][0], PhiPoint=nodePositions[node])["t"]

        else: # we have a Phi field and not a registration
            nodePositions   = PhiFromFile["fieldCoords"]
            numberOfNodes   = nodePositions.shape[0]
Edward Andò's avatar
Edward Andò committed
230
            nodesDim        = PhiFromFile["fieldDims"]
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
            nodeSpacingFile = numpy.array([numpy.unique(nodePositions[:, i])[1] - numpy.unique(nodePositions[:, i])[0] if len(numpy.unique(nodePositions[:, i])) > 1 else numpy.unique(nodePositions[:, i])[0] for i in range(3)])
            PhiField        = PhiFromFile["PhiField"]

            # In case NS is also defined, complain, but if it's the same as the loaded data, continue
            if args.NS is not None:
                # compare them
                if not numpy.allclose(numpy.array(args.NS), nodeSpacingFile, atol=0.0):
                    print(f"spam-ldic: you passed a -ns={args.NS} which contradicts the node spacing in your Phi Field TSV of {nodeSpacingFile}")
                    print(f"\thint 1: if you pass a Phi Field TSV you don't need to also define the node spacing")
                    print(f"\thint 2: if you want to use your Phi Field TSV {args.PHIFILE.name} on a finer node spacing, pass it with spam-passPhiField")
                    exit()
                else:
                    print(f"spam-ldic: passing -ns with a Phi Field TSV is not needed")

            # If it's compatible, update args.NS
            args.NS = nodeSpacingFile
247

248
    else: # No Phi file at all
249
        if args.NS is None:
250
            print("spam-ldic: I don't have a phi file or -ns defined, so don't know how to define grid...")
251
252
            exit()
        nodePositions, nodesDim = spam.DIC.makeGrid(im1.shape, args.NS)
253
        numberOfNodes = nodePositions.shape[0]
254
255
256
257
258

        PhiField = numpy.zeros((numberOfNodes, 4, 4))
        for node in range(numberOfNodes):
            PhiField[node] = numpy.eye(4)

259
260
    finishedNodes = 0

261
262
263
264
    error = numpy.zeros((numberOfNodes))
    iterations = numpy.zeros((numberOfNodes))
    returnStatus = numpy.zeros((numberOfNodes))
    deltaPhiNorm = numpy.zeros((numberOfNodes))
265

266
267
268
    if args.PROCESSES is None: args.PROCESSES = multiprocessing.cpu_count()

    print("\n\tStarting local dic (with {} process{})".format(args.PROCESSES, 'es' if args.PROCESSES > 1 else ''))
269
270
271
272

    widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
    pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes)
    pbar.start()
273
    finishedNodes = 0
274

275
276
    with multiprocessing.Pool(processes=args.PROCESSES) as pool:
        for returns in pool.imap_unordered(correlateOneNode, range(numberOfNodes)):
277
            finishedNodes += 1
278
279
280
281
282
283
284
285
286
287
288

            # Update progres bar if point is not skipped
            if returns[2] > 0:
                widgets[0] = progressbar.FormatLabel("  it={:0>3d}  dPhiNorm={:0>6.4f}  rs={:+1d} ".format(returns[4], returns[5], returns[2]))
                pbar.update(finishedNodes)
            nodeNumber = returns[0]
            PhiField[nodeNumber]     = returns[1]
            returnStatus[nodeNumber] = returns[2]
            error[nodeNumber]        = returns[3]
            iterations[nodeNumber]   = returns[4]
            deltaPhiNorm[nodeNumber] = returns[5]
289
290

    pbar.finish()
291

292
    print("\n")
293

294
295
296
297
298
299
300
301


    ## Finished! Get ready for output.
    #if args.REGSUB:
        #print("\n\tFinished correlations. Subtracting rigid-body motion from displacements of each point")
        #PhiFieldMinusRigid = PhiField.copy()
        #for node in range(nodePositions.shape[0]):
            #PhiFieldMinusRigid[node][0:3,-1] -= rigidDisp[node]
302
303
304
305
306
307
308
309

    if args.TSV:
        # Make one big array for writing:
        #   First the node number,
        #   3 node positions,
        #   F[0:3,0:2]
        #   Pixel-search CC
        #   SubPixError, SubPixIterations, SubPixelReturnStatus
310
        TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\terror\titerations\treturnStatus\tdeltaPhiNorm"
311
312
        outMatrix = numpy.array([numpy.array(range(nodePositions.shape[0])),
                                    nodePositions[:, 0], nodePositions[:, 1], nodePositions[:, 2],
313
314
315
316
317
318
                                    PhiField[:, 0, 0],   PhiField[:, 0, 1],   PhiField[:, 0, 2],    PhiField[:, 0, 3],
                                    PhiField[:, 1, 0],   PhiField[:, 1, 1],   PhiField[:, 1, 2],    PhiField[:, 1, 3],
                                    PhiField[:, 2, 0],   PhiField[:, 2, 1],   PhiField[:, 2, 2],    PhiField[:, 2, 3],
                                    error,               iterations,          returnStatus,         deltaPhiNorm]).T

        numpy.savetxt(args.OUT_DIR+"/"+args.PREFIX+"-ldic.tsv",
319
320
321
322
323
324
                        outMatrix,
                        fmt='%.7f',
                        delimiter='\t',
                        newline='\n',
                        comments='',
                        header=TSVheader)
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
        ## Hold onto that name if we need to reload
        #if args.SERIES_PHIFILE: previousPhiFile = args.OUT_DIR+"/"+args.PREFIX+".tsv"

        #if args.REGSUB:
            #outMatrix = numpy.array([numpy.array(range(nodePositions.shape[0])),
                                        #nodePositions[:, 0], nodePositions[:, 1], nodePositions[:, 2],
                                        #PhiFieldMinusRigid[:, 0, 0],         PhiFieldMinusRigid[:, 0, 1],         PhiFieldMinusRigid[:, 0, 2],    PhiFieldMinusRigid[:, 0, 3],
                                        #PhiFieldMinusRigid[:, 1, 0],         PhiFieldMinusRigid[:, 1, 1],         PhiFieldMinusRigid[:, 1, 2],    PhiFieldMinusRigid[:, 1, 3],
                                        #PhiFieldMinusRigid[:, 2, 0],         PhiFieldMinusRigid[:, 2, 1],         PhiFieldMinusRigid[:, 2, 2],    PhiFieldMinusRigid[:, 2, 3]]).T
            #if args.PS == 'on' or ( registrationSuccessful == False and args.PS == 'auto'):
                #outMatrix = numpy.hstack([outMatrix, numpy.array([PSCC]).T])

            #if args.GRID_POINT:
                #outMatrix = numpy.hstack([outMatrix, numpy.array([error,
                                                                    #iterations,
                                                                    #returnStatus,
                                                                    #deltaPhiNorm]).T])
            #numpy.savetxt(args.OUT_DIR+"/"+args.PREFIX+"-regsub.tsv",
                            #outMatrix,
                            #fmt='%.7f',
                            #delimiter='\t',
                            #newline='\n',
                            #comments='',
                            #header=TSVheader)
349
350

    if args.TIFF:
Edward Andò's avatar
Edward Andò committed
351
        if nodesDim[0] != 1:
352
353
354
355
356
357
358
359
360
361
362
363
            tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-ldic-Zdisp.tif",              PhiField[:, 0, -1].astype('<f4').reshape(nodesDim))
            #if args.REGSUB:
                #tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Zdisp-regsub.tif",   PhiFieldMinusRigid[:, 0, -1].astype('<f4').reshape(nodesDim))
        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-ldic-Ydisp.tif",                  PhiField[:, 1, -1].astype('<f4').reshape(nodesDim))
        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-ldic-Xdisp.tif",                  PhiField[:, 2, -1].astype('<f4').reshape(nodesDim))
        #if args.REGSUB:
            #tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Ydisp-regsub.tif",       PhiFieldMinusRigid[:, 1, -1].astype('<f4').reshape(nodesDim))
            #tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-Xdisp-regsub.tif",       PhiFieldMinusRigid[:, 2, -1].astype('<f4').reshape(nodesDim))
        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-ldic-error.tif",        error.astype('<f4').reshape(nodesDim))
        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-ldic-iterations.tif",   iterations.astype('<f4').reshape(nodesDim))
        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-ldic-returnStatus.tif", returnStatus.astype('<f4').reshape(nodesDim))
        tifffile.imsave(args.OUT_DIR+"/"+args.PREFIX+"-ldic-deltaPhiNorm.tif", deltaPhiNorm.astype('<f4').reshape(nodesDim))
364
365
366
367

    # Collect data into VTK output
    if args.VTK:
        cellData = {}
368
369
370
371
372
373
374
375
376
377
378
379
        cellData['displacements']         = PhiField[:, :-1, 3].reshape((nodesDim[0], nodesDim[1], nodesDim[2], 3))
        #if args.REGSUB:
            #cellData['displacements-regsub'] = PhiFieldMinusRigid[:, :-1, 3].reshape((nodesDim[0], nodesDim[1], nodesDim[2], 3))
        cellData['error']        = error.reshape(nodesDim)
        cellData['iterations']   = iterations.reshape(nodesDim)
        cellData['returnStatus'] = returnStatus.reshape(nodesDim)
        cellData['deltaPhiNorm'] = deltaPhiNorm.reshape(nodesDim)

        cellData['error'       ][numpy.logical_not(numpy.isfinite(cellData['error']))]        = 0
        cellData['iterations'  ][numpy.logical_not(numpy.isfinite(cellData['iterations']))]   = 0
        cellData['returnStatus'][numpy.logical_not(numpy.isfinite(cellData['returnStatus']))] = 0
        cellData['deltaPhiNorm'][numpy.logical_not(numpy.isfinite(cellData['deltaPhiNorm']))] = 0
380
381
382

        # Overwrite nans and infs with 0, rubbish I know
        cellData['displacements'][numpy.logical_not(numpy.isfinite(cellData['displacements']))] = 0
383
384
        #if args.REGSUB:
            #cellData['displacements-regsub'][numpy.logical_not(numpy.isfinite(cellData['displacements-regsub']))] = 0
385

386
387
388
389
        # This is perfect in the case where NS = 2xHWS, these cells will all be in the right place
        #   In the case of overlapping of under use of data, it should be approximately correct
        # If you insist on overlapping, then perhaps it's better to save each point as a cube glyph
        #   and actually *have* overlapping
390
        spam.helpers.writeStructuredVTK(origin=nodePositions[0]-args.HWS, aspectRatio=args.NS, cellData=cellData, fileName=args.OUT_DIR+"/"+args.PREFIX+"-ldic.vtk")
391
    firstCorrelation = False
392

393
394
395
    if args.SERIES_INCREMENTAL:
        # If in incremental mode, current deformed image is next reference image
        im1 = im2.copy()
396