deformationField.py 61.8 KB
Newer Older
Olga Stamati's avatar
Olga Stamati committed
1
"""
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Library of SPAM functions for dealing with fields of Phi or fields of F
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/>.
"""

# 2020-02-24 Olga Stamati and Edward Ando
20
import spam.deformation
21
import numpy
22
23
import progressbar
import spam.label
24
25
26
27
28
29
import multiprocessing

# Global number of processes
nProcessesDefault = multiprocessing.cpu_count()

def FfieldRegularQ8(displacementField, nodeSpacing, nProcesses=nProcessesDefault, verbose=True):
30
31
32
33
34
35
36
37
    """
    This function computes the transformation gradient field F from a given displacement field.
    Please note: the transformation gradient tensor: F = I + du/dx.

    This function computes du/dx in the centre of an 8-node cell (Q8 in Finite Elements terminology) using order one (linear) shape functions.

    Parameters
    ----------
38
        displacementField : 4D array of floats
39
40
41
            The vector field to compute the derivatives.
            #Its shape is (nz, ny, nx, 3)

42
        nodeSpacing : 3-component list of floats
43
44
            Length between two nodes in every direction (*i.e.,* size of a cell)

45
46
47
48
49
50
51
52
        nProcesses : integer, optional
            Number of processes for multiprocessing
            Default = number of CPUs in the system

        verbose : boolean, optional
            Print progress?
            Default = True

53
54
    Returns
    -------
55
        F : (nz-1) x (ny-1) x (nx-1) x 3x3 array of n cells
56
57
58
59
60
61
            The field of the transformation gradient tensors
    """
    # import spam.DIC.deformationFunction
    #import spam.mesh.strain

    # Define dimensions
62
    dims = list(displacementField.shape[0:3])
63
64

    # Q8 has 1 element fewer than the number of displacement points
65
    nCells = [n - 1 for n in dims]
66
67

    # Check if a 2D field is passed
68
    if dims[0] == 1:
69
        # Add a ficticious layer of nodes and cells in Z direction
70
        dims[0] += 1
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
        nCells[0] += 1
        nodeSpacing[0] += 1

        # Add a ficticious layer of equal displacements so that the strain in z is null
        displacementField = numpy.concatenate((displacementField, displacementField))

    # Transformation gradient tensor F = du/dx +I
    Ffield = numpy.zeros((nCells[0], nCells[1], nCells[2], 3, 3))

    # Define the coordinates of the Parent Element
    # we're using isoparametric Q8 elements
    lid = numpy.zeros((8, 3)).astype('<u1')  # local index
    lid[0] = [0, 0, 0]
    lid[1] = [0, 0, 1]
    lid[2] = [0, 1, 0]
    lid[3] = [0, 1, 1]
    lid[4] = [1, 0, 0]
    lid[5] = [1, 0, 1]
    lid[6] = [1, 1, 0]
    lid[7] = [1, 1, 1]

92

93
94
95
96
97
98
99
100
101
102
103
104
105
106
    # Calculate the derivatives of the shape functions
    # Since the center is equidistant from all 8 nodes, each one gets equal weighting
    SFderivative = numpy.zeros((8, 3))
    for node in range(8):
        # (local nodes coordinates) / weighting of each node
        SFderivative[node, 0] = (2.0 * (float(lid[node, 0]) - 0.5)) / 8.0
        SFderivative[node, 1] = (2.0 * (float(lid[node, 1]) - 0.5)) / 8.0
        SFderivative[node, 2] = (2.0 * (float(lid[node, 2]) - 0.5)) / 8.0

    # Compute the jacobian to go from local(Parent Element) to global base
    jacZ = 2.0 / float(nodeSpacing[0])
    jacY = 2.0 / float(nodeSpacing[1])
    jacX = 2.0 / float(nodeSpacing[2])

107
    if verbose: bar = progressbar.ProgressBar(maxval=nCells[0]*nCells[1]*nCells[2]).start()
108
109
    nodesDone = 0

110
111
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
    ## Loop over the cells
    #global computeConvexVolume
    #def computeConvexVolume(label):
        #labelI = spam.label.getLabel(lab, label, boundingBoxes = boundingBoxes, centresOfMass = centresOfMass)
        #subvol = labelI['subvol']
        #points = numpy.transpose(numpy.where(subvol))
        #try:
            #hull = scipy.spatial.ConvexHull(points)
            #deln = scipy.spatial.Delaunay(points[hull.vertices])
            #idx = numpy.stack(numpy.indices(subvol.shape), axis = -1)
            #out_idx = numpy.nonzero(deln.find_simplex(idx) + 1)
            #hullIm = numpy.zeros(subvol.shape)
            #hullIm[out_idx] = 1
            #hullVol = spam.label.volumes(hullIm)
            #return label, hullVol[-1]
        #except:
            #return label, 0

    ## Run multiprocessing
    #with multiprocessing.Pool(processes=nProcesses) as pool:
        #for returns in pool.imap_unordered(computeConvexVolume, range(1, nLabels+1)):
            #if verbose:
                #finishedNodes += 1
                #widgets[0] = progressbar.FormatLabel("  vol={:0>7.5f} ".format(returns[1]))
                #pbar.update(finishedNodes)
            #convexVolume[returns[0]] = returns[1]


138
    for kCell in range(nCells[0]):
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
        for jCell in range(nCells[1]):
            for iCell in range(nCells[2]):
                # Check for nans in one of the 8 nodes of the cell
                dCell = displacementField[kCell:kCell + 2, jCell:jCell + 2, iCell:iCell + 2]
                if not numpy.all(numpy.isfinite(dCell)):
                    Ffield[kCell, jCell, iCell, :, :] = numpy.zeros((3, 3)) * numpy.nan

                # If no nans start the strain calculation
                else:
                    # Initialise the gradient of the displacement tensor
                    dudx = numpy.zeros((3, 3))

                    # Loop over each node of the cell
                    for node in range(8):
                        # Get the displacement value
                        d = displacementField[int(kCell + lid[node, 0]), int(jCell + lid[node, 1]), int(iCell + lid[node, 2]), :]

                        # Compute the influence of each node to the displacement gradient tensor
                        dudx[0, 0] += jacZ * SFderivative[node, 0] * d[0]
                        dudx[1, 1] += jacY * SFderivative[node, 1] * d[1]
                        dudx[2, 2] += jacX * SFderivative[node, 2] * d[2]
                        dudx[1, 0] += jacY * SFderivative[node, 1] * d[0]
                        dudx[0, 1] += jacZ * SFderivative[node, 0] * d[1]
                        dudx[2, 1] += jacX * SFderivative[node, 2] * d[1]
                        dudx[1, 2] += jacY * SFderivative[node, 1] * d[2]
                        dudx[2, 0] += jacX * SFderivative[node, 2] * d[0]
                        dudx[0, 2] += jacZ * SFderivative[node, 0] * d[2]

167
168
                    # Adding a transpose to dudx, it's ugly but allows us to pass #142
                    Ffield[kCell, jCell, iCell, :, :] = numpy.eye(3)+dudx.T
169
170
171
                bar.update(nodesDone)
                nodesDone += 1
    bar.finish()
172

173
174
    return Ffield

175

