grid.py 10.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
"""
Library of SPAM functions for defining a regular grid in a reproducible way.
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/>.
"""

from __future__ import print_function
20

21
import numpy
Edward Andò's avatar
Edward Andò committed
22
import spam.DIC.correlate
Edward Andò's avatar
Edward Andò committed
23
import sys
Emmanuel Roubin's avatar
Emmanuel Roubin committed
24
import progressbar
Edward Andò's avatar
Edward Andò committed
25
26
27
28
if sys.version[0] == '2':
    import Queue as queue
else:
    import queue as queue
29

30

31
def makeGrid(imageSize, nodeSpacing):
32
33
34
35
36
37
38
39
40
41
42
43
44
    """
        Define a grid of correlation points.

        Parameters
        ----------
        imageSize : 3-item list
            Size of volume to spread the grid inside

        nodeSpacing : 3-item list or int
            Spacing between nodes

        Returns
        -------
45
        nodePositions : nPointsx3 numpy.array
46
47
            Array containing Z, Y, X positions of each point in the grid
    """
48
49
50

    if len(imageSize) != 3:
        print("\tgrid.makeGrid(): imageSize doesn't have three dimensions, exiting")
51
52
        return

53
    if type(nodeSpacing) == int or type(nodeSpacing) == float:
Emmanuel Roubin's avatar
Emmanuel Roubin committed
54
55
        nodeSpacing = [nodeSpacing] * 3
    elif len(nodeSpacing) != 3:
56
        print("\tgrid.makeGrid(): nodeSpacing is not an int or float and doesn't have three dimensions, exiting")
57
58
        return

Emmanuel Roubin's avatar
Emmanuel Roubin committed
59
60
61
62
    if imageSize[0] == 1:
        twoD = True
    else:
        twoD = False
63

64
65
    # Note: in this cheap node spacing, the first node is always at a distance of --nodeSpacing-- from the origin
    # The following could just be done once in principle...
66
67
68
    nodesMgrid = numpy.mgrid[nodeSpacing[0]:imageSize[0]:nodeSpacing[0],
                             nodeSpacing[1]:imageSize[1]:nodeSpacing[1],
                             nodeSpacing[2]:imageSize[2]:nodeSpacing[2]]
69
70
71

    # If twoD then overwrite nodesMgrid
    if twoD:
Emmanuel Roubin's avatar
Emmanuel Roubin committed
72
        nodesMgrid = numpy.mgrid[0: 1: 1,
73
74
                                 nodeSpacing[1]:imageSize[1]:nodeSpacing[1],
                                 nodeSpacing[2]:imageSize[2]:nodeSpacing[2]]
75

Emmanuel Roubin's avatar
Emmanuel Roubin committed
76
    nodesDim = (nodesMgrid.shape[1], nodesMgrid.shape[2], nodesMgrid.shape[3])
77

78
    numberOfNodes = int(nodesMgrid.shape[1] * nodesMgrid.shape[2] * nodesMgrid.shape[3])
79

80
    nodePositions = numpy.zeros((numberOfNodes, 3))
81

Emmanuel Roubin's avatar
Emmanuel Roubin committed
82
83
84
    nodePositions[:, 0] = nodesMgrid[0].ravel()
    nodePositions[:, 1] = nodesMgrid[1].ravel()
    nodePositions[:, 2] = nodesMgrid[2].ravel()
85

Edward Andò's avatar
Edward Andò committed
86
    return nodePositions, nodesDim
87

