nii_to_dicom.py

import nibabel as nib
import pydicom
import os
import glob
import numpy as np
from copy import deepcopy
from matplotlib.patches import Polygon
import warnings
from scipy.ndimage import find_objects
from scipy.ndimage.morphology import binary_fill_holes
from skimage import measure
from PIL import Image, ImageDraw
import scipy
import datetime

def convert_nii_to_dicom(dicomctdir, predictedNiiFile, predictedDicomFile, predicted_structures=[], rtstruct_colors=[], refCT = None):
    # img = nib.load(os.path.join(predniidir, patient_id, 'RTStruct.nii.gz'))
    # data = img.get_fdata()[:,:,:,1]
    # patient_list = PatientList() # initialize list of patient data
    # patient_list.list_dicom_files(os.path.join(ct_ref_path,patient,inner_ct_ref_path), 1) # search dicom files in the patient data folder, stores all files in the attributes (all CT images, dose file, struct file)
    # refCT = patient_list.list[0].CTimages[0]
    # refCT.import_Dicom_CT()

    struct = RTstruct()
    struct.load_from_nii(predictedNiiFile, predicted_structures, rtstruct_colors) #TODO add already the refCT info in here because there are fields to do that
    if not struct.Contours[0].Mask_PixelSpacing == refCT.PixelSpacing:
      struct.resample_struct(refCT.PixelSpacing)
    struct.export_Dicom(refCT, predictedDicomFile)

    # create_RT_struct(dicomctdir, data.transpose([1,0,2]).astype(int), dicomdir, predicted_structures)

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
    if newvoxelsize is not None: 
      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("This structure 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_or(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):
      
    # load the nii image 
    struct_nib = nib.load(struct_nii_path)
    struct_data = struct_nib.get_fdata()
            
    # get contourexists from header
    if len(struct_nib.header.extensions)==0:
      contoursexist = []
    else:
      contoursexist = list(struct_nib.header.extensions[0].get_content())
    
    # get number of rois in struct_data 
    # for nii with consecutive integers
    #roinumbers = np.unique(struct_data) 
    # for nii with power of 2 format
    roinumbers = list(np.arange(np.floor(np.log2(np.max(struct_data))).astype(int)+1)) # CAREFUL WITH THIS LINE, MIGHT NOT WORK ALWAYS IF WE HAVE OVERLAP OF 
    nb_rois_in_struct = len(roinumbers)
    
    # check that they match
    if len(contoursexist)!=0 and (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(contoursexist) as number of rois")
        self.NumContours = len(contoursexist)
    else:
        self.NumContours = nb_rois_in_struct
        
    # 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 len(contoursexist)!=0 and 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:
            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 PARROT 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
            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