There was a problem merging two.NIIs into one

[tanjia123456]

Hello, I have a dicom sequence, use 3D slicer to load the dicom folder, and then mark labels on it using segmentation, but now I have two mask1 and mask2, how do I merge them into one? Am I using the following code correctly?
`
import numpy as np
import os
import SimpleITK as sitk
import random
from scipy import ndimage
from os.path import join

Rseg = sitk.ReadImage("R-mark-label.nii", sitk.sitkInt8)
Rseg_array = sitk.GetArrayFromImage(Rseg)

Lseg = sitk.ReadImage("L-mark-label.nii", sitk.sitkInt8)
Lseg_array = sitk.GetArrayFromImage(Lseg)
print(" Rseg = ", Rseg_array.shape, "Lseg = ", Lseg_array.shape)

unique_elements = np.unique(Rseg_array, return_counts=False)
print("Rseg 类别数=", unique_elements, "Rseg=0:", np.sum(Rseg_array == 0), "Rseg=1:", np.sum(Rseg_array == 1), "Rseg=2:",np.sum(Rseg_array == 2))

unique_elements = np.unique(Lseg_array, return_counts=False)
print("Lseg 类别数=", unique_elements, "Lseg=0:", np.sum(Lseg_array == 0), "Rseg=1:", np.sum(Lseg_array == 1), "Rseg=2:",np.sum(Lseg_array == 2))

merge_seg = np.add(Rseg_array, Lseg_array)
unique_elements = np.unique(merge_seg, return_counts=False)

print("merge_seg 类别数=", unique_elements, "merge_seg=0:", np.sum(merge_seg == 0), "merge_seg=1:", np.sum(merge_seg == 1), "merge_seg=2:",np.sum(merge_seg == 2))
out = sitk.GetImageFromArray(merge_seg)

sitk.WriteImage(out, 'seg-12.nii')

`

[tanjia123456]

BY the way, why don't the same sequences match up?

`
import numpy as np
import os
import SimpleITK as sitk
import random
from scipy import ndimage
from os.path import join
import matplotlib.pyplot as plt

for i in range(20):
mri_path = 'process1/'+str(i)+'.nii.gz'
seg_path1 = 'process1/seg-'+str(i)+ '.nii'
seg_path2 = 'process1/seg-'+str(i)+ '.nii.gz'

mri = sitk.ReadImage(mri_path, sitk.sitkInt16)
mri_array = sitk.GetArrayFromImage(mri)

if os.path.exists(seg_path1):
    seg = sitk.ReadImage(seg_path1, sitk.sitkInt8)
    seg_array = sitk.GetArrayFromImage(seg)
else:
    seg = sitk.ReadImage(seg_path2, sitk.sitkInt8)
    seg_array = sitk.GetArrayFromImage(seg)
print("\nprocess1: seg_array=", seg_array.shape, "mri_array=", mri_array.shape)




# 找到肝脏区域开始和结束的slice，并各向外扩张
z = np.any(seg_array, axis=(1, 2))
print("z=", z.shape)
# np.where(z)[0] 输出满足条件，即true的元素位置，返回为一个列表，用0，-1取出第一个位置和最后一个位置的值
start_slice, end_slice = np.where(z)[0][[0, -1]]
print("process3: start_slice = ", start_slice, "end_slice = ", end_slice)

mid_slice = np.where(z)[0]

plt.figure(figsize=(10, 10))
plt.subplot(2, 2, 1)
plt.imshow(seg_array[start_slice + int(len(np.where(z)[0])/2), :, :], cmap='gray')
plt.subplot(2, 2, 2)
plt.imshow(mri_array[start_slice + int(len(np.where(z)[0])/2), :, :], cmap='gray')
plt.show()






# 截取保留区域
mri_array = mri_array[start_slice:end_slice + 1, :, :]
seg_array = seg_array[start_slice:end_slice + 1, :, :]
print("Preprocessed shape:", mri_array.shape, seg_array.shape)


# 保存为对应的格式
new_mri = sitk.GetImageFromArray(mri_array)
new_mri.SetDirection(mri.GetDirection())
new_mri.SetOrigin(mri.GetOrigin())
#new_mri.SetSpacing((mri.GetSpacing()[0] * int(1 / self.xy_down_scale),ct.GetSpacing()[1] * int(1 / self.xy_down_scale), self.slice_down_scale))

new_seg = sitk.GetImageFromArray(seg_array)
new_seg.SetDirection(seg.GetDirection())
new_seg.SetOrigin(seg.GetOrigin())
# new_seg.SetSpacing((ct.GetSpacing()[0] * int(1 / self.xy_down_scale),ct.GetSpacing()[1] * int(1 / self.xy_down_scale), self.slice_down_scale))

if new_mri != None and new_seg != None:
    # sitk.WriteImage(new_ct, os.path.join(self.fixed_path, 'ct', ct_file))
    sitk.WriteImage(new_mri, "process2/"+str(i)+".nii.gz")
    sitk.WriteImage(new_seg, "process2/seg-"+str(i)+".nii.gz")`