Edward Andò's avatar
Edward Andò committed
176
def FfieldRegularGeers(displacementField, nodeSpacing, neighbourRadius=1.5, mask=True, bruteForceDistance=False, nProcesses=nProcessesDefault, verbose=True):
177
178
179
180
181
182
    """
    This function computes the transformation gradient field F from a given displacement field.
    Please note: the transformation gradient tensor: F = I + du/dx.

    This function computes du/dx as a weighted function of neighbouring points.
    Here is implemented the linear model proposed in:
183
    "Computing strain fields from discrete displacement fields in 2D-solids", Geers et al., 1996
184
185
186

    Parameters
    ----------
187
        displacementField : 4D array of floats
188
189
190
            The vector field to compute the derivatives.
            Its shape is (nz, ny, nx, 3).

191
        nodeSpacing : 3-component list of floats
192
193
            Length between two nodes in every direction (*i.e.,* size of a cell)

194
        neighbourRadius : float, optional
195
            Distance in nodeSpacings to include neighbours in the strain calcuation.
196
197
198
199
200
201
202
203
204
205
            Default = 1.5*nodeSpacing which will give radius = 1.5*min(nodeSpacing)

        mask : bool, optional
            Avoid non-correlated NaN points in the displacement field?
            Default = True

        bruteForceDistance : bool, optional
            Use scipy.spatial.KDtree.query_ball_point to obtain neighbours (bruteForceDistance=False),
            or explicitly compute distance to each point (bruteForceDistance=True).
            Default = False
206

Edward Andò's avatar
Edward Andò committed
207
208
209
210
211
212
213
214
        nProcesses : integer, optional
            Number of processes for multiprocessing
            Default = number of CPUs in the system

        verbose : boolean, optional
            Print progress?
            Default = True

215
216
217
218
219
220
221
222
223
    Returns
    -------
        Ffield: nz x ny x nx x 3x3 array of n cells
            The field of the transformation gradient tensors

    Note
    ----
        Taken from the implementation in "TomoWarp2: A local digital volume correlation code", Tudisco et al., 2017
    """
224
    import scipy.spatial
225

226
    ##Define dimensions
227
228
229
    dims = displacementField.shape[0:3]
    nNodes = dims[0]*dims[1]*dims[2]
    displacementFieldFlat = displacementField.reshape(nNodes, 3)
230

231
232
233
234
    # Check if a 2D field is passed
    twoD = False
    if dims[0] == 1: twoD = True

235
236
    #Deformation gradient tensor F = du/dx +I
    Ffield = numpy.zeros((dims[0], dims[1], dims[2], 3, 3))
237
238
    FfieldFlat = numpy.zeros((nNodes, 3, 3))

239
240
241
242
243
244
245
246
    if twoD:
        fieldCoordsFlat = numpy.mgrid[0:1:1,
                                      nodeSpacing[1]:dims[1]*nodeSpacing[1]+nodeSpacing[1]:nodeSpacing[1],
                                      nodeSpacing[2]:dims[2]*nodeSpacing[2]+nodeSpacing[2]:nodeSpacing[2]].reshape(3, nNodes).T
    else:
        fieldCoordsFlat = numpy.mgrid[nodeSpacing[0]:dims[0]*nodeSpacing[0]+nodeSpacing[0]:nodeSpacing[0],
                                      nodeSpacing[1]:dims[1]*nodeSpacing[1]+nodeSpacing[1]:nodeSpacing[1],
                                      nodeSpacing[2]:dims[2]*nodeSpacing[2]+nodeSpacing[2]:nodeSpacing[2]].reshape(3, nNodes).T
247

248
    #Get non-nan displacements
249
250
251
252
253
254
255
256
257
258
    if mask:
        goodPointsMask = numpy.isfinite(displacementField[:,:,:,0].reshape(nNodes))
        badPointsMask  = numpy.isnan(   displacementField[:,:,:,0].reshape(nNodes))
        #Flattened variables
        fieldCoordsFlatGood       = fieldCoordsFlat[goodPointsMask]
        displacementFieldFlatGood = displacementFieldFlat[goodPointsMask]
        #set bad points to nan
        FfieldFlat[badPointsMask] = numpy.eye(3) * numpy.nan
    else:
        #Flattened variables
Edward Andò's avatar
Edward Andò committed
259
        goodPointsMask = numpy.ones(nNodes, dtype=bool)
260
261
        fieldCoordsFlatGood       = fieldCoordsFlat
        displacementFieldFlatGood = displacementFieldFlat
262

263
    #build KD-tree for neighbour identification
Edward Andò's avatar
Edward Andò committed
264
265
    if not bruteForceDistance:
        treeCoord = scipy.spatial.KDTree(fieldCoordsFlatGood)
266

Edward Andò's avatar
Edward Andò committed
267
    # Output array for good points
268
    FfieldFlatGood = numpy.zeros_like(FfieldFlat[goodPointsMask])
Edward Andò's avatar
Edward Andò committed
269
270
271
272

    # Function for parallel mode
    global geersOnePoint
    def geersOnePoint(goodPoint):
273
274
275
        #This is for the linear model, equation 15 in Geers
        centralNodePosition       = fieldCoordsFlatGood[goodPoint]
        centralNodeDisplacement   = displacementFieldFlatGood[goodPoint]
276
277
        sX0X0 = numpy.zeros((3, 3))
        sX0Xt = numpy.zeros((3, 3))
278
279
280
        m0    = numpy.zeros((3))
        mt    = numpy.zeros((3))

281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
        if bruteForceDistance:
            ## Option 1: Manual calculation of distance
            ## superslow alternative to KDtree
            distances = numpy.sqrt(numpy.sum(numpy.square(fieldCoordsFlatGood - centralNodePosition), axis=1))
            ind = numpy.where(distances < neighbourRadius*max(nodeSpacing))[0]

        else:
            ## Option 2: KDTree on distance
            # KD-tree will always give the current point as zero-distance
            ind = treeCoord.query_ball_point(centralNodePosition, neighbourRadius*max(nodeSpacing))

        # We know that the current point will also be included, so remove it from the index list.
        ind = numpy.array(ind)
        ind = ind[ind != goodPoint]
        nNeighbours = len(ind)
296
297
298
299
        nodalRelativePositionsRef = numpy.zeros((nNeighbours, 3)) # Delta_X_0 in paper
        nodalRelativePositionsDef = numpy.zeros((nNeighbours, 3)) # Delta_X_t in paper

        for neighbour, i in enumerate(ind):
300
301
302
            # Relative position in reference configuration
            #                                         absolute position of this neighbour node
            #                                                                  minus abs position of central node
303
304
            nodalRelativePositionsRef[neighbour, :] = fieldCoordsFlatGood[i] - centralNodePosition

305
306
307
308
309
            # Relative position in deformed configuration (i.e., plus displacements)
            #                                         absolute position of this neighbour node
            #                                                                  plus displacement of this neighbour node
            #                                                                                                minus abs position of central node
            #                                                                                                                       minus displacement of central node
310
            nodalRelativePositionsDef[neighbour, :] = fieldCoordsFlatGood[i] + displacementFieldFlatGood[i] - centralNodePosition - centralNodeDisplacement
311
312
313

            for u in range(3):
                for v in range(3):
314
315
                    #sX0X0[u, v] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsRef[neighbour, v]
                    #sX0Xt[u, v] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsDef[neighbour, v]
Edward Andò's avatar
Edward Andò committed
316
                    # Proposed solution for #142 for direction of rotation
317
318
                    sX0X0[v, u] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsRef[neighbour, v]
                    sX0Xt[v, u] += nodalRelativePositionsRef[neighbour, u] * nodalRelativePositionsDef[neighbour, v]
319
320
321
322

            m0 += nodalRelativePositionsRef[neighbour, :]
            mt += nodalRelativePositionsDef[neighbour, :]

323
324
        sX0X0 = nNeighbours*sX0X0
        sX0Xt = nNeighbours*sX0Xt
325

326
327
        A = sX0X0 - numpy.dot(m0, m0)
        C = sX0Xt - numpy.dot(m0, mt)
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
        F = numpy.zeros((3, 3))

        if twoD:
            A = A[1:, 1:]
            C = C[1:, 1:]
            try:
                F[1:, 1:] = numpy.dot(numpy.linalg.inv(A), C)
                F[0, 0] = 1.0
            except numpy.linalg.linalg.LinAlgError:
                #print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
                pass
        else:
            try:
                F = numpy.dot(numpy.linalg.inv(A), C)
            except numpy.linalg.linalg.LinAlgError:
                #print("spam.deformation.deformationField.FfieldRegularGeers(): LinAlgError: A", A)
                pass
345