88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
"""
This function handles grid-based local correlation, offering an initial rough dispalcement-only guess.
At the moment matching of windows is done with a Normalised-Correlation-Coefficient approach.

Parameters
----------
    im1 : 3D numpy array
        A 3D image of greylevels defining a reference configuration for the pixel search

    im2 : 3D numpy array
        A deformed 3D image of greylevels

    nodePositions : nPoints*3 numpy array
        Array containing Z, Y, X positions of each point in the grid, as returned by ``makeGrid`` for example
        defined at im1

    halfWindowSize : 3-item list or int
        Size of subvolumes to perform the image correlation on, as a data range taken either side of the voxel on which the node is placed.
        The subvolume will be 2*halfWindowSize + 1 pixels on each side.
        A general recommendation is to make this half the node spacing

    searchRange : dictionary
        Search range as a dictionary containing 3 keys: 'zRange', 'yRange', and 'xRange',
        Each of which contains a list with two items

    PhiField : nPoints*4*4 numpy array, optional
        Optional field of ``F`` transformation operators defined for each node.
        Currently, only the translational components of F will be taken into account.
        Default = No displacement

    minMaskCoverage : float, optional
        Minimum number of pixels in a subvolume for it to be correlated (only considered in the case of im1mask).
        Default = 0.5

    im1mask : 3D boolean numpy array, optional
        A mask for im1 which is true in the zones to correlate.
        Default = None

    greyThreshold : list of two floats, optional
        Threshold for the mean greylevel in each im1 subvolume.
        If the mean is below the first value or above the second value, the grid point is not correlated.
        Default = [ -inf, inf ]

131
132
133
    twoD : bool, optional
        Are passed imagette1 and imagette2 3D images?

134
135
136
137
138
139
140
141
142
143
144
Returns
-------
    Dictionary containing:

    Keys
        PhiField : nNodes*4*4 numpy array of floats
            For each node, the measured transformatio operator (displacement only)

        pixelSearchCC : nNodes numpy array of floats
            For each node, the NCC score obtained
"""
145

146
# WARNING ------------------------- VVVVVVVVVVV--- gets easily overwritten, pass a .copy()!
147
def getImagettes(nodePosition, Phi, searchRange, boundingBox, im1, im2, im1mask, im2mask, minMaskCoverage, greyThreshold, twoD=False):
148
149
    returnStatus = 1
    imagette1mask = None
150
    imagette2mask = None
151
    initialDisplacement = numpy.round(Phi[0:3, 3]).astype(int)
152

153
    # Catch bad bounding boxes:
154
    if numpy.all((boundingBox[1::2] - boundingBox[0::2]) > 0) or ( numpy.all(boundingBox[3::2] - boundingBox[2::2] > 0) and twoD):
155
156
        # 2020-09-25 OS and EA: Prepare startStop array for imagette 1 to be extracted with new slicePadded
        PhiNoDisp = Phi.copy()
157
158
        PhiNoDisp[0:3,-1] -= initialDisplacement
        #PhiNoDisp[0:3,-1] = 0.0
159

160
161
        # If F is not the identity, create a pad to be able to apply F to imagette 1
        if numpy.allclose(PhiNoDisp, numpy.eye(4)):
162
            applyPhiPad = (0, 0, 0)
163
164
165
        else:
            # 2020-10-06 OS and EA: Add a pad to each dimension of 25% of max(halfWindowSize) to allow space to apply F (no displacement) to imagette1
            applyPhiPad = int(0.5*numpy.ceil(max(boundingBox[1::2]-boundingBox[0::2])))
166

167
168
169
170
171
172
            if twoD: applyPhiPad = (          0, applyPhiPad, applyPhiPad)
            else:    applyPhiPad = (applyPhiPad, applyPhiPad, applyPhiPad)

        startStopIm1 = [int(boundingBox[0] - applyPhiPad[0]), int(boundingBox[1] + applyPhiPad[0] + 1),
                        int(boundingBox[2] - applyPhiPad[1]), int(boundingBox[3] + applyPhiPad[1] + 1),
                        int(boundingBox[4] - applyPhiPad[2]), int(boundingBox[5] + applyPhiPad[2] + 1)]
173
174
175
176
177
178
179
180
181
182

        # In either case, extract imagette1, now guaranteed to be the right size
        imagette1padded = spam.helpers.slicePadded(im1, startStopIm1)

        # If F is not the identity, apply F and undo crop
        if numpy.allclose(PhiNoDisp, numpy.eye(4)):
            # In this case there is is no padding (despite the name) and we can just keep going
            imagette1def = imagette1padded
        else:
            # apply F to imagette 1 padded
