Hello Dear Developers and Users
I save RT structure in nrrd format (seg.nrrd suffix). Then I run the code in ScriptRepository. I noticed that Layer and LabelValue keys do not exist in header dictionary.
In my case, also the shape of voxels array (the output of read_segmentation_info) is (14, 375, 252, 101) and uint8 type, but the shape of output_voxels array (the output of extract_segments) is (375, 252, 101) and float64 type.
Then I change this code as following:
import nrrd
input_filename = "/home/sn/Ms.Arjmandi/commands/RTSTRUCTStudy1.seg.nrrd"
output_filename = "/home/sn/Ms.Arjmandi/commands/RTSTRUCTStudy1New.seg.nrrd"
segment_list = [("skin",0),("Bladder",1),("Femour_L",2),("Femour_R",3),("Rectum",4),("CTV 1",5),("CTV 2",6),("CTV LN",7),("Penile bulb",8),("PTV 60",9),("PTV 46",10),("PTV 72",11),("Bowel bag",12),("PTV TOTAL",13)]
voxels, header = nrrd.read(input_filename)
def read_segmentation_info(input_filename):
header = nrrd.read_header(input_filename)
segmentation_info = {}
segments = []
segment_index = 0
while True:
prefix = "Segment{0}_".format(segment_index)
if not prefix + "ID" in header.keys():
break
segment = {}
segment["index"] = segment_index
segment["color"] = [float(i) for i in header[prefix + "Color"].split(" ")] # Segment0_Color:=0.501961 0.682353 0.501961
segment["colorAutoGenerated"] = int(header[prefix + "ColorAutoGenerated"]) != 0 # Segment0_ColorAutoGenerated:=1
segment["extent"] = [int(i) for i in header[prefix + "Extent"].split(" ")] # Segment0_Extent:=68 203 53 211 24 118
segment["id"] = header[prefix + "ID"] # Segment0_ID:=Segment_1
segment["name"] = header[prefix + "Name"] # Segment0_Name:=Segment_1
segment["nameAutoGenerated"] = int(header[prefix + "NameAutoGenerated"]) != 0 # Segment0_NameAutoGenerated:=1
#segment["labelValue"] = int(header[prefix + "LabelValue"]) # Segment0_LabelValue:=1
segment["labelValue"] = segment_index
#segment["layer"] = int(header[prefix + "Layer"]) # Segment0_Layer:=0
segment["layer"] = segment_index
tags = {}
tags_str = header[prefix + "Tags"].split("|")
for tag_str in tags_str:
tag_str = tag_str.strip()
if not tag_str:
continue
key, value = tag_str.split(":", maxsplit=1)
tags[key] = value
segment["tags"] = tags
segments.append(segment)
segment_index += 1
segmentation_info["segments"] = segments
return segmentation_info
segmentation_info = read_segmentation_info(input_filename)
def segment_from_name(segmentation_info, segment_name):
for segment in segmentation_info["segments"]:
if segment_name == segment["name"]:
return segment
raise KeyError('segment not found by name ' + segment_name)
def segment_names(segmentation_info):
names = []
for segment in segmentation_info["segments"]:
names.append(segment["name"])
return names
def extract_segments(voxels, header, segmentation_info, segment_list):
import numpy as np
import re
# Create empty array from last 3 dimensions (output will be flattened to a 3D array)
output_voxels = np.zeros(voxels.shape)
# Copy non-segmentation fields to the extracted header
output_header = {}
for key in header.keys():
if not re.match("^Segment[0-9]+_.+", key):
output_header[key] = header[key]
# Copy extracted segments
#dims = len(voxels.shape)
for output_segment_index, segment_name_to_label_value in enumerate(segment_list):
# Copy relabeled voxel data
segment = segment_from_name(segmentation_info, segment_name_to_label_value[0])
#print("segment_name_to_label_value[0]: ", segment_name_to_label_value[0])
input_label_value = segment["labelValue"]
print("input_label_value = ", input_label_value)
output_label_value = segment_name_to_label_value[1]
print("output_label_value = ", output_label_value)
inputLayer = segment["layer"]
print("inputLayer = ", inputLayer)
output_voxels[inputLayer,voxels[inputLayer,:,:,:] == input_label_value] = output_label_value
print("output_voxels size", output_voxels.shape)
# Copy all segment fields corresponding to this segment
for key in header.keys():
prefix = "Segment{0}_".format(segment["index"])
matched = re.match("^"+prefix+"(.+)", key)
if matched:
field_name = matched.groups()[0]
if field_name == "labelValue":
value = output_label_value
elif field_name == "layer":
# output is a single layer (3D volume)
value = 0
else:
value = header[key]
print("output_segment_index = ", output_segment_index)
output_header["Segment{0}_".format(output_segment_index) + field_name] = value
# Remove unnecessary 4th dimension (volume is collapsed into 3D)
# Remove "none" from "none (0,1,0) (0,0,-1) (-1.2999954223632812,0,0)"
output_header["space directions"] = output_header["space directions"][-3:,:]
# Remove "list" from "list domain domain domain"
output_header["kinds"] = output_header["kinds"][-3:]
output_voxels = np.uint8(output_voxels)
return output_voxels, output_header
output_voxels, output_header = extract_segments(voxels, header, segmentation_info, segment_list)
print("final voxels size = ", voxels.shape)
print("final output_voxels size = ", output_voxels.shape)
nrrd.write(output_filename, output_voxels, output_header)
After editing it, when I want to add data (RTSTRUCTStudy1New.seg.nrrd) into 3DSlicer, loading … is paused and get the following message:
Please guide me.
Best regards,
Shahrokh