Edward Andò's avatar
Edward Andò committed
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
        return goodPoint, F

    # Iterate through flat field of Fs
    if verbose: pbar = progressbar.ProgressBar(maxval=fieldCoordsFlatGood.shape[0]).start()

    finishedPoints = 0

    # Run multiprocessing filling in FfieldFlatGood, which will then update FfieldFlat
    with multiprocessing.Pool(processes=nProcesses) as pool:
        for returns in pool.imap_unordered(geersOnePoint, range(fieldCoordsFlatGood.shape[0])):
            if verbose:
                finishedPoints += 1
                pbar.update(finishedPoints)
            FfieldFlatGood[returns[0]] = returns[1]

    if verbose: pbar.finish()
362

363
    FfieldFlat[goodPointsMask] = FfieldFlatGood
364
    return FfieldFlat.reshape(dims[0], dims[1], dims[2], 3, 3)
365

366

367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
def FfieldBagi(points, connectivity, displacements):
    """
    Calculates transformation gradient function using Bagi's 1996 paper, especially equation 3 on page 174.
    Required inputs are connectivity matrix for tetrahedra (for example from spam.mesh.triangulate) and
    nodal positions in reference and deformed configurations.

    Parameters
    ----------
        points : m x 3 numpy array
            M Particles' points in reference configuration

        connectivity : n x 4 numpy array
            Delaunay triangulation connectivity generated by spam.mesh.triangulate for example

        displacements : m x 3 numpy array
            M Particles' displacement

    Returns
    -------
        Ffield: nx3x3 array of n cells
            The field of the transformation gradient tensors
    """
    from progressbar import ProgressBar
    import spam.mesh
    pbar = ProgressBar()

    Ffield = numpy.zeros([connectivity.shape[0], 3, 3], dtype='<f4')

    connectivity = connectivity.astype(numpy.uint)

    # define dimension
    D = 3.0

    # Import modules

    # Construct 4-list of 3-lists of combinations constituting a face of the tet
    combs = [[0, 1, 2],
             [1, 2, 3],
             [2, 3, 0],
             [0, 1, 3]]
    unode = [3, 0, 1, 2]

    # Precompute tetrahedron Volumes
    tetVolumes = spam.mesh.tetVolumes(points, connectivity)

    # Initialize arrays for tet strains
    # print("spam.mesh.bagiStrain(): Constructing strain from Delaunay and Displacements")

    # Loop through tetrahdra to get avec1, uPos1
    # for tet in range(connectivity.shape[0]):
    for tet in pbar(range(connectivity.shape[0])):
        # Get the list of IDs, centroids, center of tet
        tet_ids = connectivity[tet, :]
        # 2019-10-07 EA: Skip references to missing particles
        if max(tet_ids) >= points.shape[0]:
            print("spam.mesh.unstructured.bagiStrain(): this tet has node > points.shape[0], skipping")
            pass
        else:
            tetCoords = points[tet_ids, :]
            tetDisp = displacements[tet_ids, :]
            tetCen = numpy.average(tetCoords, axis=0)
            if numpy.isfinite(tetCoords).sum() + numpy.isfinite(tetDisp).sum() != 3*4*2:
                print("spam.mesh.unstructured.bagiStrain(): nans in position or displacement, skipping")
                # Compute strains
                Ffield[tet] = numpy.zeros((3,3))*numpy.nan
            else:
                # Loop through each face of tet to get avec, upos (Bagi, 1996, pg. 172)
                # aVec1 = numpy.zeros([4, 3], dtype='<f4')
                # uPos1 = numpy.zeros([4, 3], dtype='<f4')
                # uPos2 = numpy.zeros([4, 3], dtype='<f4')
                dudx = numpy.zeros((3, 3), dtype='<f4')

                for face in range(4):
                    faceNorm = numpy.cross(tetCoords[combs[face][0]] - tetCoords[combs[face][1]],
                                           tetCoords[combs[face][0]] - tetCoords[combs[face][2]])

                    # Get a norm vector to face point towards center of tet
                    faceCen = numpy.average(tetCoords[combs[face]], axis=0)
                    tmpnorm = faceNorm / (numpy.linalg.norm(faceNorm))
                    facetocen = tetCen - faceCen
                    if (numpy.dot(facetocen, tmpnorm) < 0):
                        tmpnorm = -tmpnorm

                    # Divide by 6 (1/3 for 1/Dimension; 1/2 for area from cross product)
                    # See first eqn., Bagi, 1996, pg. 172.
                    # aVec1[face] = tmpnorm*numpy.linalg.norm(faceNorm)/6

                    # Undeformed positions
                    # uPos1[face] = tetCoords[unode[face]]
                    # Deformed positions
                    # uPos2[face] = tetComs2[unode[face]]

                    dudx += numpy.tensordot(tetDisp[unode[face]], tmpnorm * numpy.linalg.norm(faceNorm) / 6, axes=0)

                dudx /= float(tetVolumes[tet])

                Ffield[tet] = numpy.eye(3) + dudx
    return Ffield

466

467
def decomposeFfield(Ffield, components, twoD=False, nProcesses=nProcessesDefault, verbose=True):
468
469
470
471
472
473
    """
    This function takes in an F field (from either FfieldRegularQ8, FfieldRegularGeers, FfieldBagi) and
    returns fields of desired transformation components.

    Parameters
    ----------
474
        Ffield : multidimensional x 3 x 3 numpy array of floats
475
476
477
478
479
            Spatial field of Fs

        components : list of strings
            These indicate the desired components consistent with spam.deformation.decomposeF or decomposePhi

480
481
482
483
        twoD : bool, optional
            Is the Ffield in 2D? This changes the strain calculation.
            Default = False

484
485
486
487
488
489
490
491
        nProcesses : integer, optional
            Number of processes for multiprocessing
            Default = number of CPUs in the system

        verbose : boolean, optional
            Print progress?
            Default = True

492
493
    Returns
    -------
494
495
496
497
498
499
        Dictionary containing appropriately reshaped fields of the transformation components requested.

        Keys:
            vol, dev, volss, devss are scalars
            r and z are 3-component vectors
            e and U are 3x3 tensors
500
501
    """
    # The last two are for the 3x3 F field
502
    fieldDimensions  = Ffield.shape[0:-2]
503
    fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
504
    FfieldFlat       = Ffield.reshape(fieldRavelLength, 3, 3)
505
506
507

    output = {}
    for component in components:
508
509
510
511
        if component == 'vol' or component == 'dev' or component == 'volss' or component == 'devss':
            output[component] = numpy.zeros(fieldRavelLength)
        if component == 'r' or component == 'z':
            output[component] = numpy.zeros((fieldRavelLength, 3))
Olga Stamati's avatar
Olga Stamati committed
512
        if component == 'U' or component == 'e':
513
            output[component] = numpy.zeros((fieldRavelLength, 3, 3))
514

515
516
517
    # Function for parallel mode
    global decomposeOneF
    def decomposeOneF(n):
518
        F = FfieldFlat[n]
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
        if numpy.isfinite(F).sum() == 9:
            decomposedF = spam.deformation.decomposeF(F, twoD=twoD)
            return n, decomposedF
        else:
            return n, {'r': numpy.array([numpy.nan] * 3),
                      'z': numpy.array([numpy.nan] * 3),
                      'U': numpy.eye(3) * numpy.nan,
                      'e': numpy.eye(3) * numpy.nan,
                      'vol': numpy.nan,
                      'dev': numpy.nan,
                      'volss': numpy.nan,
                      'devss': numpy.nan}

    # Iterate through flat field of Fs
    if verbose: pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()

    finishedPoints = 0

    # Run multiprocessing
    with multiprocessing.Pool(processes=nProcesses) as pool:
        for returns in pool.imap_unordered(decomposeOneF, range(fieldRavelLength)):
            if verbose:
                finishedPoints += 1
                pbar.update(finishedPoints)
            for component in components:
                output[component][returns[0]] = returns[1][component]

    if verbose: pbar.finish()