[tanjia123456]

SimpleITK 2.3.1

[zivy]

Hello @tanjia123456,

It is very hard to follow your code (long and not straightforward), please provide a minimal working example illustrating the problem you are having combining what looks to be two binary segmentation images.

Below are two options for combining two binary segmentations, depending on your needs.

import SimpleITK as sitk

#binary segmentation1
segmentation1 = sitk.Image([3,3], sitk.sitkUInt8)
segmentation1[0:2,0:2] = 1

# binary segmentation2
segmentation2 = sitk.Image([10,10], sitk.sitkUInt8)
segmentation2[1:3,1:3] = 1

# combine into binary segmentation
combined_binary_segmentation = (segmentation1 + segmentation2) != 0
print(sitk.GetArrayFromImage(combined_binary_segmentation))

# combine into segmentation with three arbitrary labels
segmentation1_label = 10
segmentation2_label = 20
intersection_label = 30

combined_multi_label_segmentation = segmentation1*segmentation1_label + segmentation2*segmentation2_label
combined_multi_label_segmentation[combined_multi_label_segmentation == segmentation1_label + segmentation2_label] = intersection_label
print(sitk.GetArrayFromImage(combined_multi_label_segmentation))

Finally, for questions, please post on the ITK discourse. The issue tracker is primarily intended for bug reports.

[tanjia123456]

Hello @tanjia123456,

It is very hard to follow your code (long and not straightforward), please provide a minimal working example illustrating the problem you are having combining what looks to be two binary segmentation images.

Below are two options for combining two binary segmentations, depending on your needs.

import SimpleITK as sitk

#binary segmentation1
segmentation1 = sitk.Image([3,3], sitk.sitkUInt8)
segmentation1[0:2,0:2] = 1

# binary segmentation2
segmentation2 = sitk.Image([10,10], sitk.sitkUInt8)
segmentation2[1:3,1:3] = 1

# combine into binary segmentation
combined_binary_segmentation = (segmentation1 + segmentation2) != 0
print(sitk.GetArrayFromImage(combined_binary_segmentation))

# combine into segmentation with three arbitrary labels
segmentation1_label = 10
segmentation2_label = 20
intersection_label = 30

combined_multi_label_segmentation = segmentation1*segmentation1_label + segmentation2*segmentation2_label
combined_multi_label_segmentation[combined_multi_label_segmentation == segmentation1_label + segmentation2_label] = intersection_label
print(sitk.GetArrayFromImage(combined_multi_label_segmentation))

Finally, for questions, please post on the ITK discourse. The issue tracker is primarily intended for bug reports.

thanks for your code, it's work!

[tanjia123456]

Hi, @zivy
Since the resolution of each subject is different, and some are even slightly crooked, I want to crop the middle cerebral artery (in the red box), then what needs to be done to make the mass cropping more accurate? The cropped code is shown below
'import os.path

import nibabel as nib
from nibabel.viewers import OrthoSlicer3D

def crop_save(nii_file):
img = nib.load(nii_file)
data = img.get_fdata()

x_min, x_max = 30, 30+256  # 替换为你的x轴边界
y_min, y_max = 64, 64+256  # 替换为你的y轴边界
z_min, z_max = 60, 120  # 替换为你的z轴边界

cropped_data = data[x_min:x_max, y_min:y_max, z_min:z_max]
print(data.shape, cropped_data.shape)

cropped_img = nib.Nifti1Image(cropped_data, img.affine)
nib.save(cropped_img, os.path.join('crop', nii_file.split('/')[1]))

for i in range(1,10):
# 加载NII文件
nii_file1 = 'skull_seg/skull_seg_PA' + str(i) + '_SE1'+ '.nii'
nii_file2 = 'skull_seg/regis_PA' + str(i) + '_SE0' + '.nii'
print("\n", nii_file1, nii_file2)
crop_save(nii_file1)
crop_save(nii_file2)
'

[zivy]

Hello @tanjia123456,

You are mixing and matching nibabel and SimpleITK for operations that likely exist in both toolkits. Highly recommend that you use one toolkit and familiarize yourself with all its functionality. For SimpleITK, possibly skim the SimpleITK notebooks or go over the SimpleITK tutorial.

With respect to automated cropping, the following code will do it when added to the code above:

label_shape_filter = sitk.LabelShapeStatisticsImageFilter()
label_shape_filter.Execute(combined_multi_label_segmentation)
# The bounding box's first "dim=2 or 3" entries are the starting index and last "dim=2 or 3" entries the size
bounding_box = label_shape_filter.GetBoundingBox(intersection_label)
cropped_multi_label = sitk.RegionOfInterest(
        combined_multi_label_segmentation,
        bounding_box[int(len(bounding_box) / 2) :],
        bounding_box[0 : int(len(bounding_box) / 2)],
    )
print(sitk.GetArrayFromImage(cropped_multi_label))