libraries/Process/RTstruct.py

import pydicom
import numpy as np
import nibabel as nib
#from matplotlib.path import Path
from PIL import Image, ImageDraw
import scipy
import datetime
import SimpleITK as sitk


def Taubin_smoothing(contour):
    """ Here, we do smoothing in 2D contours!
        Parameters:
            a Nx2 numpy array containing the contour to smooth
        Returns:
            a Nx2 numpy array containing the smoothed contour """
    smoothingloops = 5
    smoothed = [np.empty_like(contour) for i in range(smoothingloops+1)]
    smoothed[0] = contour
    for i in range(smoothingloops):
        # loop over all elements in the contour
        for vertex_i in range(smoothed[0].shape[0]):
            if vertex_i == 0:
                vertex_prev = smoothed[i].shape[0]-1
                vertex_next = vertex_i+1
            elif vertex_i == smoothed[i].shape[0]-1:
                vertex_prev = vertex_i-1
                vertex_next = 0
            else:
                vertex_prev = vertex_i -1
                vertex_next = vertex_i +1
            neighbours_x = np.array([smoothed[i][vertex_prev,0], smoothed[i][vertex_next,0]])
            neighbours_y = np.array([smoothed[i][vertex_prev,1], smoothed[i][vertex_next,1]])
            smoothed[i+1][vertex_i,0] = smoothed[i][vertex_i,0] - 0.3*(smoothed[i][vertex_i,0] - np.mean(neighbours_x))
            smoothed[i+1][vertex_i,1] = smoothed[i][vertex_i,1] - 0.3*(smoothed[i][vertex_i,1] - np.mean(neighbours_y))

    return np.round(smoothed[smoothingloops],3)

