File size: 19,482 Bytes
7f24887
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
2
3
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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
131
132
133
134
135
136
137
138
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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
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
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
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
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
#-- Import libraries
import os
import argparse
import json

# Numerical and Data Handling
import numpy as np
import pandas as pd

# Medical Imaging
import SimpleITK as sitk
import radiomics
from radiomics import featureextractor

# Machine Learning & Clustering
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
import skfuzzy as fuzz

# Image Processing & Segmentation
import scipy.ndimage as ndimage
from skimage.filters import threshold_otsu
from skimage.segmentation import watershed
from skimage.feature import peak_local_max
from skimage import morphology
from scipy.ndimage import distance_transform_edt, binary_erosion




def make_bold(text):
    return f"\033[1m{text}\033[0m"

def load_itk_image(filename):
    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)
    numpyOrigin = itkimage.GetOrigin()
    numpySpacing = itkimage.GetSpacing()
    return numpyImage, numpyOrigin, numpySpacing

def normalize_image_to_uint8(image, lower_bound=-1000, upper_bound=100):
    clipped_img = np.clip(image, lower_bound, upper_bound)
    normalized_img = ((clipped_img - lower_bound) / (upper_bound - lower_bound)) * 255.0
    normalized_img = normalized_img.astype(np.uint8)
    return normalized_img


