myfirstorder.py 9.95 KB
Newer Older
1
2
import numpy
from scipy import ndimage
3
from oct2py import Oct2Py, Oct2PyError
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

from radiomics.firstorder import RadiomicsFirstOrder
from .mybase import MyRadiomicsFeaturesBase


class MyRadiomicsFirstOrder(MyRadiomicsFeaturesBase):
    def __init__(self, inputImage, inputMask, **kwargs):
        super().__init__(inputImage, inputMask, **kwargs)

        self.pixelSpacing = inputImage.GetSpacing()
        self.voxelArrayShift = kwargs.get('voxelArrayShift', 0)

    def _initCalculation(self):
        return super()._initSegmentBasedCalculation()

        self.targetVoxelArray = self.imageArray[self.labelledVoxelCoordinates].astype(
            'float'
        )
        self.discretizedTargetVoxelArray = None  # Lazy instantiation

        self.logger.debug('First order feature class initialized')

    def _moment(self, a, moment=1, axis=0):
        return RadiomicsFirstOrder._moment(self, a, moment, axis)

    def _getDiscretizedTargetVoxelArray(self):
        return RadiomicsFirstOrder._getDiscretizedTargetVoxelArray(self)

    def getSuvMaxFeatureValue(self):
        r"""
        **20. SUV max**

        Implemented as one of the intensity features extracted by Vallieres et al.
        Maximum SUV of the tumour region. 

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        return numpy.max(ROIPet[~numpy.isnan(ROIPet)])

    def getSuvMeanFeatureValue(self):
        r"""
        **21. SUV mean**

        Implemented as one of the intensity features extracted by valieres et al.
        Average SUV of the tumour region. 

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        return numpy.mean(ROIPet[~numpy.isnan(ROIPet)])

64
65
    def getTLGFeatureValue(self):
        r"""
66
        **22. TLG**
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

        Total lesion glycolysis. 
        Defined as 
        
        .. math::
        \textit{SUVMean}\times \textit{total volume of the tumour region}
    
        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        z, y, x = self.pixelSpacing
        Np = len(self.labelledVoxelCoordinates[0])
        volume = Np * (z * x * y)

        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        return numpy.mean(ROIPet[~numpy.isnan(ROIPet)]) * volume

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
    def getAUCCSHFeatureValue(self):
        r"""
        **23. Inactive volume **

        Area under the curve of the cumulative SUV-volume histogram describing the percentage 
        of total tumour volume above a percentage threshold of maximum SUV, as defined by van Velden et al.

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        nBins = 1000  # By default
        ROIPet = self.imageArray
        mask = self.maskArray

        ROIPet[~mask] = numpy.nan

        outliers = numpy.where(
            ROIPet
            > (
                numpy.mean(ROIPet[~numpy.isnan(ROIPet)])
                + 3 * numpy.std(ROIPet[~numpy.isnan(ROIPet)])
            )
        )[0]

        good_voxels = numpy.where(
            ROIPet
            <= (
                numpy.mean(ROIPet[~numpy.isnan(ROIPet)])
                + 3 * numpy.std(ROIPet[~numpy.isnan(ROIPet)])
            )
        )[0]

        ROIPet[outliers] = numpy.mean(ROIPet[good_voxels])
        ROIPet = ROIPet - numpy.min(ROIPet[~numpy.isnan(ROIPet)])
        ROIPet = ROIPet / numpy.max(ROIPet[~numpy.isnan(ROIPet)])
        volume = ROIPet[~numpy.isnan(ROIPet)]
        nVoxel = len(volume)

        bins, _ = numpy.histogram(volume, nBins)

        csh = numpy.flipud(numpy.cumsum(numpy.flipud(bins)) / nVoxel)
        return numpy.sum(csh / nBins)

Alessia Marcolini's avatar
Alessia Marcolini committed
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
    def getgETUFeatureValue(self):
        r"""
        **24. gETU**

        Generalized effective total uptake, with parameter a = 0.25 as defined by Rahim et al.

        - input: 3D array representing the PET volume in SUV format, with 
                voxels outside the ROI set to NaNs.

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        a = 0.25
        n_voxels = numpy.sum(ROIPet[~numpy.isnan(ROIPet)])
        ROIPet = ROIPet ** a

        return numpy.sum(ROIPet[~numpy.isnan(ROIPet)] / n_voxels) ** (1 / a)

152
153
154
155
156
157
158
159
160
161
162
163
164
    def getSuvPeakFeatureValue(self):
        r"""
        **25. SuvPeak**

        Average of the voxel with maximum SUV within the tumour region and its 26 connected neighbours. 
    
        - input: 3D array representing the PET volume in SUV format
    
        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray
165

