CT image acquisition human body part MASK (removal machine)

When making statistics on the mean standard deviation of CT images, we hope to make statistics only for the human body part, so we need to obtain the MASK mask of the human body part from the CT image, and the human body is often connected with the machine tool part, so we also need to remove the machine tool part

Simple implementation method

Ideally, the human body can be completely separated from the machine part through a lower threshold (- 150). In this case, the mask of the human body can be obtained by taking the most general domain.

code:

import numpy as np
import os
import SimpleITK as sitk


def getmaxcomponent(mask_array, num_limit=10):
    # sitk method to obtain connected domain
    cca = sitk.ConnectedComponentImageFilter()
    cca.SetFullyConnected(True)
    cca.FullyConnectedOff()
    _input = sitk.GetImageFromArray(mask_array.astype(np.uint8))
    output_ex = cca.Execute(_input)
    labeled_img = sitk.GetArrayFromImage(output_ex)
    num = cca.GetObjectCount()
    max_label = 0
    max_num = 0
    # It is not necessary to traverse all connected domains. Generally, there is a label corresponding to the whole body mask in front, which reduces the calculation time
    for i in range(1, num_limit):  	
        if np.sum(labeled_img == i) < 1e6:		# The number of voxels in the whole body mask must be large, and those less than the set value will not be considered
            continue
        if np.sum(labeled_img == i) > max_num:
            max_num = np.sum(labeled_img == i)
            max_label = i
    maxcomponent = np.array((labeled_img == max_label)).astype(np.uint8)
    print(str(max_label) + ' num:' + str(max_num))  	# Look, which one is the biggest
    return maxcomponent


def get_body(CT_nii_array):
    """
    card CT Threshold acquisition body (ideally enough, but in most cases it will be included in the machine part)
    """
    # The threshold is binarized to obtain the largest 3d connected domain
    CT_array = np.copy(CT_nii_array)
    threshold_all = -200  # Card threshold, card out the whole body and machine part
    CT_array[CT_array >= threshold_all] = 1
    CT_array[CT_array <= threshold_all] = 0
    body_mask1 = getmaxcomponent(CT_array, 10)
    return body_mask1.astype(np.uint8)

Solutions to undesirable situations

The unsatisfactory situation means that the human body is connected to the machine tool, so the mask obtained after taking the connected domain actually contains the organic bed part. This situation also commonly occurs in practice. Therefore, it is necessary to remove the machine tool part. The following two methods are tried:

1. Morphological method:
Idea: disconnect the operation and corrosion, hoping to disconnect the connection between the human body and the machine tool, and then take the most general domain again
Disadvantages: this method has certain effect when the connection between human and machine tool is not too serious, but it is easy to cause some problems after morphological processing

2. (the last method adopted) first obtain the mask of the machine tool:
Idea: in the case of multiple CT data, after each CT image obtains the largest general domain (possibly including the machine tool), we overlay all connected domains. Because the positions of the human body and the machine tool are different, and the position of the machine tool is fixed, the overlapping times of the connected places are not higher than that of the machine tool, and we reserve the area where the overlapping times are greater than the set threshold, Then remove the most common domain (i.e. human body part) at this time. Finally, after a simple morphological processing, the mask of the machine tool part can be obtained

Get the code of the machine mask (incomplete, you need to write it yourself to read the data):

body_mask_final = np.zeros([374, 512, 512], dtype=np.uint8)
for i in train_id:
    i = str(i)
    print('\r'+i, end='')
    CT_nii_file = data_folder + sep + i + '_CT.nii.gz'
    # SimpleITK reads nii
    CT_nii = sitk.ReadImage(CT_nii_file)
    CT_a = sitk.GetArrayFromImage(CT_nii)

    body_mask = get_body(CT_a)
    body_mask_final += body_mask
   
# To the most general domain
maxcomponent = getmaxcomponent(machine_mask)

# Remove the largest body part
machine_mask = machine_mask - maxcomponent

# Expand 3 * 3 and take the most general domain
machine_mask_d = dilate_ct(machine_mask, 3, other_axis=False)
machine_mask_s = getmaxcomponent(machine_mask_d)

save_path = r".\machine_mask.nii.gz"		# preservation
save_sitk = sitk.GetImageFromArray(machine_mask_s.astype(np.uint8))
sitk.WriteImage(save_sitk, save_path)

Then cooperate with the initial method to obtain the mask with only the human body

code:

def cut_machine(CT_nii_array, CT_nii, nii_folder_path, id):
    """
    from CT The human body part after removing the machine part is obtained in the image
    """
    body_mask_file = nii_folder_path + sep + id + r'_body_mask.nii.gz'
    machine_mask_file = r'.\machine_mask.nii.gz'  # Machine mask file
    # The threshold is binarized to obtain the largest 3d connected domain
    body_mask1 = get_body(CT_nii_array)
    # Post processing
    # Remove the machine part
    machine_mask = sitk.ReadImage(machine_mask_file)
    machine_mask_nii = sitk.GetArrayFromImage(machine_mask)
    body_mask2 = body_mask1 - machine_mask_nii
    body_mask2[body_mask2 == 255] = 0
    body_mask2[body_mask2 < 0] = 0
    body_mask3 = getmaxcomponent(body_mask2, 150)
    # fill
    body_mask_fill = np.zeros(body_mask3.shape, dtype=np.uint8)
    for j in range(body_mask2.shape[0]):
        body_mask_fill[j, :, :] = fill_inter_bone(body_mask3[j, :, :])
    save_mask = body_mask_fill
    save_nii(save_mask, CT_nii, body_mask_file)
    # plt.imshow(save_nii[195, :, :])
    # plt.show()
    return save_mask.astype(np.uint8)

def save_nii(save_nii, CT_nii, save_path, save_mask=True):
    """
    preservation nii
    :param save_nii: Need to save nii Pictorial array
    :param CT_nii: Registered nii Image, used to obtain the same information
    :param save_path: Save path
    :param save_mask: Saved Yes No mask，The default is True
    :return:
    """
    if save_mask:
        # Save mask_nii
        save_sitk = sitk.GetImageFromArray(save_nii.astype(np.uint8))
    else:
        # Save img_nii
        save_sitk = sitk.GetImageFromArray(save_nii.astype(np.float))
    save_sitk.CopyInformation(CT_nii)
    sitk.WriteImage(save_sitk, save_path)
    print(save_path + ' processing successfully!')

def fill_inter_bone(mask):
    # Fill a hole in an image and read in a layer
    mask = mask_fill = mask.astype(np.uint8)
    if np.sum(mask[:]) != 0:  # That is, the read in layer has a value
        contours, hierarchy = cv2.findContours(mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
        len_contour = len(contours)
        contour_list = []
        for i in range(len_contour):
            drawing = np.zeros_like(mask, np.uint8)  # create a black image
            img_contour = cv2.drawContours(drawing, contours, i, (1, 1, 1), -1)
            contour_list.append(img_contour)
        mask_fill = sum(contour_list)
        mask_fill[mask_fill >= 1] = 1
    return mask_fill.astype(np.uint8)