def segment_nodule_kmeans(ct_image, bbox_center, bbox_whd, margin=5, n_clusters=2):
    """
    Segments a nodule in a 3D CT image using k-means clustering with a margin around the bounding box.

    Parameters:
    - ct_image: 3D NumPy array representing the CT image.
    - bbox_center: Tuple of (x, y, z) coordinates for the center of the bounding box.
    - bbox_whd: Tuple of (w, h, d) representing the width, height, and depth of the bounding box.
    - margin: Margin to add around the bounding box (default is 5).
    - n_clusters: Number of clusters to use in k-means (default is 2).

    Returns:
    - segmented_image: 3D NumPy array with the segmented nodule.
    """

    x_center, y_center, z_center = bbox_center
    w, h, d = bbox_whd

    # Calculate the bounding box with margin
    x_start, x_end = max(0, x_center - w//2 - margin), min(ct_image.shape[0], x_center + w//2 + margin)
    y_start, y_end = max(0, y_center - h//2 - margin), min(ct_image.shape[1], y_center + h//2 + margin)
    z_start, z_end = max(0, z_center - d//2 - margin), min(ct_image.shape[2], z_center + d//2 + margin)

    bbox_region = ct_image[x_start:x_end, y_start:y_end, z_start:z_end]

    # Reshape the region for k-means clustering
    flat_region = bbox_region.reshape(-1, 1)

    # Perform k-means clustering
    kmeans = KMeans(n_clusters=n_clusters, n_init=10,random_state=0).fit(flat_region)
    labels = kmeans.labels_

    # Reshape the labels back to the original bounding box shape
    clustered_region = labels.reshape(bbox_region.shape)

    # Assume the nodule is in the cluster with the highest mean intensity
    nodule_cluster = np.argmax(kmeans.cluster_centers_)

    # Create a binary mask for the nodule
    nodule_mask = (clustered_region == nodule_cluster)

    # Apply morphological operations to refine the segmentation
    nodule_mask = ndimage.binary_closing(nodule_mask, structure=np.ones((3, 3, 3)))
    nodule_mask = ndimage.binary_opening(nodule_mask, structure=np.ones((2, 2, 2)))

    # Initialize the segmented image
    segmented_image = np.zeros_like(ct_image, dtype=np.uint8)

    # Place the nodule mask in the correct position in the segmented image
    segmented_image[x_start:x_end, y_start:y_end, z_start:z_end] = nodule_mask

    return segmented_image


def segment_nodule_gmm(ct_image, bbox_center, bbox_whd, margin=5, n_components=2):
    """
    Segments a nodule in a 3D CT image using a Gaussian Mixture Model with a margin around the bounding box.

    Parameters:
    - ct_image: 3D NumPy array representing the CT image.
    - bbox_center: Tuple of (x, y, z) coordinates for the center of the bounding box.
    - bbox_whd: Tuple of (w, h, d) representing the width, height, and depth of the bounding box.
    - margin: Margin to add around the bounding box (default is 5).
    - n_components: Number of components to use in the Gaussian Mixture Model (default is 2).

    Returns:
    - segmented_image: 3D NumPy array with the segmented nodule.
    """

    x_center, y_center, z_center = bbox_center
    w, h, d = bbox_whd

    # Calculate the bounding box with margin
    x_start, x_end = max(0, x_center - w//2 - margin), min(ct_image.shape[0], x_center + w//2 + margin)
    y_start, y_end = max(0, y_center - h//2 - margin), min(ct_image.shape[1], y_center + h//2 + margin)
    z_start, z_end = max(0, z_center - d//2 - margin), min(ct_image.shape[2], z_center + d//2 + margin)

    bbox_region = ct_image[x_start:x_end, y_start:y_end, z_start:z_end]

    # Reshape the region for GMM
    flat_region = bbox_region.reshape(-1, 1)

    # Perform GMM
    gmm = GaussianMixture(n_components=n_components, random_state=0).fit(flat_region)
    labels = gmm.predict(flat_region)

    # Reshape the labels back to the original bounding box shape
    clustered_region = labels.reshape(bbox_region.shape)

    # Assume the nodule is in the component with the highest mean intensity
    nodule_component = np.argmax(gmm.means_)

    # Create a binary mask for the nodule
    nodule_mask = (clustered_region == nodule_component)

    # Apply morphological operations to refine the segmentation
    nodule_mask = ndimage.binary_closing(nodule_mask, structure=np.ones((3, 3, 3)))
    nodule_mask = ndimage.binary_opening(nodule_mask, structure=np.ones((3, 3, 3)))

    # Initialize the segmented image
    segmented_image = np.zeros_like(ct_image, dtype=np.uint8)

    # Place the nodule mask in the correct position in the segmented image
    segmented_image[x_start:x_end, y_start:y_end, z_start:z_end] = nodule_mask

    return segmented_image


def segment_nodule_fcm(ct_image, bbox_center, bbox_whd, margin=5, n_clusters=2):
    """
    Segments a nodule in a 3D CT image using Fuzzy C-means clustering with a margin around the bounding box.

    Parameters:
    - ct_image: 3D NumPy array representing the CT image.
    - bbox_center: Tuple of (x, y, z) coordinates for the center of the bounding box.
    - bbox_whd: Tuple of (w, h, d) representing the width, height, and depth of the bounding box.
    - margin: Margin to add around the bounding box (default is 5).
    - n_clusters: Number of clusters to use in Fuzzy C-means (default is 2).

    Returns:
    - segmented_image: 3D NumPy array with the segmented nodule.
    """

    x_center, y_center, z_center = bbox_center
    w, h, d = bbox_whd

    # Calculate the bounding box with margin
    x_start, x_end = max(0, x_center - w//2 - margin), min(ct_image.shape[0], x_center + w//2 + margin)
    y_start, y_end = max(0, y_center - h//2 - margin), min(ct_image.shape[1], y_center + h//2 + margin)
    z_start, z_end = max(0, z_center - d//2 - margin), min(ct_image.shape[2], z_center + d//2 + margin)

    bbox_region = ct_image[x_start:x_end, y_start:y_end, z_start:z_end]

    # Reshape the region for FCM
    flat_region = bbox_region.reshape(-1, 1)

    # Perform FCM clustering
    cntr, u, u0, d, jm, p, fpc = fuzz.cluster.cmeans(flat_region.T, n_clusters, 2, error=0.005, maxiter=1000, init=None)

    # Assign each voxel to the cluster with the highest membership
    labels = np.argmax(u, axis=0)

    # Reshape the labels back to the original bounding box shape
    clustered_region = labels.reshape(bbox_region.shape)

    # Assume the nodule is in the cluster with the highest mean intensity
    nodule_cluster = np.argmax(cntr)

    # Create a binary mask for the nodule
    nodule_mask = (clustered_region == nodule_cluster)

    # Apply morphological operations to refine the segmentation
    nodule_mask = ndimage.binary_closing(nodule_mask, structure=np.ones((3, 3, 3)))
    nodule_mask = ndimage.binary_opening(nodule_mask, structure=np.ones((3, 3, 3)))

    # Initialize the segmented image
    segmented_image = np.zeros_like(ct_image, dtype=np.uint8)

    # Place the nodule mask in the correct position in the segmented image
    segmented_image[x_start:x_end, y_start:y_end, z_start:z_end] = nodule_mask

    return segmented_image



def segment_nodule_otsu(ct_image, bbox_center, bbox_whd, margin=5):
    """
    Segments a nodule in a 3D CT image using Otsu's thresholding with a margin around the bounding box.

    Parameters:
    - ct_image: 3D NumPy array representing the CT image.
    - bbox_center: Tuple of (x, y, z) coordinates for the center of the bounding box.
    - bbox_whd: Tuple of (w, h, d) representing the width, height, and depth of the bounding box.
    - margin: Margin to add around the bounding box (default is 5).

    Returns:
    - segmented_image: 3D NumPy array with the segmented nodule.
    """

    x_center, y_center, z_center = bbox_center
    w, h, d = bbox_whd

    # Calculate the bounding box with margin
    x_start, x_end = max(0, x_center - w//2 - margin), min(ct_image.shape[0], x_center + w//2 + margin)
    y_start, y_end = max(0, y_center - h//2 - margin), min(ct_image.shape[1], y_center + h//2 + margin)
    z_start, z_end = max(0, z_center - d//2 - margin), min(ct_image.shape[2], z_center + d//2 + margin)

    bbox_region = ct_image[x_start:x_end, y_start:y_end, z_start:z_end]

    # Flatten the region for thresholding
    flat_region = bbox_region.flatten()

    # Calculate the Otsu threshold
    otsu_threshold = threshold_otsu(flat_region)

    # Apply the threshold to create a binary mask
    nodule_mask = bbox_region >= otsu_threshold

    # Apply morphological operations to refine the segmentation
    nodule_mask = ndimage.binary_closing(nodule_mask, structure=np.ones((3, 3, 3)))
    nodule_mask = ndimage.binary_opening(nodule_mask, structure=np.ones((3, 3, 3)))

    # Initialize the segmented image
    segmented_image = np.zeros_like(ct_image, dtype=np.uint8)

    # Place the nodule mask in the correct position in the segmented image
    segmented_image[x_start:x_end, y_start:y_end, z_start:z_end] = nodule_mask

    return segmented_image



def expand_mask_by_distance(segmented_nodule_gmm, spacing, expansion_mm):
    """
    Expands the segmentation mask by a given distance in mm in all directions by directly updating pixel values.

    Parameters:
    segmented_nodule_gmm (numpy array): 3D binary mask of the nodule (1 for nodule, 0 for background).
    spacing (tuple): Spacing of the image in mm for each voxel, given as (spacing_x, spacing_y, spacing_z).
    expansion_mm (float): Distance to expand the mask in millimeters.

    Returns:
    numpy array: Expanded segmentation mask.
    """
    # Reorder spacing to match the numpy array's (z, y, x) format
    spacing_reordered = (spacing[2], spacing[1], spacing[0])  # (spacing_z, spacing_y, spacing_x)

    # Calculate the number of pixels to expand in each dimension
    expand_pixels = np.array([int(np.round(expansion_mm / s)) for s in spacing_reordered])

    # Create a new expanded mask with the same shape
    expanded_mask = np.zeros_like(segmented_nodule_gmm)

    # Get the coordinates of all white pixels in the original mask
    white_pixel_coords = np.argwhere(segmented_nodule_gmm == 1)

    # Expand each white pixel by adding the specified number of pixels in each direction
    for coord in white_pixel_coords:
        z, y, x = coord  # Extract the z, y, x coordinates of each white pixel

        # Define the range to expand for each coordinate
        z_range = range(max(0, z - expand_pixels[0]), min(segmented_nodule_gmm.shape[0], z + expand_pixels[0] + 1))
        y_range = range(max(0, y - expand_pixels[1]), min(segmented_nodule_gmm.shape[1], y + expand_pixels[1] + 1))
        x_range = range(max(0, x - expand_pixels[2]), min(segmented_nodule_gmm.shape[2], x + expand_pixels[2] + 1))

        # Update the new mask by setting all pixels in this range to 1
        for z_new in z_range:
            for y_new in y_range:
                for x_new in x_range:
                    expanded_mask[z_new, y_new, x_new] = 1

    return expanded_mask


def find_nodule_lobe(cccwhd, lung_mask, class_map):
    """
    Determine the lung lobe where a nodule is located based on a 3D mask and bounding box.

    Parameters:
    cccwhd (list or tuple): Bounding box in the format [center_x, center_y, center_z, width, height, depth].
    lung_mask (numpy array): 3D array representing the lung mask with different lung regions.
    class_map (dict): Dictionary mapping lung region labels to their names.

    Returns:
    str: Name of the lung lobe where the nodule is located.
    """
    center_x, center_y, center_z, width, height, depth = cccwhd

    # Calculate the bounding box limits
    start_x = int(center_x - width // 2)
    end_x = int(center_x + width // 2)
    start_y = int(center_y - height // 2)
    end_y = int(center_y + height // 2)
    start_z = int(center_z - depth // 2)
    end_z = int(center_z + depth // 2)

    # Ensure the indices are within the mask dimensions
    start_x = max(0, start_x)
    end_x = min(lung_mask.shape[0], end_x)
    start_y = max(0, start_y)
    end_y = min(lung_mask.shape[1], end_y)
    start_z = max(0, start_z)
    end_z = min(lung_mask.shape[2], end_z)

    # Extract the region of interest (ROI) from the mask
    roi = lung_mask[start_x:end_x, start_y:end_y, start_z:end_z]

    # Count the occurrences of each lobe label within the ROI
    unique, counts = np.unique(roi, return_counts=True)
    label_counts = dict(zip(unique, counts))

    # Exclude the background (label 0)
    if 0 in label_counts:
        del label_counts[0]

    # Find the label with the maximum count
    if label_counts:
        nodule_lobe = max(label_counts, key=label_counts.get)
    else:
        nodule_lobe = None

    # Map the label to the corresponding lung lobe
    if nodule_lobe is not None:
        nodule_lobe_name = class_map["lungs"][nodule_lobe]
    else:
        nodule_lobe_name = "Undefined"

    return nodule_lobe_name


def find_nodule_lobe_and_distance(cccwhd, lung_mask, class_map,spacing):
    """
    Determine the lung lobe where a nodule is located and measure its distance from the lung wall.

    Parameters:
    cccwhd (list or tuple): Bounding box in the format [center_x, center_y, center_z, width, height, depth].
    lung_mask (numpy array): 3D array representing the lung mask with different lung regions.
    class_map (dict): Dictionary mapping lung region labels to their names.

    Returns:
    tuple: (Name of the lung lobe, Distance from the lung wall)
    """
    center_x, center_y, center_z, width, height, depth = cccwhd

    # Calculate the bounding box limits
    start_x = int(center_x - width // 2)
    end_x = int(center_x + width // 2)
    start_y = int(center_y - height // 2)
    end_y = int(center_y + height // 2)
    start_z = int(center_z - depth // 2)
    end_z = int(center_z + depth // 2)

    # Ensure the indices are within the mask dimensions
    start_x = max(0, start_x)
    end_x = min(lung_mask.shape[0], end_x)
    start_y = max(0, start_y)
    end_y = min(lung_mask.shape[1], end_y)
    start_z = max(0, start_z)
    end_z = min(lung_mask.shape[2], end_z)

    # Extract the region of interest (ROI) from the mask
    roi = lung_mask[start_x:end_x, start_y:end_y, start_z:end_z]

    # Count the occurrences of each lobe label within the ROI
    unique, counts = np.unique(roi, return_counts=True)
    label_counts = dict(zip(unique, counts))

    # Exclude the background (label 0)
    if 0 in label_counts:
        del label_counts[0]

    # Find the label with the maximum count
    if label_counts:
        nodule_lobe = max(label_counts, key=label_counts.get)
    else:
        nodule_lobe = None

    # Map the label to the corresponding lung lobe
    if nodule_lobe is not None:
        nodule_lobe_name = class_map["lungs"][nodule_lobe]
    else:
        nodule_lobe_name = "Undefined"

    # Calculate the distance from the nodule centroid to the nearest lung wall
    nodule_centroid = np.array([center_x, center_y, center_z])
    
    # Create a binary lung mask where lung region is 1 and outside lung is 0
    lung_binary_mask = lung_mask > 0

    # Create the lung wall mask by finding the outer boundary
    # Use binary erosion to shrink the lung mask, then subtract it from the original mask to get the boundary
    lung_eroded    = binary_erosion(lung_binary_mask)
    lung_wall_mask = lung_binary_mask & ~lung_eroded  # Lung wall mask is the outermost boundary (contour)

    # Compute the distance transform from the lung wall
    distance_transform = distance_transform_edt(~lung_wall_mask)  # Compute distance to nearest lung boundary

    
    
    # Get the distance from the nodule centroid to the nearest lung wall in voxel units
    voxel_distance_to_lung_wall = distance_transform[center_x, center_y, center_z]

    # Convert voxel distance to real-world distance in mm
    physical_distance_to_lung_wall = voxel_distance_to_lung_wall * np.sqrt(
        spacing[0]**2 + spacing[1]**2 + spacing[2]**2
    )



    return nodule_lobe_name, voxel_distance_to_lung_wall,physical_distance_to_lung_wall


# +
def expand_mask_by_distance(segmented_nodule_gmm, spacing, expansion_mm):
    """
    Expands the segmentation mask by a given distance in mm in all directions by directly updating pixel values.

    Parameters:
    segmented_nodule_gmm (numpy array): 3D binary mask of the nodule (1 for nodule, 0 for background).
    spacing (tuple): Spacing of the image in mm for each voxel, given as (spacing_x, spacing_y, spacing_z).
    expansion_mm (float): Distance to expand the mask in millimeters.

    Returns:
    numpy array: Expanded segmentation mask.
    """
    # Reorder spacing to match the numpy array's (z, y, x) format
    spacing_reordered = (spacing[2], spacing[1], spacing[0])  # (spacing_z, spacing_y, spacing_x)

    # Calculate the number of pixels to expand in each dimension
    expand_pixels = np.array([int(np.round(expansion_mm / s)) for s in spacing_reordered])

    # Create a new expanded mask with the same shape
    expanded_mask = np.zeros_like(segmented_nodule_gmm)

    # Get the coordinates of all white pixels in the original mask
    white_pixel_coords = np.argwhere(segmented_nodule_gmm == 1)

    # Expand each white pixel by adding the specified number of pixels in each direction
    for coord in white_pixel_coords:
        z, y, x = coord  # Extract the z, y, x coordinates of each white pixel

        # Define the range to expand for each coordinate
        z_range = range(max(0, z - expand_pixels[0]), min(segmented_nodule_gmm.shape[0], z + expand_pixels[0] + 1))
        y_range = range(max(0, y - expand_pixels[1]), min(segmented_nodule_gmm.shape[1], y + expand_pixels[1] + 1))
        x_range = range(max(0, x - expand_pixels[2]), min(segmented_nodule_gmm.shape[2], x + expand_pixels[2] + 1))

        # Update the new mask by setting all pixels in this range to 1
        for z_new in z_range:
            for y_new in y_range:
                for x_new in x_range:
                    expanded_mask[z_new, y_new, x_new] = 1

    return expanded_mask



# Function to plot the contours of a mask
def plot_contours(ax, mask, color, linewidth=1.5):
    contours = measure.find_contours(mask, level=0.5)  # Find contours at a constant level
    for contour in contours:
        ax.plot(contour[:, 1], contour[:, 0], color=color, linewidth=linewidth)