Add Add SUVPeak and fix PercentageInactive

pyradiomics_radler/myfirstorder.py

 import numpy
 from scipy import ndimage
+from oct2py import Oct2Py, Oct2PyError
 
 from radiomics.firstorder import RadiomicsFirstOrder
 from .mybase import MyRadiomicsFeaturesBase
@@ -148,3 +149,190 @@ class MyRadiomicsFirstOrder(MyRadiomicsFeaturesBase):
 
         return numpy.sum(ROIPet[~numpy.isnan(ROIPet)] / n_voxels) ** (1 / a)
 
+    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
+        ROIPet[~mask] = numpy.nan
+
+        oc = Oct2Py()
+        ROIPet = oc.double(ROIPet.tolist())
+
+        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)])
+                indMax = numpy.where(ROIPet[~numpy.isnan(ROIPet)] == SUVmax)[0]
+
+                [indMaxX, indMaxY, indMaxZ] = numpy.unravel_index(
+                    indMax[0], ROIPet.shape, 'F'
+                )
+                # connectivity = oc.getneighbors(
+                #     oc.strel(
+                #         'arbitrary',
+                #         oc.double(numpy.ones((3, 3, 3), dtype=numpy.int).tolist()),
+                #     )
+                # )
+
+                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,
+                    ]
+                # import ipdb
+
+                # ipdb.set_trace()
+            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
+
+        ROIPet[~mask] = 0
+
+        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
+                new_mask[mask == 0] = 10
+                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)
+                b = numpy.zeros((1, connObjects.NumObjects))
+
+                for i in range(connObjects.NumObjects):
+                    a = oc.find(
+                        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
+
+                [row, col] = oc.find(b > 0)
+
+                sumInactive = 0
+
+                for i in range(len(col)):
+                    sumInactive = sumInactive + len(
+                        connObjects.PixelIdxList[row[i], col[i]]
+                    )
+
+                sumVolume = numpy.sum(mask)
+                import ipdb
+
+                ipdb.set_trace()
+
+            except Oct2PyError as e:
+                print(e)
+                oc.exit()
+
+        return sumInactive / sumVolume * 100