183
184
185
            if twoD:    imagette1paddedDef = spam.DIC.applyPhiPython(imagette1padded, PhiNoDisp)
            else:       imagette1paddedDef = spam.DIC.applyPhi(imagette1padded, PhiNoDisp)

186
            # undo padding
187
188
189
190
191
192
193
194
            if twoD:
                imagette1def = imagette1paddedDef[              :               ,
                                                  applyPhiPad[1]:-applyPhiPad[1],
                                                  applyPhiPad[2]:-applyPhiPad[2]]
            else:
                imagette1def = imagette1paddedDef[applyPhiPad[0]:-applyPhiPad[0],
                                                  applyPhiPad[1]:-applyPhiPad[1],
                                                  applyPhiPad[2]:-applyPhiPad[2]]
195
196
197
198
199
200
201
202
203
204
205
206
207

        ### Check mask
        if im1mask is None:
            # no mask1 --> always pas this test (e.g., labelled image)
            maskVolumeCondition = True
            imagette1mask = None
        else:
            imagette1mask = spam.helpers.slicePadded(im1mask, boundingBox + numpy.array([0, 1, 0, 1, 0, 1])) != 0
            maskVolumeCondition = imagette1mask.mean() >= minMaskCoverage


        # Make sure imagette is not 0-dimensional in any dimension
        # Check minMaskVolume
208
        if numpy.all(numpy.array(imagette1def.shape) > 0) or (twoD and numpy.all(numpy.array(imagette1def.shape[1:3]) > 0)):
209
210
211
212
213
214
215
216
217
218
219
220
221
222
            # ------------ Grey threshold low --------------- and -------------- Grey threshold high -----------
            if numpy.nanmean(imagette1def) > greyThreshold[0] and numpy.nanmean(imagette1def) < greyThreshold[1]:
                if maskVolumeCondition:
                    # Slice for image 2
                    ## 2020-09-25 OS and EA: Prepare startStop array for imagette 2 to be extracted with new slicePadded
                    ## Extract it...
                    startStopIm2 = [int(boundingBox[0] + initialDisplacement[0] + searchRange[0]    ),
                                    int(boundingBox[1] + initialDisplacement[0] + searchRange[1] + 1),
                                    int(boundingBox[2] + initialDisplacement[1] + searchRange[2]    ),
                                    int(boundingBox[3] + initialDisplacement[1] + searchRange[3] + 1),
                                    int(boundingBox[4] + initialDisplacement[2] + searchRange[4]    ),
                                    int(boundingBox[5] + initialDisplacement[2] + searchRange[5] + 1)]
                    imagette2 = spam.helpers.slicePadded(im2, startStopIm2)

223
224
                    if im2mask is not None: imagette2mask = spam.helpers.slicePadded(im2mask, startStopIm2)
                        
225
                # Failed minMaskVolume condition
226
                else:
227
                    returnStatus = -5
228
                    imagette1def = None
229
                    imagette2 = None
Edward Andò's avatar
Edward Andò committed
230

231
            # Failed greylevel condition
232
233
            else:
                returnStatus = -5
234
                imagette1def = None
235
                imagette2 = None
236

237
        # Failed 0-dimensional imagette test
238
        else:
239
240
            returnStatus = -5
            imagette1def = None
241
            imagette2 = None
242

243
    # Failed bounding box test
244
    else:
245
246
247
248
249
250
251
        returnStatus = -7
        imagette1def = None
        imagette2 = None

    return {'imagette1':            imagette1def,
            'imagette1mask':        imagette1mask,
            'imagette2':            imagette2,
252
            'imagette2mask':        imagette2mask,
253
254
            'returnStatus':         returnStatus,
            'pixelSearchOffset':    searchRange[0::2]
255
            }