166
167
168
169
170
171
172
173
174
175
176
177
        ROIPet[~mask] = numpy.nan

        oc = Oct2Py()

        with Oct2Py() as oc:
            try:
                oc.eval('pkg load image')
                ROIPet = oc.padarray(
                    ROIPet, oc.double(numpy.array([1, 1, 1]).tolist()), numpy.nan
                )

                SUVmax = numpy.max(ROIPet[~numpy.isnan(ROIPet)])
178
                indMax = numpy.where(ROIPet == SUVmax)
179

180
                indMaxX, indMaxY, indMaxZ = indMax[0][0], indMax[1][0], indMax[2][0]
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289

                connectivity = numpy.array(
                    [
                        [-1, -1, -1],
                        [0, -1, -1],
                        [1, -1, -1],
                        [-1, 0, -1],
                        [0, 0, -1],
                        [1, 0, -1],
                        [-1, 1, -1],
                        [0, 1, -1],
                        [1, 1, -1],
                        [-1, -1, 0],
                        [0, -1, 0],
                        [1, -1, 0],
                        [-1, 0, 0],
                        [0, 0, 0],
                        [1, 0, 0],
                        [-1, 1, 0],
                        [0, 1, 0],
                        [1, 1, 0],
                        [-1, -1, 1],
                        [0, -1, 1],
                        [1, -1, 1],
                        [-1, 0, 1],
                        [0, 0, 1],
                        [1, 0, 1],
                        [-1, 1, 1],
                        [0, 1, 1],
                        [1, 1, 1],
                    ]
                )
                nPeak = len(connectivity)
                neighborsMax = numpy.zeros((1, nPeak))

                for i in range(nPeak):
                    neighborsMax[0, i] = ROIPet[
                        connectivity[i, 0] + indMaxX - 1,
                        connectivity[i, 1] + indMaxY - 1,
                        connectivity[i, 2] + indMaxZ - 1,
                    ]

            except Oct2PyError as e:
                print(e)
                oc.exit()

        return numpy.mean(neighborsMax[~numpy.isnan(neighborsMax)])

    def getInactiveVolumeFeatureValue(self):
        r"""
        **26. Inactive volume **

        Percentage of the tumour region that is inactive. Ast defined by Valieres et al
        a threshold of 0.01 × (SUVmax)^2 followed by closing and opening morphological operations were used 
        to differentiate active and inactive regions on FDG-PET scans
        
        - input: 3D array representing the PET volume in SUV format

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray

        thresh = 0.01 * (numpy.max(ROIPet[~numpy.isnan(ROIPet)]) ** 2)

        mask_inactive = ROIPet > thresh

        # MORPHOLOGICAL OPERATIONS
        conn = numpy.zeros([5, 5, 5])
        conn[:, :, 0] = numpy.array(
            [
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
            ]
        )
        conn[:, :, 4] = conn[:, :, 0]

        conn[:, :, 1] = numpy.array(
            [
                [0, 0, 0, 0, 0],
                [0, 0, 1, 0, 0],
                [0, 1, 1, 1, 0],
                [0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0],
            ]
        )
        conn[:, :, 3] = conn[:, :, 1]

        conn[:, :, 2] = numpy.array(
            [
                [0, 0, 1, 0, 0],
                [0, 1, 1, 1, 0],
                [1, 1, 1, 1, 1],
                [0, 1, 1, 1, 0],
                [0, 0, 1, 0, 0],
            ]
        )

        with Oct2Py() as oc:
            try:
                oc.eval('pkg load image')
                perimeter = oc.bwperim(mask, 26)
                mask_inactive = oc.imclose(mask_inactive, conn)
                mask_inactive = oc.imopen(mask_inactive, conn)
                new_mask = mask_inactive + perimeter
Alessia Marcolini's avatar
Alessia Marcolini committed
290

291
292
293
294
295
296
297
                new_mask[new_mask == 1] = 10
                new_mask[new_mask == 2] = 10
                new_mask[new_mask == 0] = 1
                new_mask[new_mask == 10] = 0

                connObjects = oc.bwconncomp(new_mask, 26)

Alessia Marcolini's avatar
Alessia Marcolini committed
298
299
300
301
                b = numpy.zeros((1, int(connObjects.NumObjects)))

                for i in range(int(connObjects.NumObjects)):
                    a = numpy.where(
302
303
304
305
306
307
308
                        len(connObjects.PixelIdxList[i]) >= 15
                    )  # If the number of of pixel forming and object is lower than 15, reject it
                    if a.size == 0:
                        b[i] = 0
                    else:
                        b[i] = 1

Alessia Marcolini's avatar
Alessia Marcolini committed
309
                row, col = numpy.where(b > 0)
310
311
312
313
314
315
316
317
318
319
320
321
322
323
                sumInactive = 0

                for i in range(len(col)):
                    sumInactive = sumInactive + len(
                        connObjects.PixelIdxList[row[i], col[i]]
                    )

                sumVolume = numpy.sum(mask)

            except Oct2PyError as e:
                print(e)
                oc.exit()

        return sumInactive / sumVolume * 100