547
548
549
550
551
552
553

    # Reshape on the output
    for component in components:
        if component == 'vol' or component == 'dev' or component == 'volss' or component == 'devss':
            output[component] = numpy.array(output[component]).reshape(fieldDimensions)

        if component == 'r' or component == 'z':
554
            output[component] = numpy.array(output[component]).reshape(Ffield.shape[0:-1])
555

Olga Stamati's avatar
Olga Stamati committed
556
        if component == 'U' or component == 'e':
557
            output[component] = numpy.array(output[component]).reshape(Ffield.shape)
558

559
560
    return output

561

562
def decomposePhiField(PhiField, components, twoD=False, nProcesses=nProcessesDefault, verbose=True):
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
    """
    This function takes in a Phi field (from readCorrelationTSV?) and
    returns fields of desired transformation components.

    Parameters
    ----------
        PhiField : multidimensional x 4 x 4 numpy array of floats
            Spatial field of Phis

        components : list of strings
            These indicate the desired components consistent with spam.deformation.decomposePhi

        twoD : bool, optional
            Is the PhiField in 2D? This changes the strain calculation.
            Default = False

579
580
581
582
583
584
585
586
        nProcesses : integer, optional
            Number of processes for multiprocessing
            Default = number of CPUs in the system

        verbose : boolean, optional
            Print progress?
            Default = True

587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
    Returns
    -------
        Dictionary containing appropriately reshaped fields of the transformation components requested.

        Keys:
            vol, dev, volss, devss are scalars
            t, r, and z are 3-component vectors
            e and U are 3x3 tensors
    """
    # The last two are for the 4x4 Phi field
    fieldDimensions  = PhiField.shape[0:-2]
    fieldRavelLength = numpy.prod(numpy.array(fieldDimensions))
    PhiFieldFlat     = PhiField.reshape(fieldRavelLength, 4, 4)

    output = {}
    for component in components:
603
604
605
606
        if component == 'vol' or component == 'dev' or component == 'volss' or component == 'devss':
            output[component] = numpy.zeros(fieldRavelLength)
        if component == 't' or component == 'r' or component == 'z':
            output[component] = numpy.zeros((fieldRavelLength, 3))
Olga Stamati's avatar
Olga Stamati committed
607
        if component == 'U' or component == 'e':
608
            output[component] = numpy.zeros((fieldRavelLength, 3, 3))
609

610
611
612
    # Function for parallel mode
    global decomposeOnePhi
    def decomposeOnePhi(n):
613
        Phi = PhiFieldFlat[n]
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
        if numpy.isfinite(Phi).sum() == 16:
            decomposedPhi = spam.deformation.decomposePhi(Phi, twoD=twoD)
            return n, decomposedPhi
        else:
            return n, {'t': numpy.array([numpy.nan] * 3),
                       'r': numpy.array([numpy.nan] * 3),
                       'z': numpy.array([numpy.nan] * 3),
                       'U': numpy.eye(3) * numpy.nan,
                       'e': numpy.eye(3) * numpy.nan,
                       'vol': numpy.nan,
                       'dev': numpy.nan,
                       'volss': numpy.nan,
                       'devss': numpy.nan}

    # Iterate through flat field of Fs
    if verbose: pbar = progressbar.ProgressBar(maxval=fieldRavelLength).start()

    finishedPoints = 0

    # Run multiprocessing
    with multiprocessing.Pool(processes=nProcesses) as pool:
        for returns in pool.imap_unordered(decomposeOnePhi, range(fieldRavelLength)):
            if verbose:
                finishedPoints += 1
                pbar.update(finishedPoints)
            for component in components:
                output[component][returns[0]] = returns[1][component]

    if verbose: pbar.finish()
643
644
645
646
647
648
649
650
651

    # Reshape on the output
    for component in components:
        if component == 'vol' or component == 'dev' or component == 'volss' or component == 'devss':
            output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2])

        if component == 't' or component == 'r' or component == 'z':
            output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3)

Olga Stamati's avatar
Olga Stamati committed
652
        if component == 'U' or component == 'e':
653
654
655
656
657
            output[component] = numpy.array(output[component]).reshape(*PhiField.shape[0:-2], 3, 3)

    return output


658
def correctPhiField(fileName=None, fieldCoords=None, fieldValues=None, fieldRS=None, fieldDPhi=None, fieldPSCC=None, fieldIT=None, fieldBinRatio=1.0, ignoreBadPoints=False, ignoreBackGround=False, correctBadPoints=False, deltaPhiNormMin=0.001, pixelSearchCCmin=0.98, nNeighbours=12, filterPoints=False, filterPointsRadius=3, verbose=False, saveFile=False, saveFileName=None):
659
660
661
    """
    This function corrects a field of deformation functions **Phi** calculated at a number of points.
    This is typically the output of the DICdiscrete and DICregularGrid clients.
Edward Andò's avatar
Edward Andò committed
662
    The correction is done based on the `returnStatus` and `deltaPhiNorm` of the correlated points.
663
664
665
666
667
668
669
    It takes as an input either a tsv file containing the result of the correlation or
    6 separate arrays:

        1 the coordinates of the points

        2 the PhiField

Edward Andò's avatar
Edward Andò committed
670
        3 the `returnStatus`
671

Edward Andò's avatar
Edward Andò committed
672
        4 the `deltaPhiNorm`
673
674
675

        5 the `PSCC`

Edward Andò's avatar
Edward Andò committed
676
        6 the `iterations`
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691


    Parameters
    ----------
        fileName : string, optional
            name of the file

        fieldCoords : 2D array, optional
            nx3 array of n points coordinates (ZYX)
            centre where each deformation function **Phi** has been calculated

        fieldValues : 3D array, optional
            nx4x4 array of n points deformation functions

        fieldRS : 1D array, optional
Edward Andò's avatar
Edward Andò committed
692
            nx1 array of n points `returnStatus` from the correlation
693
694

        fieldDf : 1D array, optional
Edward Andò's avatar
Edward Andò committed
695
            nx1 array of n points `deltaPhiNorm` from the correlation
696
697

        fieldIT : 1D array, optional
Edward Andò's avatar
Edward Andò committed
698
            nx1 array of n points `iterations` from the correlation
699
700

        fieldPSCC : 1D array, optional
701
            nx1 array of n points `PScorrelationCoefficient` from the correlation
702
703
704
705
706
707
708
709

        fieldBinRatio : int, optional
            if the input field is referred to a binned version of the image
            *e.g.*, if `fieldBinRatio = 2` the fileName values have been calculated for an image half the size of what the returned PhiField is referring to.
            Default = 1.0

        ignoreBadPoints : bool, optional
            if True it will replace the **Phi** matrices of the badly correlated points with nans.
Edward Andò's avatar
Edward Andò committed
710
            Bad points are set according to `returnStatus` and `deltaPhiNorm` or `PSCC` of the correlation.
711
712
713
714
            Default = False

         ignoreBackGround : bool, optional
            if True it will replace the **Phi** matrices of the back ground points with nans.
Edward Andò's avatar
Edward Andò committed
715
            Back ground points are set according to `returnStatus` (<-4) of the correlation.
716
717
718
719
            Default = False

        correctBadPoints : bool, optional
            if True it will replace the **Phi** matrices of the badly correlated points with the weighted function of the k nearest good points.
Edward Andò's avatar
Edward Andò committed
720
            Bad points are set according to `returnStatus` and `deltaPhiNorm` or `PSCC` of the correlation
721
            The number of the nearest good neighbours can be defined (see `nNeighbours` below).
722
723
724
            Default = False

        deltaPhiNormMin: float, optional
Edward Andò's avatar
Edward Andò committed
725
            minimum value of subpixel change in Phi to consider a point with `returnStatus` = 1 as good or bad.
726
727
728
729
730
731
            Default = 0.001

        picelSearchCCmin: float, optional
            minimum value of pixel search correlation coefficient to consider a point as good or bad.
            Default = 0.98

732
        nNeighbours : int, optional
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
            if `correctBadPoints` is activated, it specifies the number of the nearest neighbours to consider.
            If == 1, the nearest neighbour is used, if >1 neighbours are weighted according to distance.
            Default = 12

        filterPoints : bool, optional
            if True: a median filter will be applied on the **Phi** of each point.
            Default = False

        filterPointsRadius : int, optional
            Radius of median filter.
            Size of cubic structuring element is 2*filterPointsRadius+1.
            Default = 3

        verbose : bool, optional
            follow the progress of the function.
            Default = False

        saveFile : bool, optional
            save the corrected file into a tsv
            Default = False

        saveFileName : string, optional
            The file name for output.
            Default = 'spam'

    Returns
    -------
        PhiField : nx4x4 array
            n points deformation functions **Phi** after the correction

    """
    import os
    # read the input arguments
    if fileName:
        if not os.path.isfile(fileName):
            print("\n\tdeformationFunction.correctPhiField():{} is not a file. Exiting.".format(fileName))
            return
        else:
            import spam.helpers.tsvio