class RTstruct:

  def __init__(self):
    self.SeriesInstanceUID = ""
    self.PatientInfo = {}
    self.StudyInfo = {}
    self.CT_SeriesInstanceUID = ""
    self.DcmFile = ""
    self.isLoaded = 0
    self.Contours = []
    self.NumContours = 0
    
    
  def print_struct_info(self, prefix=""):
    print(prefix + "Struct: " + self.SeriesInstanceUID)
    print(prefix + "   " + self.DcmFile)
    
    
  def print_ROINames(self):
    print("RT Struct UID: " + self.SeriesInstanceUID)
    count = -1
    for contour in self.Contours:
      count += 1
      print('  [' + str(count) + ']  ' + contour.ROIName)
    
  def resample_struct(self, newvoxelsize):
    # Rescaling to the newvoxelsize if given in parameter

    for i, Contour in enumerate(self.Contours):
      source_shape = Contour.Mask_GridSize
      voxelsize = Contour.Mask_PixelSpacing
      VoxelX_source = Contour.Mask_Offset[0] + np.arange(source_shape[0])*voxelsize[0]
      VoxelY_source = Contour.Mask_Offset[1] + np.arange(source_shape[1])*voxelsize[1]
      VoxelZ_source = Contour.Mask_Offset[2] + np.arange(source_shape[2])*voxelsize[2]

      target_shape = np.ceil(np.array(source_shape).astype(float)*np.array(voxelsize).astype(float)/newvoxelsize).astype(int)
      VoxelX_target = Contour.Mask_Offset[0] + np.arange(target_shape[0])*newvoxelsize[0]
      VoxelY_target = Contour.Mask_Offset[1] + np.arange(target_shape[1])*newvoxelsize[1]
      VoxelZ_target = Contour.Mask_Offset[2] + np.arange(target_shape[2])*newvoxelsize[2]

      contour = Contour.Mask
      
      if(all(source_shape == target_shape) and np.linalg.norm(np.subtract(voxelsize, newvoxelsize) < 0.001)):
        print("! Image does not need filtering")
      else:
        # anti-aliasing filter
        sigma = [0, 0, 0]
        if(newvoxelsize[0] > voxelsize[0]): sigma[0] = 0.4 * (newvoxelsize[0]/voxelsize[0])
        if(newvoxelsize[1] > voxelsize[1]): sigma[1] = 0.4 * (newvoxelsize[1]/voxelsize[1])
        if(newvoxelsize[2] > voxelsize[2]): sigma[2] = 0.4 * (newvoxelsize[2]/voxelsize[2])
        
        if(sigma != [0, 0, 0]):
            contour = scipy.ndimage.gaussian_filter(contour.astype(float), sigma)
            #come back to binary
            contour[np.where(contour>=0.5)] = 1
            contour[np.where(contour<0.5)] = 0
            
        xi = np.array(np.meshgrid(VoxelX_target, VoxelY_target, VoxelZ_target))
        xi = np.rollaxis(xi, 0, 4)
        xi = xi.reshape((xi.size // 3, 3))
        
        # get resized ct
        contour = scipy.interpolate.interpn((VoxelX_source,VoxelY_source,VoxelZ_source), contour, xi, method='nearest', fill_value=0, bounds_error=False).astype(bool).reshape(target_shape).transpose(1,0,2)
      Contour.Mask_PixelSpacing = newvoxelsize
      Contour.Mask_GridSize = list(contour.shape) 
      Contour.NumVoxels = Contour.Mask_GridSize[0] * Contour.Mask_GridSize[1] * Contour.Mask_GridSize[2]
      Contour.Mask = contour
      self.Contours[i]=Contour
        
  
  def import_Dicom_struct(self, CT):
    if(self.isLoaded == 1):
      print("Warning: RTstruct " + self.SeriesInstanceUID + " is already loaded")
      return 
    dcm = pydicom.dcmread(self.DcmFile)
    
    self.CT_SeriesInstanceUID = CT.SeriesInstanceUID
    
    for dcm_struct in dcm.StructureSetROISequence:  
      ReferencedROI_id = next((x for x, val in enumerate(dcm.ROIContourSequence) if val.ReferencedROINumber == dcm_struct.ROINumber), -1)
      dcm_contour = dcm.ROIContourSequence[ReferencedROI_id]
    
      Contour = ROIcontour()
      Contour.SeriesInstanceUID = self.SeriesInstanceUID
      Contour.ROIName = dcm_struct.ROIName
      Contour.ROIDisplayColor = dcm_contour.ROIDisplayColor
    
      print("Import contour " + str(len(self.Contours)) + ": " + Contour.ROIName)
    
      Contour.Mask = np.zeros((CT.GridSize[0], CT.GridSize[1], CT.GridSize[2]), dtype=np.bool)
      Contour.Mask_GridSize = CT.GridSize
      Contour.Mask_PixelSpacing = CT.PixelSpacing
      Contour.Mask_Offset = CT.ImagePositionPatient
      Contour.Mask_NumVoxels = CT.NumVoxels   
      Contour.ContourMask = np.zeros((CT.GridSize[0], CT.GridSize[1], CT.GridSize[2]), dtype=np.bool)
      
      SOPInstanceUID_match = 1
      
      if not hasattr(dcm_contour, 'ContourSequence'):
          print("The structure [ ", dcm_struct.ROIName ," ] has no attribute ContourSequence. Skipping ...")
          continue

      for dcm_slice in dcm_contour.ContourSequence:
        Slice = {}
      
        # list of Dicom coordinates
        Slice["XY_dcm"] = list(zip( np.array(dcm_slice.ContourData[0::3]), np.array(dcm_slice.ContourData[1::3]) ))
        Slice["Z_dcm"] = float(dcm_slice.ContourData[2])
      
        # list of coordinates in the image frame
        Slice["XY_img"] = list(zip( ((np.array(dcm_slice.ContourData[0::3]) - CT.ImagePositionPatient[0]) / CT.PixelSpacing[0]), ((np.array(dcm_slice.ContourData[1::3]) - CT.ImagePositionPatient[1]) / CT.PixelSpacing[1]) ))
        Slice["Z_img"] = (Slice["Z_dcm"] - CT.ImagePositionPatient[2]) / CT.PixelSpacing[2]
        Slice["Slice_id"] = int(round(Slice["Z_img"]))
      
        # convert polygon to mask (based on matplotlib - slow)
        #x, y = np.meshgrid(np.arange(CT.GridSize[0]), np.arange(CT.GridSize[1]))
        #points = np.transpose((x.ravel(), y.ravel()))
        #path = Path(Slice["XY_img"])
        #mask = path.contains_points(points)
        #mask = mask.reshape((CT.GridSize[0], CT.GridSize[1]))
      
        # convert polygon to mask (based on PIL - fast)
        img = Image.new('L', (CT.GridSize[0], CT.GridSize[1]), 0)
        if(len(Slice["XY_img"]) > 1): ImageDraw.Draw(img).polygon(Slice["XY_img"], outline=1, fill=1)
        mask = np.array(img)
        Contour.Mask[:,:,Slice["Slice_id"]] = np.logical_xor(Contour.Mask[:,:,Slice["Slice_id"]], mask)
        
        # do the same, but only keep contour in the mask
        img = Image.new('L', (CT.GridSize[0], CT.GridSize[1]), 0)
        if(len(Slice["XY_img"]) > 1): ImageDraw.Draw(img).polygon(Slice["XY_img"], outline=1, fill=0)
        mask = np.array(img)
        Contour.ContourMask[:,:,Slice["Slice_id"]] = np.logical_or(Contour.ContourMask[:,:,Slice["Slice_id"]], mask)
            
        Contour.ContourSequence.append(Slice)
      
        # check if the contour sequence is imported on the correct CT slice:
        if(hasattr(dcm_slice, 'ContourImageSequence') and CT.SOPInstanceUIDs[Slice["Slice_id"]] != dcm_slice.ContourImageSequence[0].ReferencedSOPInstanceUID):
          SOPInstanceUID_match = 0
      
      if SOPInstanceUID_match != 1:
        print("WARNING: some SOPInstanceUIDs don't match during importation of " + Contour.ROIName + " contour on CT image")
      
      self.Contours.append(Contour)
      self.NumContours += 1
    #print("self.NumContours",self.NumContours, len(self.Contours))
    self.isLoaded = 1

  def load_from_nii(self, struct_nii_path, rtstruct_labels, rtstruct_colors, data_format, flip):
      
    # data_format can be either "integers" or "one-hot" (i.e. 2**i)
          
    #flip if needed
    if flip == True:
        #struct_nib = nib.orientations.flip_axis(struct_nib, axis=0)
        img = sitk.ReadImage(struct_nii_path)
        img2 = sitk.PermuteAxes(img, [1,0,2])
        sitk.WriteImage(img2, struct_nii_path)
        
    # load the nii image 
    struct_nib = nib.load(struct_nii_path)
    struct_data = struct_nib.get_fdata()
        
    # load the file according to data_format
    if data_format == 'integers':
        roinumbers = np.unique(struct_data)
        roinumbers = roinumbers[roinumbers > 0] # assumes we have a background = 0
        print(roinumbers)
    elif data_format == 'one-hot':
        roinumbers = list(np.arange(np.floor(np.log2(np.max(struct_data))).astype(int)+1)) # CAREFUL WITH THIS LINE, MIGHT NOT WORK IF WE HAVE OVERLAP OF STRUCTURES
    
    # get contourexists from header or compute it from labels and data values
    if hasattr(struct_nib.header, 'extensions') and len(struct_nib.header.extensions) != 0: 
        contoursexist = list(struct_nib.header.extensions[0].get_content())
    else:
        if data_format == 'integers':
            contoursexist = []
            for i in range(len(rtstruct_labels)):
                if i+1 in roinumbers:
                    contoursexist.append(1)
                else:
                    contoursexist.append(0)
        elif data_format == 'one-hot': # AGAIN, THIS ONLY WORKS IF WE DONT HAVE OVERLAP
            contoursexist = []
            for i in range(len(rtstruct_labels)):
                if 2**i in roinumbers:
                    contoursexist.append(1)
                else:
                    contoursexist.append(0)
    
    # get number of rois in struct_data    
    nb_rois_in_struct = len(roinumbers)
    self.NumContours = len(rtstruct_labels)
    # check that they match
    if not len(rtstruct_labels) == len(contoursexist) == nb_rois_in_struct:
        #raise TypeError("The number or struct labels, contoursexist, and  masks in struct.nii.gz is not the same")
        raise Warning("The number or struct labels, contoursexist, and estimated masks in struct.nii.gz is not the same. Taking len(rtstruct_labels) as number of rois")

    # fill in contours
    #TODO fill in ContourSequence and ContourData to be faster later in writeDicomRTstruct
    for c in range(self.NumContours):
        
        Contour = ROIcontour()
        Contour.SeriesInstanceUID = self.SeriesInstanceUID
        Contour.ROIName = rtstruct_labels[c]
        if rtstruct_colors[c] == None:
            Contour.ROIDisplayColor = [0, 0, 255] # default color is blue
        else:
            Contour.ROIDisplayColor = rtstruct_colors[c] 
        if contoursexist[c] == 0:
            Contour.Mask = np.zeros((struct_nib.header['dim'][1], struct_nib.header['dim'][2], struct_nib.header['dim'][3]), dtype=np.bool)
        else:
            if data_format == 'integers':
                Contour.Mask = (struct_data == c+1).astype(bool)
            elif data_format == 'one-hot':
                Contour.Mask = np.bitwise_and(struct_data.astype(int), 2 ** c).astype(bool)
        #TODO enable option for consecutive integers masks?
        Contour.Mask_GridSize = [struct_nib.header['dim'][1], struct_nib.header['dim'][2], struct_nib.header['dim'][3]]
        Contour.Mask_PixelSpacing = [struct_nib.header['pixdim'][1], struct_nib.header['pixdim'][2], struct_nib.header['pixdim'][3]]
        Contour.Mask_Offset = [struct_nib.header['qoffset_x'], struct_nib.header['qoffset_y'], struct_nib.header['qoffset_z']]
        Contour.Mask_NumVoxels = struct_nib.header['dim'][1].astype(int) * struct_nib.header['dim'][2].astype(int) * struct_nib.header['dim'][3].astype(int) 
        # Contour.ContourMask --> this should be only the contour, so far we don't need it so I'll skip it
      
        # apend to self
        self.Contours.append(Contour)
        

  def export_Dicom(self, refCT, outputFile):   
                 
    # meta data
    
    # generate UID
    #uid_base = '' #TODO define one for us if we want? Siri is using: uid_base='1.2.826.0.1.3680043.10.230.',
    # personal UID, applied for via https://www.medicalconnections.co.uk/FreeUID/
    
    SOPInstanceUID = pydicom.uid.generate_uid() #TODO verify this! Siri was using a uid_base, this line is taken from OpenTPS writeRTPlan
    #SOPInstanceUID = pydicom.uid.generate_uid('1.2.840.10008.5.1.4.1.1.481.3.') # siri's version
    
    meta = pydicom.dataset.FileMetaDataset()
    meta.MediaStorageSOPClassUID = '1.2.840.10008.5.1.4.1.1.481.3' # UID class for RTSTRUCT
    meta.MediaStorageSOPInstanceUID = SOPInstanceUID
    # meta.ImplementationClassUID = uid_base + '1.1.1' # Siri's
    meta.ImplementationClassUID =  '1.2.250.1.59.3.0.3.5.0' # from OpenREGGUI
    meta.TransferSyntaxUID = '1.2.840.10008.1.2' # Siri's and OpenREGGUI
    meta.FileMetaInformationGroupLength = 188 # from Siri
    # meta.ImplementationVersionName = 'DCIE 2.2' # from Siri
    
    
    # Main data elements - only required fields, optional fields like StudyDescription are not included for simplicity
    ds = pydicom.dataset.FileDataset(outputFile, {}, file_meta=meta, preamble=b"\0" * 128) # preamble is taken from this example https://pydicom.github.io/pydicom/dev/auto_examples/input_output/plot_write_dicom.html#sphx-glr-auto-examples-input-output-plot-write-dicom-py
    
    # Patient info - will take it from the referenced CT image
    ds.PatientName = refCT.PatientInfo.PatientName
    ds.PatientID = refCT.PatientInfo.PatientID
    ds.PatientBirthDate = refCT.PatientInfo.PatientBirthDate
    ds.PatientSex = refCT.PatientInfo.PatientSex
    
    # General Study 
    dt = datetime.datetime.now()
    ds.StudyDate = dt.strftime('%Y%m%d')
    ds.StudyTime = dt.strftime('%H%M%S.%f')
    ds.AccessionNumber = '1' # A RIS/PACS (Radiology Information System/picture archiving and communication system) generated number that identifies the order for the Study.
    ds.ReferringPhysicianName = 'NA'
    ds.StudyInstanceUID = refCT.StudyInfo.StudyInstanceUID # get from reference CT to indicate that they belong to the same study
    ds.StudyID = refCT.StudyInfo.StudyID # get from reference CT to indicate that they belong to the same study
    
    # RT Series
    #ds.SeriesDate # optional
    #ds.SeriesTime # optional
    ds.Modality = 'RTSTRUCT'
    ds.SeriesDescription = 'AI-predicted' + dt.strftime('%Y%m%d') + dt.strftime('%H%M%S.%f')
    ds.OperatorsName = 'MIRO AI team'
    ds.SeriesInstanceUID = pydicom.uid.generate_uid() # if we have a uid_base --> pydicom.uid.generate_uid(uid_base)
    ds.SeriesNumber = '1'
    
    # General Equipment
    ds.Manufacturer = 'MIRO lab'
    #ds.InstitutionName = 'MIRO lab' # optional
    #ds.ManufacturerModelName = 'nnUNet' # optional, but can be a good tag to insert the model information or label
    #ds.SoftwareVersions # optional, but can be used to insert the version of the code in ECHARP or the version of the model
    
    # Frame of Reference
    ds.FrameOfReferenceUID = refCT.FrameOfReferenceUID
    ds.PositionReferenceIndicator = '' # empty if unknown - info here https://dicom.innolitics.com/ciods/rt-structure-set/frame-of-reference/00201040
    
    # Structure Set
    ds.StructureSetLabel = 'AI predicted' # do not use - or spetial characters or the Dicom Validation in Raystation will give a warning
    #ds.StructureSetName # optional
    #ds.StructureSetDescription # optional
    ds.StructureSetDate = dt.strftime('%Y%m%d')
    ds.StructureSetTime = dt.strftime('%H%M%S.%f')
    ds.ReferencedFrameOfReferenceSequence = pydicom.Sequence()# optional
    # we assume there is only one, the CT
    dssr = pydicom.Dataset()
    dssr.FrameOfReferenceUID = refCT.FrameOfReferenceUID
    dssr.RTReferencedStudySequence = pydicom.Sequence()
    # fill in sequence
    dssr_refStudy = pydicom.Dataset()
    dssr_refStudy.ReferencedSOPClassUID = '1.2.840.10008.3.1.2.3.1' # Study Management Detached
    dssr_refStudy.ReferencedSOPInstanceUID = refCT.StudyInfo.StudyInstanceUID
    dssr_refStudy.RTReferencedSeriesSequence = pydicom.Sequence()
    #initialize
    dssr_refStudy_series = pydicom.Dataset()
    dssr_refStudy_series.SeriesInstanceUID = refCT.SeriesInstanceUID
    dssr_refStudy_series.ContourImageSequence = pydicom.Sequence()
    # loop over slices of CT
    for slc in range(len(refCT.SOPInstanceUIDs)):
        dssr_refStudy_series_slc = pydicom.Dataset()
        dssr_refStudy_series_slc.ReferencedSOPClassUID = refCT.SOPClassUID
        dssr_refStudy_series_slc.ReferencedSOPInstanceUID = refCT.SOPInstanceUIDs[slc]
        # append
        dssr_refStudy_series.ContourImageSequence.append(dssr_refStudy_series_slc)
    
    # append
    dssr_refStudy.RTReferencedSeriesSequence.append(dssr_refStudy_series)
    # append
    dssr.RTReferencedStudySequence.append(dssr_refStudy)
    #append
    ds.ReferencedFrameOfReferenceSequence.append(dssr)
    #
    ds.StructureSetROISequence = pydicom.Sequence()   
    # loop over the ROIs to fill in the fields
    for iroi in range(self.NumContours):
        # initialize the Dataset        
        dssr = pydicom.Dataset()
        dssr.ROINumber = iroi + 1 # because iroi starts at zero and ROINumber cannot be zero
        dssr.ReferencedFrameOfReferenceUID = ds.FrameOfReferenceUID # coming from refCT
        dssr.ROIName = self.Contours[iroi].ROIName
        #dssr.ROIDescription # optional 
        dssr.ROIGenerationAlgorithm = 'AUTOMATIC' # can also be 'SEMIAUTOMATIC' OR 'MANUAL', info here https://dicom.innolitics.com/ciods/rt-structure-set/structure-set/30060020/30060036
        #TODO enable a function to tell us which type of GenerationAlgorithm we have    
        ds.StructureSetROISequence.append(dssr)

    # delete to remove space
    del dssr
    
    #TODO merge all loops into one to be faster, although like this the code is easier to follow I find
    
    # ROI Contour
    ds.ROIContourSequence = pydicom.Sequence()
    # loop over the ROIs to fill in the fields
    for iroi in range(self.NumContours):
        # initialize the Dataset
        dssr = pydicom.Dataset()
        dssr.ROIDisplayColor = self.Contours[iroi].ROIDisplayColor
        dssr.ReferencedROINumber = iroi + 1 # because iroi starts at zero and ReferencedROINumber cannot be zero
        dssr.ContourSequence = pydicom.Sequence() 
        # mask to polygon
        polygonMeshList = self.Contours[iroi].getROIContour()
        # get z vector
        z_coords = list(np.arange(self.Contours[iroi].Mask_Offset[2],self.Contours[iroi].Mask_Offset[2]+self.Contours[iroi].Mask_GridSize[2]*self.Contours[iroi].Mask_PixelSpacing[2], self.Contours[iroi].Mask_PixelSpacing[2]))
        # loop over the polygonMeshList to fill in ContourSequence
        for polygon in polygonMeshList:

            # initialize the Dataset
            dssr_slc = pydicom.Dataset()
            dssr_slc.ContourGeometricType = 'CLOSED_PLANAR' # can also be 'POINT', 'OPEN_PLANAR', 'OPEN_NONPLANAR', info here https://dicom.innolitics.com/ciods/rt-structure-set/roi-contour/30060039/30060040/30060042
            #TODO enable the proper selection of the ContourGeometricType

            # fill in contour points and data
            dssr_slc.NumberOfContourPoints = len(polygon[0::3])
            #dssr_slc.ContourNumber # optional
            # Smooth contour
            smoothed_array_2D = Taubin_smoothing(np.transpose(np.array([polygon[0::3],polygon[1::3]])))
            # fill in smoothed contour
            polygon[0::3] = smoothed_array_2D[:,0]
            polygon[1::3] = smoothed_array_2D[:,1]
            dssr_slc.ContourData = polygon
            
            #get slice
            polygon_z = polygon[2]
            slc = z_coords.index(polygon_z)
            # fill in ContourImageSequence
            dssr_slc.ContourImageSequence = pydicom.Sequence() # Sequence of images containing the contour
            # in our case, we assume we only have one, the reference CT (refCT)
            dssr_slc_ref = pydicom.Dataset()
            dssr_slc_ref.ReferencedSOPClassUID = refCT.SOPClassUID
            dssr_slc_ref.ReferencedSOPInstanceUID = refCT.SOPInstanceUIDs[slc]
            dssr_slc.ContourImageSequence.append(dssr_slc_ref)
              
            # append Dataset to Sequence
            dssr.ContourSequence.append(dssr_slc)
            
        # append Dataset
        ds.ROIContourSequence.append(dssr)
    
    # RT ROI Observations
    ds.RTROIObservationsSequence = pydicom.Sequence()
    # loop over the ROIs to fill in the fields
    for iroi in range(self.NumContours):
        # initialize the Dataset
        dssr = pydicom.Dataset()
        dssr.ObservationNumber = iroi + 1 # because iroi starts at zero and ReferencedROINumber cannot be zero
        dssr.ReferencedROINumber = iroi + 1 ## because iroi starts at zero and ReferencedROINumber cannot be zero
        dssr.ROIObservationLabel = self.Contours[iroi].ROIName #optional
        dssr.RTROIInterpretedType = 'ORGAN' # we can have many types, see here https://dicom.innolitics.com/ciods/rt-structure-set/rt-roi-observations/30060080/300600a4
        # TODO enable a better fill in of the RTROIInterpretedType
        dssr.ROIInterpreter = '' # empty if unknown
        # append Dataset
        ds.RTROIObservationsSequence.append(dssr)
    
    # Approval
    ds.ApprovalStatus = 'UNAPPROVED'#'APPROVED' 
    # if ds.ApprovalStatus = 'APPROVED', then we need to fill in the reviewer information
    #ds.ReviewDate = dt.strftime('%Y%m%d')
    #ds.ReviewTime = dt.strftime('%H%M%S.%f')
    #ds.ReviewerName = 'MIRO AI team'
    
    # SOP common
    ds.SpecificCharacterSet = 'ISO_IR 100' # conditionally required - see info here https://dicom.innolitics.com/ciods/rt-structure-set/sop-common/00080005
    #ds.InstanceCreationDate # optional
    #ds.InstanceCreationTime # optional
    ds.SOPClassUID = '1.2.840.10008.5.1.4.1.1.481.3' #RTSTRUCT file
    ds.SOPInstanceUID = SOPInstanceUID# Siri's --> pydicom.uid.generate_uid(uid_base)
    #ds.InstanceNumber # optional
    
    # save dicom file
    print("Export dicom RTSTRUCT: " + outputFile)
    ds.save_as(outputFile)



      
class ROIcontour:

    def __init__(self):
      self.SeriesInstanceUID = ""
      self.ROIName = ""
      self.ContourSequence = []
      
    def getROIContour(self): # this is from new version of OpenTPS, I(ana) have adapted it to work with old version of self.Contours[i].Mask
    
        try:
            from skimage.measure import label, find_contours
            from skimage.segmentation import find_boundaries
        except:
            print('Module skimage (scikit-image) not installed, ROIMask cannot be converted to ROIContour')
            return 0
    
        polygonMeshList = []
        for zSlice in range(self.Mask.shape[2]):
    
            labeledImg, numberOfLabel = label(self.Mask[:, :, zSlice], return_num=True)
    
            for i in range(1, numberOfLabel + 1):
    
                singleLabelImg = labeledImg == i
                contours = find_contours(singleLabelImg.astype(np.uint8), level=0.6)
    
                if len(contours) > 0:
    
                    if len(contours) == 2:
    
                        ## use a different threshold in the case of an interior contour
                        contours2 = find_contours(singleLabelImg.astype(np.uint8), level=0.4)
    
                        interiorContour = contours2[1]
                        polygonMesh = []
                        for point in interiorContour:
    
                            #xCoord = np.round(point[1]) * self.Mask_PixelSpacing[1] + self.Mask_Offset[1] # original Damien in OpenTPS
                            #yCoord = np.round(point[0]) * self.Mask_PixelSpacing[0] + self.Mask_Offset[0] # original Damien in OpenTPS
                            xCoord = np.round(point[1]) * self.Mask_PixelSpacing[0] + self.Mask_Offset[0] #AB
                            yCoord = np.round(point[0]) * self.Mask_PixelSpacing[1] + self.Mask_Offset[1] #AB
                            zCoord = zSlice * self.Mask_PixelSpacing[2] + self.Mask_Offset[2]
    
                            #polygonMesh.append(yCoord) # original Damien in OpenTPS
                            #polygonMesh.append(xCoord) # original Damien in OpenTPS
                            polygonMesh.append(xCoord) # AB
                            polygonMesh.append(yCoord) # AB
                            polygonMesh.append(zCoord)
    
                        polygonMeshList.append(polygonMesh)
    
                    contour = contours[0]
    
                    polygonMesh = []
                    for point in contour:
    
                        #xCoord = np.round(point[1]) * self.Mask_PixelSpacing[1] + self.Mask_Offset[1] # original Damien in OpenTPS
                        #yCoord = np.round(point[0]) * self.Mask_PixelSpacing[0] + self.Mask_Offset[0] # original Damien in OpenTPS
                        xCoord = np.round(point[1]) * self.Mask_PixelSpacing[0] + self.Mask_Offset[0] #AB
                        yCoord = np.round(point[0]) * self.Mask_PixelSpacing[1] + self.Mask_Offset[1] #AB
                        zCoord = zSlice * self.Mask_PixelSpacing[2] + self.Mask_Offset[2]
    
                        polygonMesh.append(xCoord) # AB
                        polygonMesh.append(yCoord) # AB
                        #polygonMesh.append(yCoord) # original Damien in OpenTPS
                        #polygonMesh.append(xCoord) # original Damien in OpenTPS
                        polygonMesh.append(zCoord)
    
                    polygonMeshList.append(polygonMesh)
    
        ## I (ana) will comment this part since I will not use the class ROIContour for simplicity ###
        #from opentps.core.data._roiContour import ROIContour  ## this is done here to avoir circular imports issue
        #contour = ROIContour(name=self.ROIName, displayColor=self.ROIDisplayColor)
        #contour.polygonMesh = polygonMeshList
    
        #return contour
        
        # instead returning the polygonMeshList directly
        return polygonMeshList