772
            fi = spam.helpers.tsvio.readCorrelationTSV(fileName, fieldBinRatio=fieldBinRatio, readConvergence=True, readPSCC=True)
773
774
775
            PhiField     = fi["PhiField"]
            fieldCoords  = fi["fieldCoords"]
            fieldDims    = fi["fieldDims"]
Edward Andò's avatar
Edward Andò committed
776
777
778
            RS           = fi["returnStatus"]
            deltaPhiNorm = fi["deltaPhiNorm"]
            iterations   = fi["iterations"]
779
            PSCC         = fi["PSCC"]
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
    elif fieldCoords is not None and fieldValues is not None and fieldRS is not None and fieldDPhi is not None and fieldPSCC is not None and fieldIT is not None:
        fieldCoords = fieldCoords
        fieldDims = numpy.array([len(numpy.unique(fieldCoords[:, 0])), len(numpy.unique(fieldCoords[:, 1])), len(numpy.unique(fieldCoords[:, 2]))])
        PhiField = fieldValues
        RS = fieldRS
        deltaPhiNorm = fieldDPhi
        PSCC = fieldPSCC
        iterations = fieldIT
    else:
        print("\tdeformationFunction.correctPhiField(): Not enough arguments given. Exiting.")
        return

    # check if it is a subPixel field or a pixel search field
    if numpy.nansum(PSCC) > 0 and iterations.sum() == 0:
        pixelSearch = True
        subPixel = False
    else:
        pixelSearch = False
        subPixel = True

Edward Andò's avatar
Edward Andò committed
800
    # define good and bad correlation points according to `returnStatus` and `deltaPhiNorm` or `PSCC`conditions
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
    if ignoreBackGround is False:
        if subPixel:
            goodPoints = numpy.where(numpy.logical_or(RS == 2, numpy.logical_and(RS == 1, deltaPhiNorm <= deltaPhiNormMin)))
            badPoints = numpy.where(numpy.logical_or(RS <= 0, numpy.logical_and(RS == 1, deltaPhiNorm > deltaPhiNormMin)))
        if pixelSearch:
            goodPoints = numpy.where(PSCC >= pixelSearchCCmin)
            badPoints = numpy.where(PSCC < pixelSearchCCmin)
    else:
        if subPixel:
            goodPoints = numpy.where(numpy.logical_or(RS == 2, numpy.logical_and(RS == 1, deltaPhiNorm <= deltaPhiNormMin)))
            badPoints = numpy.where(numpy.logical_or(numpy.logical_and(RS <= 0, RS >= -4), numpy.logical_and(RS == 1, deltaPhiNorm > deltaPhiNormMin)))
            backGroundPoints = numpy.where(RS < -4)
        if pixelSearch:
            goodPoints = numpy.where(numpy.logical_and(RS >= -4, PSCC >= pixelSearchCCmin))
            badPoints = numpy.where(numpy.logical_and(RS >= -4, PSCC < pixelSearchCCmin))
        backGroundPoints = numpy.where(RS < -4)
        PhiField[backGroundPoints] = numpy.nan

    # if asked, ignore the bad correlation points by setting their Phi to identity matrix
    if ignoreBadPoints:
        PhiField[badPoints] = numpy.eye(4) * numpy.nan

    # if asked, replace the bad correlation points with the weighted influence of the k nearest good neighbours
    if correctBadPoints:
        # create the k-d tree of the coordinates of good points, we need this to search for the k nearest neighbours easily
        #   for details see: https://en.wikipedia.org/wiki/K-d_tree &
        #   https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.query.html

        from scipy.spatial import KDTree
        treeCoord = KDTree(fieldCoords[goodPoints])

        # extract the Phi matrices of the bad points
        fieldBad = numpy.zeros_like(PhiField[badPoints])
        fieldBad[:, -1, :] = numpy.array([0, 0, 0, 1])

        # check if we have asked only for the closest neighbour
837
        if nNeighbours == 1:
838
839
840
841
842
843
844
845

            # loop over each bad point
            for badPoint in range(badPoints[0].shape[0]):
                if verbose:
                    print("\rWorking on bad point: {} of {}".format(badPoint + 1, badPoints[0].shape[0]), end='')
                # call tree.query to calculate:
                #   {ind}: the index of the nearest neighbour (as neighbours we consider only good points)
                #   {distnace}: distance (Minkowski norm 2, which is  the usual Euclidean distance) of the bad point to the nearest neighbour
846
                distance, ind = treeCoord.query(fieldCoords[badPoints][badPoint], k=nNeighbours)
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863

                # replace bad point's Phi with the Phi of the nearest good point
                fieldBad[badPoint][:-1] = PhiField[goodPoints][ind][:-1].copy()

            # replace the corrected Phi field
            PhiField[badPoints] = fieldBad

        # if we have asked for more neighbours
        else:

            # loop over each bad point
            for badPoint in range(badPoints[0].shape[0]):
                if verbose:
                    print("\rWorking on bad point: {} of {}".format(badPoint + 1, badPoints[0].shape[0]), end='')
                # call tree.query to calculate:
                #   {ind}: k nearest neighbours (as neighbours we consider only good points)
                #   {distnace}: distance (Minkowski norm 2, which is  the usual Euclidean distance) of the bad point to each of the ith nearest neighbour
864
                distance, ind = treeCoord.query(fieldCoords[badPoints][badPoint], k=nNeighbours)
865
866
867
868
869

                # compute the "Inverse Distance Weighting" since the nearest points should have the major influence
                weightSumInv = sum(1 / distance)

                # loop over each good neighbour point:
870
                for neighbour in range(nNeighbours):
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
                    # calculate its weight
                    weightInv = (1 / distance[neighbour]) / float(weightSumInv)

                    # replace the Phi components of the bad point with the weighted Phi components of the ith nearest good neighbour
                    fieldBad[badPoint][:-1] += PhiField[goodPoints][ind[neighbour]][:-1] * weightInv

            # replace the corrected Phi field
            PhiField[badPoints] = fieldBad
        # overwrite RS to the corrected
        RS[badPoints] = 2

    # if asked, apply a median filter of a specific size in the Phi field
    if filterPoints:
        if verbose:
            print("\nFiltering...")
        import scipy.ndimage
        filterPointsRadius = int(filterPointsRadius)

        PhiField[:, 0, 0] = scipy.ndimage.generic_filter(PhiField[:, 0, 0].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 1, 0] = scipy.ndimage.generic_filter(PhiField[:, 1, 0].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 2, 0] = scipy.ndimage.generic_filter(PhiField[:, 2, 0].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()

        PhiField[:, 0, 1] = scipy.ndimage.generic_filter(PhiField[:, 0, 1].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 1, 1] = scipy.ndimage.generic_filter(PhiField[:, 1, 1].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 2, 1] = scipy.ndimage.generic_filter(PhiField[:, 2, 1].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()

        PhiField[:, 0, 2] = scipy.ndimage.generic_filter(PhiField[:, 0, 2].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 1, 2] = scipy.ndimage.generic_filter(PhiField[:, 1, 2].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 2, 2] = scipy.ndimage.generic_filter(PhiField[:, 2, 2].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()

        PhiField[:, 0, -1] = scipy.ndimage.generic_filter(PhiField[:, 0, -1].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 1, -1] = scipy.ndimage.generic_filter(PhiField[:, 1, -1].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()
        PhiField[:, 2, -1] = scipy.ndimage.generic_filter(PhiField[:, 2, -1].reshape(fieldDims), numpy.nanmedian, size=(2 * filterPointsRadius + 1)).ravel()

        if ignoreBackGround:
            PhiField[backGroundPoints] = numpy.nan

    if saveFile:
        # if asked, write the corrected PhiField into a TSV
        if fileName:
            outDir = os.path.dirname(fileName)
            prefix = os.path.splitext(os.path.basename(fileName))[0]
            saveFileName = outDir+"/"+prefix
        elif saveFileName is None and fileName is None:
            saveFileName = "spam"

917
        TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\treturnStatus\tdeltaPhiNorm\titerations\tPSCC"
918
919
920
921
922
923
924
925
        outMatrix = numpy.array([numpy.array(range(PhiField.shape[0])),
                                 fieldCoords[:, 0], fieldCoords[:, 1], fieldCoords[:, 2],
                                 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],
                                 RS, deltaPhiNorm, iterations, PSCC]).T

        if filterPoints:
926
            title = "{}-corrected-N{}-filteredRad{}.tsv".format(saveFileName, nNeighbours, filterPointsRadius)
927
        else:
928
            title = "{}-corrected-N{}.tsv".format(saveFileName, nNeighbours)
929
930
931
932
933
934
935
936
937
938
939
        numpy.savetxt(title,
                      outMatrix,
                      fmt='%.7f',
                      delimiter='\t',
                      newline='\n',
                      comments='',
                      header=TSVheader)

    return PhiField


940
def mergeRegularGridAndDiscrete(regularGrid=None, discrete=None, labelledImage=None, binningLabelled=1, alwaysLabel=True, fileName=None):
941
    """
942
943
944
945
    This function corrects a Phi field from the spam-ldic script (measured on a regular grid)
    by looking into the results from one or more spam-ddic results (measured on individual labels)
    by applying the discrete measurements to the grid points.

946
947
    This can be useful where there are large flat zones in the image that cannot
    be correlated with small correlation windows, but can be identified and
948
    tracked with a spam-ddic computation (concrete is a good example).
949
950
951

    Parameters
    -----------
952
953
        regularGrid : dictionary
            Dictionary containing readCorrelationTSV of regular grid correlation script, `spam-ldic`.
954
955
            Default = None

956
957
        discrete : dictionary or list of dictionaries
            Dictionary (or list thereof) containing readCorrelationTSV of discrete correlation script, `spam-ddic`.
958
959
960
            File name of TSV from DICdiscrete client, or list of filenames
            Default = None

961
        labelledImage : 3D numpy array of ints, or list of numpy arrays
962
963
964
            Labelled volume used for discrete computation
            Default = None

Edward Andò's avatar
Edward Andò committed
965
966
967
968
969
        binningLabelled : int
            Are the labelled images and their PhiField at a different bin level than
            the regular field?
            Default = 1

970
971
972
        alwaysLabel : bool
            If regularGrid point falls inside the label, should we use the
            label displacement automatically?
973
            Otherwise if the regularGrid point has converged should we use that?
974
975
            Default = True (always use Label displacement)

976
977
        fileName : string
            Output filename, if None return dictionary as from spam.helpers.readCorrelationTSV()
978
979
980
981
            Default = None

    Returns
    --------
982
983
        Either dictionary or TSV file
            Output matrix, with number of rows equal to spam-ldic (the node spacing of the regular grid) and with columns:
Edward Andò's avatar
Edward Andò committed
984
            "NodeNumber", "Zpos", "Ypos", "Xpos", "Zdisp", "Ydisp", "Xdisp", "deltaPhiNorm", "returnStatus", "iterations"
985
986
987
988
    """
    import spam.helpers

    # If we have a list of input discrete files, we also need a list of labelled images
989
990
    if type(discrete) == list:
        if type(labelledImage) != list:
991
            print("spam.deformation.deformationFunction.mergeRegularGridAndDiscrete(): if you pass a list of discreteTSV you must also pass a list of labelled images")
992
            return
993
994
        if len(discrete) != len(labelledImage):
            print("spam.deformation.deformationFunction.mergeRegularGridAndDiscrete(): len(discrete) must be equal to len(labelledImage)")
995
            return
996
        nDiscrete = len(discrete)
997
998
999

    # We have only one TSV and labelled image, it should be a number array
    else:
1000
1001
        if type(labelledImage) != numpy.ndarray:
            print("spam.deformation.deformationFunction.mergeRegularGridAndDiscrete(): with a single discrete TSV file, labelledImage must be a numpy array")
1002
1003
            return
        discrete = [discrete]
1004
        labelledImage = [labelledImage]
1005
1006
        nDiscrete = 1

1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
    output = {}

    # Regular grid is the master, and so we copy dimensions and positions
    output['fieldDims']   = regularGrid['fieldDims']
    output['fieldCoords'] = regularGrid['fieldCoords']

    output['PhiField']     = numpy.zeros_like(regularGrid['PhiField'])
    output['iterations']   = numpy.zeros_like(regularGrid['iterations'])
    output['deltaPhiNorm'] = numpy.zeros_like(regularGrid['deltaPhiNorm'])
    output['returnStatus'] = numpy.zeros_like(regularGrid['returnStatus'])
    output['mergeSource']  = numpy.zeros_like(regularGrid['iterations'])
1018
1019
1020
1021
1022

    # from progressbar import ProgressBar
    # pbar = ProgressBar()

    # For each point on the regular grid...
1023
1024
1025
    #for n, gridPoint in pbar(enumerate(regularGrid['fieldCoords'].astype(int))):
    for n, gridPoint in enumerate(regularGrid['fieldCoords'].astype(int)):
        # Find labels corresponding to this grid position for the labelledImage images
1026
1027
        labels = []
        for m in range(nDiscrete):
1028
1029
1030
            labels.append(int(labelledImage[m][int(gridPoint[0] / float(binningLabelled)),
                                               int(gridPoint[1] / float(binningLabelled)),
                                               int(gridPoint[2] / float(binningLabelled))]))
1031
1032
1033
        labels = numpy.array(labels)

        # Is the point inside a discrete label?
1034
1035
        if (labels == 0).all() or (not alwaysLabel and regularGrid['returnStatus'][n] == 2):
            ### Use the REGULAR GRID MEASUREMENT
1036
            # If we're not in a label, copy the results from DICregularGrid
1037
1038
1039
1040
            output['PhiField'][n]     = regularGrid['PhiField'][n]
            output['deltaPhiNorm'][n] = regularGrid['deltaPhiNorm'][n]
            output['returnStatus'][n] = regularGrid['returnStatus'][n]
            output['iterations'][n]   = regularGrid['iterations'][n]
1041
        else:
1042
            ### Use the DISCRETE MEASUREMENT
1043
1044
1045
1046
            # Give precedence to earliest non-zero-labelled discrete field, conflicts not handled
            m = numpy.where(labels != 0)[0][0]
            label = labels[m]
            #print("m,label = ", m, label)
1047
1048
1049
1050
            tmp = discrete[m]['PhiField'][label].copy()
            tmp[0:3, -1] *= float(binningLabelled)
            translation = spam.deformation.decomposePhi(tmp, PhiCentre=discrete[m]['fieldCoords'][label] * float(binningLabelled), PhiPoint=gridPoint)['t']
            # This is the Phi we will save for this point -- take the F part of the labelled's Phi
1051
            phi = discrete[m]['PhiField'][label].copy()
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
            # ...and add the computed displacement as applied to the grid point
            phi[0:3, -1] = translation

            output['PhiField'][n]     = phi
            output['deltaPhiNorm'][n] = discrete[m]['deltaPhiNorm'][label]
            output['returnStatus'][n] = discrete[m]['returnStatus'][label]
            output['iterations'][n]   = discrete[m]['iterations'][label]
            output['mergeSource'][n]  = m+1

    if fileName is not None:
        TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp"
        outMatrix = numpy.array([numpy.array(range(output['fieldCoords'].shape[0])),
                                 output['fieldCoords'][:, 0],     output['fieldCoords'][:, 1],    output['fieldCoords'][:, 2],
                                 output['PhiField'][:, 0, 0],     output['PhiField'][:, 0, 1],    output['PhiField'][:, 0, 2],    output['PhiField'][:, 0, 3],
                                 output['PhiField'][:, 1, 0],     output['PhiField'][:, 1, 1],    output['PhiField'][:, 1, 2],    output['PhiField'][:, 1, 3],
                                 output['PhiField'][:, 2, 0],     output['PhiField'][:, 2, 1],    output['PhiField'][:, 2, 2],    output['PhiField'][:, 2, 3]]).T

        outMatrix = numpy.hstack([outMatrix, numpy.array([output['iterations'],
                                                          output['returnStatus'],
                                                          output['deltaPhiNorm'],
                                                          output['mergeSource']]).T])
        TSVheader = TSVheader+"\titerations\treturnStatus\tdeltaPhiNorm\tmergeSource"

        numpy.savetxt(fileName,
                      outMatrix,
                      fmt='%.7f',
                      delimiter='\t',
                      newline='\n',
                      comments='',
                      header=TSVheader)
    else:
        return output
1084
1085


1086
def interpolatePhiField(fieldCoords, fieldValues, interpCoords, fieldInterpBinRatio=1.0):
1087
1088
1089
1090
1091
    """
    Interpolate a field of deformation functions (Phi).

    Parameters
    ----------
1092
1093
        fieldCoords : nPointsField x 3 array
            Z Y X coordinates of points where ``fieldValues`` are defined
1094

1095
1096
        fieldValues : nPointsField x 4 x 4 array
            Phi defined at ``fieldCoords``
1097

1098
1099
        interpCoords : nPointsInterpolate x 3
            Z Y X coordinates of points to interpolate Phi for
1100

1101
1102
1103
1104
        fieldInterpBinRatio : int, optional
            If the ``fieldCoords`` and ``fieldValues`` matrices refer to a binned version of the new coordintes.
            `e.g.`, if ``fieldInterpBinRatio = 2`` then ``fieldCoords`` and ``fieldValues`` have been calculated on an
            image half the size of what ``interpCoords`` are referring to.
1105
1106
1107

    Returns
    -------
1108
        interpValues : nPointsInterpolate x 4 x 4 array of Phis
1109
    """
1110
1111
1112
1113
    import scipy.ndimage

    assert(type(interpCoords)==numpy.ndarray), "spam.deformation.deformationField.interpolatePhiField(): Coordinates to interpolate should be numpy array"
    assert(interpCoords.dtype==float), "spam.deformation.deformationField.interpolatePhiField(): Coordinates to interpolate should be float"
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
    # This version of the function will use scipy.ndimage.interpolation.map_coordinates().
    # It takes in a field, which means that our fieldValues Phi field MUST be regularly spaced.
    # Furthermore it takes points at integer values (voxels), so we have to convert from
    # positions in the Phi field, and the "real" voxel coordinates.
    # e.g., Our first measurement point is 12,12,12 and the node spacing is 20 pixels.
    # map_coordinates will access this first Phi 12,12,12 at a position [0,0,0] in the matrix of Phi values in space
    # The next Phi 32,12,12 at a position [1,0,0]
    # Define the output array
    output = numpy.zeros((interpCoords.shape[0], 4, 4))

    # 1. calculate node spacing and position of first point
    # Measure node spacing in all three directions:
    zUnique = numpy.unique(fieldCoords[:, 0])
    yUnique = numpy.unique(fieldCoords[:, 1])
    xUnique = numpy.unique(fieldCoords[:, 2])

    zSpacing = zUnique[1] - zUnique[0]
    ySpacing = yUnique[1] - yUnique[0]
    xSpacing = xUnique[1] - xUnique[0]

    if zSpacing == ySpacing and zSpacing == xSpacing:
        nodeSpacing = zSpacing

        # TopPoint -- Ask ER -- Elizabeth Regina -- Honni soit qui mal y pense
        taupPoihunt = [zUnique[0], yUnique[0], xUnique[0]]
        # print "Top point:", taupPoihunt

1141
        dims = [  int(1 + (zUnique[-1] - zUnique[0]) / zSpacing),
1142
1143
1144
1145
1146
1147
1148
1149
                  int(1 + (yUnique[-1] - yUnique[0]) / ySpacing),
                  int(1 + (xUnique[-1] - xUnique[0]) / xSpacing)]
    else:
        print("Not all node spacings are the same, help! {} {} {} ".format(
            zSpacing, ySpacing, xSpacing))
        return "moinsUn"

    # 2. reshape fieldValues into an Z*Y*X*Phiy*Phix array for map_coordinates
1150
    fieldValues = fieldValues.reshape([dims[0], dims[1], dims[2], 4, 4])
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164

    # 3. Convert interpCoords into positions in reshaped Phi array
    # If we have a non-zero bin, scale coordinates
    interpCoords /= fieldInterpBinRatio

    # Remove top corner coords...
    interpCoords -= taupPoihunt
    # And divide by node spacing, now coords are in 0->1 format
    interpCoords /= float(nodeSpacing)

    # 4. Call map_coordinates and return
    # Loop over each component of Phi, so they are not interpolated together.
    for Fy in range(4):
        for Fx in range(4):
1165
            output[:, Fy, Fx] = scipy.ndimage.interpolation.map_coordinates(fieldValues[:, :, :, Fy, Fx], interpCoords.T, order=1)
1166
1167
1168
1169

    # 5. Scale transformation by binning value
    output[:, 0:3, 3] *= fieldInterpBinRatio
    return output
1170

1171
1172
1173
def getDisplacementFromNeighbours(labIm, DVC, fileName, method = 'getLabel', centresOfMass = None, previousDVC = None):
    """
    This function computes the displacement as the mean displacement from the neighbours, for non-converged grains using a TSV file obtained from `spam-ddic` script. Returns a new TSV file with the new Phi (composed only by the displacement part). 
1174

1175
1176
1177
1178
1179
1180
    The generated TSV can be used as an input for `spam-ddic`. 

    Parameters
    -----------
        lab : 3D array of integers
            Labelled volume, with lab.max() labels
1181

1182
        DVC : dictionary
1183
            Dictionary with deformation field, obtained from `spam-ddic` script, and read using `spam.helpers.tsvio.readCorrelationTSV()` with `readConvergence=True, readPSCC=True, readLabelDilate=True`
1184

1185
1186
        fileName : string
            FileName including full path and .tsv at the end to write
1187

1188
1189
1190
1191
1192
        method : string, optional
            Method to compute the neighbours using `spam.label.getNeighbours()`.
            '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'
1193

1194
1195
1196
        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``
1197

1198
1199
1200
1201
1202
1203
1204
1205
1206
        previousDVC : dictionary, optional
            Dictionary with deformation field, obtained from `spam-ddic` script, and read using `spam.helpers.tsvio.readCorrelationTSV()` for the previous step.
            This allows the to compute only the displacement increment from the neighbours, while using the F tensor from a previous (converged) step.
            If `previousDVS = None`, then the resulting Phi would be composed only by the displacement of the neighbours.
            Default = None

    Returns
    --------
        Dictionary
1207
            TSV file with the same columns as the input
1208
    """
1209

1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
    # Compute centreOfMass if needed
    if centresOfMass == None:
        centresOfMass = spam.label.centresOfMass(labIm)
    # Get number of labels
    numberOfLabels = (labIm.max() + 1).astype('u4')
    # Create Phi field
    PhiField = numpy.zeros((numberOfLabels, 4, 4), dtype='<f4')
    # Rest of arrays
    try:
        iterations = DVC['iterations']
        returnStatus = DVC['returnStatus']
        deltaPhiNorm = DVC['deltaPhiNorm']
        PSCC = DVC['PSCC']
        labelDilateList = DVC['LabelDilate']
        error = DVC['error']
1225

1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
        # Get the problematic labels
        probLab = numpy.where(DVC['returnStatus']!= 2)[0]
        # Remove the 0 from the wrongLab list
        probLab = probLab[~numpy.in1d(probLab, 0)]
        # Get neighbours
        neighbours = spam.label.getNeighbours(labIm, probLab, method = method)
        # Solve first the converged particles - make a copy of the PhiField
        for i in range(numberOfLabels):
            PhiField[i] = DVC['PhiField'][i]
        # Work on the problematic labels
        widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
        pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(probLab))
        pbar.start()
        for i in range(0, len(probLab), 1):
            wrongLab = probLab[i]
            neighboursLabel = neighbours[i]
            t = []
            for n in neighboursLabel:
                # Check the return status of each neighbour 
                if DVC['returnStatus'][n] == 2:
                    # We dont have a previous DVC file loaded
                    if previousDVC is None:
                        # Get translation, rotation and zoom from Phi
                        nPhi = spam.deformation.deformationFunction.decomposePhi(DVC['PhiField'][n])
                        # Append the results
                        t.append(nPhi['t'])
                    # We have a previous DVC file loaded
                    else:
                        # Get translation, rotation and zoom from Phi at t=step
                        nPhi = spam.deformation.deformationFunction.decomposePhi(DVC['PhiField'][n])
                        # Get translation, rotation and zoom from Phi at t=step-1
                        nPhiP = spam.deformation.deformationFunction.decomposePhi(previousDVC['PhiField'][n])
                        # Append the incremental results
                        t.append(nPhi['t']-nPhiP['t'])
            # Transform list to array
            if not t:
                # This is a non-working label, take care of it
                Phi = spam.deformation.computePhi({'t': [0,0,0]})
                PhiField[wrongLab] = Phi
            else:
                t = numpy.asarray(t)
                # Compute mean
                meanT = numpy.mean(t, axis = 0)
                # Reconstruct 
                transformation = {'t': meanT}
                Phi = spam.deformation.computePhi(transformation)
                # Save
                if previousDVC is None:
                    PhiField[wrongLab] = Phi
                else:
                    PhiField[wrongLab] = previousDVC['PhiField'][wrongLab]
                    # Add the incremental displacement
                    PhiField[wrongLab][:-1,-1] += Phi[:-1,-1]
1279
1280
1281
            # Update the progressbar
            widgets[0] = progressbar.FormatLabel("{}/{} ".format(i, numberOfLabels))
            pbar.update(i)
1282
1283
1284
1285
1286
1287
1288
1289
1290
        # Save
        outMatrix = numpy.array([numpy.array(range(numberOfLabels)), 
                                centresOfMass[:, 0], centresOfMass[:, 1], centresOfMass[:, 2],
                                PhiField[:, 0, 3], PhiField[:, 1, 3], PhiField[:, 2, 3], 
                                PhiField[:, 0, 0], PhiField[:, 0, 1], PhiField[:, 0, 2], 
                                PhiField[:, 1, 0], PhiField[:, 1, 1], PhiField[:, 1, 2], 
                                PhiField[:, 2, 0], PhiField[:, 2, 1], PhiField[:, 2, 2], 
                                PSCC, error, iterations, returnStatus, deltaPhiNorm, labelDilateList]).T
        numpy.savetxt(fileName, 
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
                      outMatrix, 
                      fmt='%.7f', 
                      delimiter='\t',
                      newline='\n', 
                      comments='', 
                      header="Label\tZpos\tYpos\tXpos\t" + 
                              "Zdisp\tYdisp\tXdisp\t" + 
                              "Fzz\tFzy\tFzx\t" + 
                              "Fyz\tFyy\tFyx\t" + 
                              "Fxz\tFxy\tFxx\t" + 
                              "PSCC\terror\titerations\treturnStatus\tdeltaPhiNorm\tLabelDilate") 
1302
1303
1304
    except:
        print('spam.deformation.deformationField.getDisplacementFromNeighbours(): Missing information in the input TSV. Make sure you are reading iterations, returnStatus, deltaPhiNorm, PSCC, LabelDilate, and error.')
        print('spam.deformation.deformationField.getDisplacementFromNeighbours(): Aborting') 
1305

1306
def mergeRegistrationAndDiscreteFields(regTSV, discreteTSV, fileName, regTSVBinRatio = 1):
1307
1308
1309
    """
    This function merges a registration TSV with a discrete TSV.
    Can be used to create the first guess for `spam-ddic`, using the registration over the whole file, and a previous result from `spam-ddic`.
1310

1311

1312
1313
    Parameters
    -----------
1314

1315
1316
        regTSV : dictionary
            Dictionary with deformation field, obtained from a registrion, usually from the whole sample, and read using `spam.helpers.tsvio.readCorrelationTSV()`
1317

1318
1319
        discreteTSV : dictionary
            Dictionary with deformation field, obtained from `spam-ddic` script, and read using `spam.helpers.tsvio.readCorrelationTSV()`
1320

1321
1322
        fileName : string
            FileName including full path and .tsv at the end to write
1323

1324
1325
        regTSVBinRatio : float, optional
            Change translations from regTSV, if it's been calculated on a differently-binned image. Default = 1
1326

1327
1328
1329
1330
1331
    Returns
    --------
        Dictionary
            TSV file with the same columns as the input
    """
1332

1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
    # Create a first guess
    phiGuess = discreteTSV['PhiField'].copy()
    # Main loop
    for lab in range(discreteTSV['numberOfLabels']):
        # Initial position of a grain
        iniPos = discreteTSV['fieldCoords'][lab]
        # Position of the label at T+1
        deformPos = iniPos + discreteTSV['PhiField'][lab][:-1,-1]
        # Compute the extra displacement and rotation
        extraDisp = spam.deformation.decomposePhi(regTSV['PhiField'][0], 
                                              PhiCentre = regTSV['fieldCoords'][0], 
                                              PhiPoint = deformPos)['t']
        # Add the extra disp to the phi guess
1346
        phiGuess[lab][:-1,-1] += extraDisp*regTSVBinRatio
1347

1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
    # Save
    outMatrix = numpy.array([numpy.array(range(discreteTSV['numberOfLabels'])),
                            discreteTSV['fieldCoords'][:, 0], 
                            discreteTSV['fieldCoords'][:, 1], 
                            discreteTSV['fieldCoords'][:, 2], 
                            phiGuess[:, 0, 3], phiGuess[:, 1, 3], phiGuess[:, 2, 3], 
                            phiGuess[:, 0, 0], phiGuess[:, 0, 1], phiGuess[:, 0, 2], 
                            phiGuess[:, 1, 0], phiGuess[:, 1, 1], phiGuess[:, 1, 2], 
                            phiGuess[:, 2, 0], phiGuess[:, 2, 1], phiGuess[:, 2, 2], 
                            numpy.zeros(((discreteTSV['numberOfLabels'])), dtype='<f4'), 
                            discreteTSV['iterations'], 
                            discreteTSV['returnStatus'], 
                            discreteTSV['deltaPhiNorm']]).T
    numpy.savetxt(fileName, 
                  outMatrix, 
                  fmt='%.7f', 
                  delimiter='\t', 
                  newline='\n', 
                  comments='', 
                  header="Label\tZpos\tYpos\tXpos\t" + 
                         "Zdisp\tYdisp\tXdisp\t" + 
                         "Fzz\tFzy\tFzx\t" + 
                         "Fyz\tFyy\tFyx\t" + 
                         "Fxz\tFxy\tFxx\t" + "PSCC\titerations\treturnStatus\tdeltaPhiNorm")