<jupyter_start><jupyter_text>Definitions<jupyter_code># Path to dicom files
PathDicom = 'DICOM/HNSCC'
# Dataset (tfrecord) output path
OUTPUT_PATH = './'

# Dataset Full-Resolution size
Objective_size = np.array([256,256,512])
# Number of down-samples dataset to be created.
# The smallest dataset will have a size of Objective_size/(2**(DOWNSAMPLE_SCALES-1))
DOWNSAMPLE_SCALES = 8

# Fraction of dataset used for validation
FRAC_VALIDATION = 0.10
# Random selection of validation samples (Hard-Coded for reproducibility)
# Leave empty to randomly pick the samples.
validation_patients_ID = ['HNSCC-01-0024', 
                          'HNSCC-01-0097', 
                          'HNSCC-01X-0169', 
                          'HNSCC-01-0032', 
                          'HNSCC-01-0062', 
                          'HNSCC-01-0148', 
                          'HNSCC-01-0137']


# Set different thresholds for the segmentation of the CT images
# Thresholds are in Hounsfield Units
HU_Air_limit = -1000
HU_Lung_upper_limit = -500
HU_Fluids_Fat_lower_limit = -125
HU_Fluids_Fat_upper_limit = 10.0
HU_Parenchyma_upper_limit = 90.0
HU_Parenchyma_lower_limit = 10.0
HU_Bone_low_limit = 90.0
HU_Bone_upper_limit = 1300<jupyter_output><empty_output><jupyter_text>List samples<jupyter_code>ruta, datasetName = os.path.split(PathDicom)

sample_old = ''
name_old = ''
current_type = ''

no_attenuation_corrected_names = ['PETNAC', 'no ac', 'no_ac', 'not ac', 'not_ac']
CT_names = ['ct atten cor', 'CT Soft', 'CT Atten', 'CT Body']

# Create dataframe
col_names =  ['PatientID', 'Study', 'Type', 'Slices', 'Axial_size','Location']
HNSCC_dataframe = pd.DataFrame(columns=col_names) 

study_count = -1
for dirName, subdirList, fileList in os.walk(PathDicom):

    # If it is the base directory, we get the samples names
    ruta, aux_name = os.path.split(dirName)
    if aux_name == datasetName:
        sample_name_list = subdirList
        print('Found %d samples'%len(sample_name_list))

    if len(subdirList) != 0:
        continue

    # Get name of current sample
    base_aux, name_sample = os.path.split(dirName)
    base_aux, name_study = os.path.split(base_aux)
    if 'simulation' in name_study.lower():
        continue
    base_aux, sample_act = os.path.split(base_aux)

    if sample_act != sample_old:
        print('Scanning sample: %s'%sample_act)
        sample_old = sample_act


    if name_study != name_old:
        print('\tStudy: %s'%name_study)
        name_old = name_study


    if any(nac_name.lower() in name_sample.lower() for nac_name in no_attenuation_corrected_names):
        sys.stdout.write('\t\tNAC PET sample:')
        current_type = 'PET_NO_AC'
    elif any(ct_name.lower() in name_sample.lower() for ct_name in CT_names):
        sys.stdout.write('\t\tCorr CT sample:')
        current_type = 'CT'
    else:
        continue

    axial_size = DH.get_sample_axial_FOV_limits(dirName)

    # Print sample and sample size
    print('\t\t\t%s -- %d'%(name_sample, len(fileList)))
    print('\t\t\t\t\t\tAxial_FoV -- [%0.2f ; %0.2f]'%(axial_size[0], axial_size[1]))


    # Add to dataframe
    HNSCC_dataframe.loc[len(HNSCC_dataframe)] = [sample_act, name_study, current_type, len(fileList), axial_size, dirName]<jupyter_output><empty_output><jupyter_text>Clean dataset from non-matched samples<jupyter_code>patient_list = HNSCC_dataframe.PatientID.unique()

# Hand removed studies
study_ban_list = ['06-03-2002-PETCT HEAD  NECK CA-86406']


for patient in patient_list:

    # get patient dataframe
    patient_df = HNSCC_dataframe.loc[HNSCC_dataframe['PatientID'] == patient]

    # get all studies of this patient
    study_list = patient_df.Study.unique()

    for study_id in study_list:

        kill_data_point = False

        # get study dataframe
        study_df = patient_df.loc[patient_df['Study'] == study_id]

        # If the number of unique image types is not at least 2, erase.
        if len(study_df.Type.unique()) < 2:
            kill_data_point = True

        elif any(ban_name.lower() in study_id.lower() for ban_name in study_ban_list):
            kill_data_point = True

        else:          
            # Check the number of study types
            pet_df = study_df.loc[study_df['Type'] == 'PET_NO_AC']
            ct_df = study_df.loc[study_df['Type'] == 'CT']

            num_ct = len(ct_df.index)
            num_pet = len(pet_df.index)

            # Get the axial FoV limits of each image
            ct_axial_limits = np.zeros((num_ct,2))
            for idx_ct in range(num_ct):
                ct_axial_limits[idx_ct,:] = ct_df.Axial_size.values[idx_ct]
            pet_axial_limits = np.zeros((num_pet,2))
            for idx_pet in range(num_pet):
                pet_axial_limits[idx_pet,:] = pet_df.Axial_size.values[idx_pet]


            # get the FoV likelihood
            simil = np.zeros((num_ct,num_pet))
            for idx_ct in range(num_ct):
                for idx_pet in range(num_pet):
                    simil[idx_ct,idx_pet] = (ct_axial_limits[idx_ct,0]-pet_axial_limits[idx_pet,0])**2 + (ct_axial_limits[idx_ct,1]-pet_axial_limits[idx_pet,1])**2

            # Get the most similar studies
            ct_min, pet_min = np.unravel_index(np.argmin(simil), simil.shape)

            # Get the minimal overlap limit
            lim_inf = np.max([ct_axial_limits[ct_min,0],pet_axial_limits[pet_min,0]])
            lim_sup = np.min([ct_axial_limits[ct_min,1],pet_axial_limits[pet_min,1]])

            # If they are not conected delete entry...
            if lim_inf>lim_sup:
                kill_data_point = True
            else:
                # Kill other studies datapoints 
                HNSCC_dataframe = HNSCC_dataframe.drop(HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'CT') & (HNSCC_dataframe.Location != ct_df.Location.values[ct_min]) ])
                HNSCC_dataframe = HNSCC_dataframe.drop(HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'PET_NO_AC') & (HNSCC_dataframe.Location != pet_df.Location.values[pet_min]) ])


        if kill_data_point:
            HNSCC_dataframe = HNSCC_dataframe.drop(HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) ])<jupyter_output><empty_output><jupyter_text>Read in samplesDuring this process we complete a new database with the volume information before converting and normalizing volumes.<jupyter_code>col_names =  ['PatientID',
              'Study',
              'CT_Vol_Shape_X', 
              'CT_Vol_Shape_Y', 
              'CT_Vol_Shape_Z', 
              'CT_Voxel_size_X', 
              'CT_Voxel_size_Y', 
              'CT_Voxel_size_Z', 
              'CT_FoV_Size_X', 
              'CT_FoV_Size_Y', 
              'CT_FoV_Size_Z', 
              'CT_FoV_X_min',
              'CT_FoV_X_max',
              'CT_FoV_Y_min',
              'CT_FoV_Y_max',
              'CT_FoV_Z_min',
              'CT_FoV_Z_max',
              'PET_Vol_Shape_X', 
              'PET_Vol_Shape_Y',
              'PET_Vol_Shape_Z',
              'PET_Voxel_size_X', 
              'PET_Voxel_size_Y',
              'PET_Voxel_size_Z',
              'PET_FoV_Size_X',
              'PET_FoV_Size_Y',
              'PET_FoV_Size_Z',
              'PET_FoV_X_min',
              'PET_FoV_X_max',
              'PET_FoV_Y_min',
              'PET_FoV_Y_max',
              'PET_FoV_Z_min',
              'PET_FoV_Z_max']


DICOM_dataframe = pd.DataFrame(columns=col_names)


# For each patient
patient_list = HNSCC_dataframe.PatientID.unique()
idx_pat = 0
for patient in patient_list:

    idx_pat = idx_pat + 1
    print('(%d/%d) Loading patient: %s'%(idx_pat,len(patient_list), patient))


    # and each study
    patient_df = HNSCC_dataframe.loc[HNSCC_dataframe['PatientID'] == patient]
    # get all studies of this patient
    study_list = patient_df.Study.unique()

    for study_id in study_list:
        # And create one entry...

        ct_path = HNSCC_dataframe.Location.loc[HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'CT')  ]].values[0]
        pet_path = HNSCC_dataframe.Location.loc[HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'PET_NO_AC')  ]].values[0]

        # Get all slice paths
        list_paths_ct = sorted(os.listdir(ct_path))
        list_paths_pet = sorted(os.listdir(pet_path))


        # Read first slice and allocate memory
        RefDs = dicom.read_file(os.path.join(ct_path,list_paths_ct[0]))
        ct_cols = RefDs.Columns
        ct_rows = RefDs.Rows
        ct_slices = len(list_paths_ct)
        ct_x_size = RefDs.PixelSpacing[0]
        ct_y_size = RefDs.PixelSpacing[1]
        ct_z_size = RefDs.SliceThickness
        ct_FOV_X = ct_cols*ct_x_size
        ct_FOV_Y = ct_rows*ct_y_size
        ct_FOV_Z = ct_slices*ct_z_size
        ct_FOV_X_min = -ct_FOV_X/2.0
        ct_FOV_X_max = ct_FOV_X/2.0
        ct_FOV_Y_min = -ct_FOV_Y/2.0
        ct_FOV_Y_max = ct_FOV_Y/2.0
#         ct_FOV_Z_min = -ct_FOV_Z/2.0
#         ct_FOV_Z_max = ct_FOV_Z/2.0
        axial_aux = HNSCC_dataframe.Axial_size.loc[HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'CT')  ]].values[0]
        ct_FOV_Z_min = axial_aux[0]
        ct_FOV_Z_max = axial_aux[1]



        RefDs = dicom.read_file(os.path.join(pet_path,list_paths_pet[0]))
        pet_cols = RefDs.Columns
        pet_rows = RefDs.Rows
        pet_slices = len(list_paths_pet)
        pet_x_size = RefDs.PixelSpacing[0]
        pet_y_size = RefDs.PixelSpacing[1]
        pet_z_size = ct_z_size
        pet_FOV_X = pet_cols*pet_x_size
        pet_FOV_Y = pet_rows*pet_y_size
        pet_FOV_Z = pet_slices*pet_z_size
        pet_FOV_X_min = -pet_FOV_X/2.0
        pet_FOV_X_max = pet_FOV_X/2.0
        pet_FOV_Y_min = -pet_FOV_Y/2.0
        pet_FOV_Y_max = pet_FOV_Y/2.0
#         pet_FOV_Z_min = -pet_FOV_Z/2.0
#         pet_FOV_Z_max = pet_FOV_Z/2.0
        axial_aux = HNSCC_dataframe.Axial_size.loc[HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'PET_NO_AC')  ]].values[0]
        pet_FOV_Z_min = axial_aux[0]
        pet_FOV_Z_max = axial_aux[1]

        # Add to dataframe
        DICOM_dataframe.loc[len(DICOM_dataframe)] = [patient, 
                                                     study_id,
                                                     ct_cols,
                                                     ct_rows, 
                                                     ct_slices,
                                                     ct_x_size, 
                                                     ct_y_size, 
                                                     ct_z_size,
                                                     ct_FOV_X,
                                                     ct_FOV_Y,
                                                     ct_FOV_Z, 
                                                     ct_FOV_X_min,
                                                     ct_FOV_X_max,
                                                     ct_FOV_Y_min,
                                                     ct_FOV_Y_max,
                                                     ct_FOV_Z_min,
                                                     ct_FOV_Z_max,
                                                     pet_cols, 
                                                     pet_rows, 
                                                     pet_slices, 
                                                     pet_x_size, 
                                                     pet_y_size, 
                                                     pet_z_size,
                                                     pet_FOV_X,
                                                     pet_FOV_Y,
                                                     pet_FOV_Z,
                                                     pet_FOV_X_min,
                                                     pet_FOV_X_max,
                                                     pet_FOV_Y_min,
                                                     pet_FOV_Y_max,
                                                     pet_FOV_Z_min,
                                                     pet_FOV_Z_max]


print('CT FoV Limits:\n \t(%0.2f , %0.2f)\n \t(%0.2f , %0.2f)\n \t(%0.2f , %0.2f)'%(DICOM_dataframe.CT_FoV_X_min.min(),DICOM_dataframe.CT_FoV_X_max.max(),
                                                                        DICOM_dataframe.CT_FoV_Y_min.min(),DICOM_dataframe.CT_FoV_Y_max.max(),
                                                                        DICOM_dataframe.CT_FoV_Z_min.min(),DICOM_dataframe.CT_FoV_Z_max.max()))

print('PET FoV Limits:\n \t(%0.2f , %0.2f)\n \t(%0.2f , %0.2f)\n \t(%0.2f , %0.2f)'%(DICOM_dataframe.PET_FoV_X_min.min(),DICOM_dataframe.PET_FoV_X_max.max(),
                                                                        DICOM_dataframe.PET_FoV_Y_min.min(),DICOM_dataframe.PET_FoV_Y_max.max(),
                                                                        DICOM_dataframe.PET_FoV_Z_min.min(),DICOM_dataframe.PET_FoV_Z_max.max()))<jupyter_output><empty_output><jupyter_text>Load Dataset and Normalization<jupyter_code># Get absolute FoV limits
Lim_X_min = min([DICOM_dataframe.CT_FoV_X_min.min(),DICOM_dataframe.PET_FoV_X_min.min()])
Lim_Y_min = min([DICOM_dataframe.CT_FoV_Y_min.min(),DICOM_dataframe.PET_FoV_Y_min.min()])
Lim_Z_min = min([DICOM_dataframe.CT_FoV_Z_min.min(),DICOM_dataframe.PET_FoV_Z_min.min()])

Lim_X_max = max([DICOM_dataframe.CT_FoV_X_max.max(),DICOM_dataframe.PET_FoV_X_max.max()])
Lim_Y_max = max([DICOM_dataframe.CT_FoV_Y_max.max(),DICOM_dataframe.PET_FoV_Y_max.max()])
Lim_Z_max = max([DICOM_dataframe.CT_FoV_Z_max.max(),DICOM_dataframe.PET_FoV_Z_max.max()])

# Calculate voxel size
Voxel_size = [(Lim_X_max-Lim_X_min)/float(Objective_size[0]),
              (Lim_Y_max-Lim_Y_min)/float(Objective_size[1]),
              (Lim_Z_max-Lim_Z_min)/float(Objective_size[2])]

# Transform dataset



print('Volume size:')
print(Objective_size)
print('Voxel size:')
print(Voxel_size)

low_res_post_scaling = 2<jupyter_output><empty_output><jupyter_text>Divide dataset<jupyter_code>unique_patients_list = DICOM_dataframe.PatientID.unique()
num_uniques = len(unique_patients_list)
print("%d unique patients."%num_uniques)

num_validation = np.ceil(num_uniques*FRAC_VALIDATION)

# If no validation patient list was selected, they are randomly chosen
if len(validation_patients_ID) == 0:
    # Draw the valiation patients IDs at random
    validation_patients_ID = np.random.choice(unique_patients_list, num_validation)

    

print("\nvaliation Patients IDs:")
print(validation_patients_ID)

# Create the partition entry
DICOM_dataframe["Partition"] = ""
DICOM_dataframe["Partition"] = "Train"

# Replace the valiation patients with the valiation flag
DICOM_dataframe.loc[DICOM_dataframe['PatientID'].isin(validation_patients_ID), "Partition"] = "Validation"

# Print info
num_sample_train = np.sum(DICOM_dataframe.Partition == "Train")
num_sample_valiation = np.sum(DICOM_dataframe.Partition == "Validation")
print("\n%d train patients. %d samples."%(num_uniques-num_validation, num_sample_train))
print("%d validation patients. %d samples."%(num_validation, num_sample_valiation))
print("\nvalidation/train fraction = %f."%(num_sample_valiation/num_sample_train))<jupyter_output><empty_output><jupyter_text>Preprocess<jupyter_code># The couch removal was not perfect, some samples must be fine tunned:
couch_removal_skip = list()
list_conf_1 = ['HNSCC-01-0027','HNSCC-01-0107']
list_conf_2 = ['HNSCC-01-0121','HNSCC-01-0095','HNSCC-01-0107','HNSCC-01-0081','HNSCC-01-0118','HNSCC-01-0212','HNSCC-01-0137','HNSCC-01-0091','HNSCC-01-0093']
list_conf_3 = ['HNSCC-01-0091']
# Create all sub-sampled datasets
write_train_list = list()
write_validation_list = list()


OUTPUT_NAME_TRAIN = os.path.join(OUTPUT_PATH, 'Train_Dataset')
OUTPUT_NAME_VALIDATION = os.path.join(OUTPUT_PATH, 'Validation_Dataset')

for idx_downSample in range(DOWNSAMPLE_SCALES):

    scale_factor = 2**idx_downSample

    NAME_APPEND = "_%dx%dx%d"%(Objective_size[0]/scale_factor,
                                Objective_size[1]/scale_factor,
                                Objective_size[2]/scale_factor)

    write_train = tf.io.TFRecordWriter(OUTPUT_NAME_TRAIN+NAME_APPEND+'.tfrecord')
    write_validation = tf.io.TFRecordWriter(OUTPUT_NAME_VALIDATION+NAME_APPEND+'.tfrecord')

    write_train_list.append(write_train)
    write_validation_list.append(write_validation)
order_zoom_use = 0

tic = time.time()

col_names =  ['PatientID',
              'Location',
              'CT_min',
              'CT_max',
              'CT_5_lim',
              'PET_min',
              'PET_max',
              'PET_5_lim',
              'PET_bone_mean',
              'PET_parenchyma_mean',
              'PET_fluids-fat_mean',
              'PET_air_mean',
              'X_size',
              'Y_size',
              'Z_size']


DATASET_frame = pd.DataFrame(columns=col_names)

# Limits margin
HU_Lung_upper_limit = HU_Lung_upper_limit
HU_Bone_upper_limit = HU_Bone_upper_limit
HU_Fluids_Fat_lower_limit = HU_Fluids_Fat_lower_limit * 1.2

# Used to clean the image histograms
PORC_HIST_LIM_PET = 0.01
PORC_HIST_LIM_CT = 0.01

# CT dynamic range limits
LIM_SUP_CT = (((HU_Bone_upper_limit/1000.0)+1)*DH.mu_agua_120kev/10.0)
LIM_INF_CT = (((HU_Fluids_Fat_lower_limit/1000.0)+1)*DH.mu_agua_120kev/10.0) 

# Histogram analysis
Lim_bines = 3
Lim_bines_test = 35
LIM_SUP_CT

# For each patient
patient_list = DICOM_dataframe.PatientID.unique()
idx_pat = 0
for patient in patient_list:

    idx_pat = idx_pat + 1
    print('(%d/%d) Loading patient: %s'%(idx_pat,len(patient_list), patient))


    # and each study
    patient_df = DICOM_dataframe.loc[DICOM_dataframe['PatientID'] == patient]
    # get all studies of this patient
    study_list = patient_df.Study.unique()

    study_number = 0
    for study_id in study_list:
        # And create one entry...

        tac = time.time()

        study_df = patient_df.loc[patient_df['Study'] == study_id]

        ct_path = HNSCC_dataframe.Location.loc[HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'CT')  ]].values[0]
        pet_path = HNSCC_dataframe.Location.loc[HNSCC_dataframe.index[(HNSCC_dataframe.PatientID == patient) & (HNSCC_dataframe.Study == study_id) & (HNSCC_dataframe.Type == 'PET_NO_AC')  ]].values[0]

        # Get all slice paths
        list_paths_ct = sorted(os.listdir(ct_path))
        list_paths_pet = sorted(os.listdir(pet_path))

        # Allocate memory for input
        ct_input = np.zeros((study_df.CT_Vol_Shape_X.values[0],
                             study_df.CT_Vol_Shape_Y.values[0],
                             study_df.CT_Vol_Shape_Z.values[0]), dtype=np.float32)
        pet_input = np.zeros((study_df.PET_Vol_Shape_X.values[0],
                              study_df.PET_Vol_Shape_Y.values[0],
                              study_df.PET_Vol_Shape_Z.values[0]), dtype=np.float32)

        # Load volumes
        locations_CT = np.zeros((study_df.CT_Vol_Shape_Z.values[0]))
        for idx, filename in enumerate(list_paths_ct):
            if ".dcm" in filename.lower():  # check whether the file's DICOM
                RefDs = dicom.read_file(os.path.join(ct_path,filename))
                ct_input[:,:,idx] = (RefDs.pixel_array*RefDs.RescaleSlope)+RefDs.RescaleIntercept
                locations_CT[idx] = RefDs.SliceLocation

        # Reorder data to its actual position in space
        orden_fix = np.argsort(locations_CT)
        ct_input = ct_input[:,:,orden_fix]

        # Remove coach from CT
        if patient in couch_removal_skip or (patient == 'HNSCC-01-0091' and study_number == 1):
            print_ct_orig = np.copy(ct_input[:,:,int(ct_input.shape[2]/2)])
            print_ct_final = np.copy(ct_input[:,:,int(ct_input.shape[2]/2)])
            couch_mask = np.zeros(print_ct_orig.shape)
            print('No couch removal possible.')
        else:
            B_n_d = 4
            B_n_e = 2
            B_d_d = 15
            B_d_e = 13
            limits_ct_use = [0,0]
            corr_disk = 6
            if patient == list_conf_1:
                B_n_d = 0
                B_n_e = 0
            if patient in list_conf_2:
                B_n_d = 2
                B_n_e = 2
            if patient in list_conf_3:
                limits_ct_use[1] = int(ct_input.shape[2]/2)
            if patient == 'HNSCC-01-0114' and study_number == 1:
                corr_disk = 3


            couch_mask = DH.CT_couch_removal_mask_multicore(ct_input, 
                                                               samples_use = 20, 
                                                               n_corr_disk = corr_disk, 
                                                               limits_CT = limits_ct_use,
                                                               B_n_dilatation_disk_rad = B_n_d,
                                                               B_n_erosion_disk_rad = B_n_e,
                                                               B_d_dilatation_disk_rad = B_d_d,
                                                               B_d_erosion_disk_rad = B_d_e) 
            print_ct_orig = np.copy(ct_input[:,:,int(ct_input.shape[2]/2)])
            for z_idx in range(0,ct_input.shape[2]):
                ct_input[:,:,z_idx][couch_mask] = HU_Air_limit
            print_ct_final = np.copy(ct_input[:,:,int(ct_input.shape[2]/2)])

        # Convert from HU tu linear attenuation coeficient
        ct_input = ((ct_input/1000.0)+1)*DH.mu_agua_120kev/10.0
        ct_input[ct_input<0]=0

        locations_PET = np.zeros((study_df.PET_Vol_Shape_Z.values[0]))
        time_reference = np.zeros((study_df.PET_Vol_Shape_Z.values[0]))
        for idx, filename in enumerate(list_paths_pet):
            if ".dcm" in filename.lower():  # check whether the file's DICOM
                RefDs = dicom.read_file(os.path.join(pet_path,filename))
                if RefDs.Units == 'PROPCNTS':
                    pet_input[:,:,idx] = ((RefDs.pixel_array*RefDs.RescaleSlope)+RefDs.RescaleIntercept) / (RefDs.ActualFrameDuration/1000.0)
                else:
                    pet_input[:,:,idx] = ((RefDs.pixel_array*RefDs.RescaleSlope)+RefDs.RescaleIntercept) #*RefDs.DecayFactor
                locations_PET[idx] = RefDs.SliceLocation
                time_reference[idx] = RefDs.FrameReferenceTime

        # PET data to SUVlbm
        if RefDs.PatientSex == 'M':
            # SUVlbm(Janma): males: 9.27E3 * weight / (6.68E3 + 216 * weight / (height^2))
            pet_input = pet_input/(9.27e3 * RefDs.PatientWeight / (6.68e3 + 216 * RefDs.PatientWeight / (RefDs.PatientSize**2)))
        else:
            # SUVlbm(Janma): females: 9.27E3 * weight / (8.78E3 + 244 * weight / (height^2))
            pet_input = pet_input/(9.27e3 * RefDs.PatientWeight / (8.75e3 + 244 * RefDs.PatientWeight / (RefDs.PatientSize**2)))
        
        # Reorder data to its actual position in space
        orden_fix = np.argsort(locations_PET)
        pet_input = pet_input[:,:,orden_fix]

        pet_input[pet_input<0]=0
        
        
        # Get voxels size
        locations_PET = np.sort(locations_PET)
        locations_CT = np.sort(locations_CT)
        Voxel_size_INPUT = np.zeros((2,3), dtype=np.float32)
        Voxel_size_INPUT[0,0] = study_df.CT_Voxel_size_X.values[0]
        Voxel_size_INPUT[0,1] = study_df.CT_Voxel_size_Y.values[0]
        Voxel_size_INPUT[0,2] = (locations_CT[1]-locations_CT[0])#study_df.CT_Voxel_size_Z.values[0]
        Voxel_size_INPUT[1,0] = study_df.PET_Voxel_size_X.values[0]
        Voxel_size_INPUT[1,1] = study_df.PET_Voxel_size_Y.values[0]
        Voxel_size_INPUT[1,2] = (locations_PET[1]-locations_PET[0])#study_df.PET_Voxel_size_Z.values[0]

        
        

        # Normalize volume size
        ct_output, pet_output, limits_Z = DH.normalize_volume_size(ct_input, 
                                                                      pet_input, 
                                                                      locations_CT, 
                                                                      locations_PET, 
                                                                      Voxel_size_INPUT, 
                                                                      Voxel_size,
                                                                      Objective_size)


        # Calculate dataframe info

        # Get CT histogram at tiesue data location           
        bines_ct = np.linspace(LIM_INF_CT,LIM_SUP_CT*2,128)
        vals_hst_ct, _ = np.histogram(ct_output.flatten(), bins=bines_ct)
        ct_low_lim = bines_ct[np.argwhere(vals_hst_ct[Lim_bines:] < (PORC_HIST_LIM_CT*vals_hst_ct[Lim_bines:].max()))[0]][0]
        bines_ct = np.linspace(LIM_INF_CT,LIM_SUP_CT,128)
        vals_hst_ct, _ = np.histogram(ct_output.flatten(), bins=bines_ct)

        # Get PET histogram
        # Analize only those voxels where there is body mass
        # First get a binary mask of the body
        binary_mask = np.zeros(ct_output.shape)
        binary_mask[ct_output > LIM_INF_CT] = 1
        # Dilate binary mask
        binary_mask = ndimage.morphology.binary_dilation(binary_mask,structure=morphology.ball(3))
        # Get voxels to be analized
        pet_vector_hist = pet_output[binary_mask > 0].flatten()
        # get histogram
        bines_pet = np.linspace(0,pet_vector_hist.max(),4096*8)
        vals_hst_pet, _ = np.histogram(pet_vector_hist, bins=bines_pet)
        # Look for most repeated value
        bin_max_rep_pet = vals_hst_pet[Lim_bines_test:].max()
        # Select a upper limit abobe this point
        pet_thrs_lim = bines_pet[Lim_bines+np.argmax(vals_hst_pet[Lim_bines:])+np.argwhere(vals_hst_pet[Lim_bines+np.argmax(vals_hst_pet[Lim_bines:]):] < bin_max_rep_pet*PORC_HIST_LIM_PET)[0]][0]
        # Create a new histogram
        bines_pet = np.linspace(bines_pet[Lim_bines],pet_thrs_lim,128, endpoint=False)
        vals_hst_pet, _ = np.histogram(pet_vector_hist, bins=bines_pet)
        # Generate datapoint
        pet_low_lim = pet_thrs_lim

        # Generate input and labels
        (input_net, labels_out) = DH.preprocess_sample(ct_output, 
                                                          pet_output, 
                                                          pet_low_lim, 
                                                          LIM_INF_CT, 
                                                          HU_Lung_upper_limit, 
                                                          HU_Parenchyma_lower_limit,
                                                          HU_Parenchyma_upper_limit, 
                                                          HU_Bone_low_limit, 
                                                          HU_Bone_upper_limit,
                                                          HU_Fluids_Fat_lower_limit,
                                                          HU_Fluids_Fat_upper_limit,
                                                          USE_FAT_THRESHOLD = True,
                                                          USE_HIST_THRESHOLD = False)

        # Clean labels from interpolation errors
        labels_out[labels_out <= 0.3] = 0
        labels_out[labels_out > 0.3] = 1
        # Remove double labeling of pixels whith priority order Bones -> parenchyma -> fuid-fat -> Air
        labels_out[:,:,:,2][(labels_out[:,:,:,0]+labels_out[:,:,:,1]) > 0] = 0
        labels_out[:,:,:,1][labels_out[:,:,:,0] > 0] = 0
        labels_out[:,:,:,3][labels_out[:,:,:,0] > 0] = 0
        labels_out[:,:,:,3][labels_out[:,:,:,1] > 0] = 0
        labels_out[:,:,:,3][labels_out[:,:,:,2] > 0] = 0

        # Get mean bone, parenchyma, fluid-fat and air mean PET activity
        PET_bone_mean = np.mean(pet_output[labels_out[:,:,:,0]==1])
        PET_parenchyma_mean = np.mean(pet_output[labels_out[:,:,:,1]==1])
        PET_fluid_fat_mean = np.mean(pet_output[labels_out[:,:,:,2]==1])
        PET_air_mean = np.mean(pet_output[labels_out[:,:,:,3]==1])

        # Clamp PET at double of bone activity
        clamp_mult = 10.0
        input_net[input_net > (clamp_mult*PET_bone_mean)]  = clamp_mult*PET_bone_mean

        # Clamp CT volume
        ct_out_scale = np.copy(ct_output)
        ct_out_scale[ct_out_scale < LIM_INF_CT]  = LIM_INF_CT
        ct_out_scale[ct_out_scale > LIM_SUP_CT]  = LIM_SUP_CT

        # Normalize to range 0:1 
        input_net = (input_net - input_net.min()) / (input_net.max() - input_net.min())
        ct_out_scale = (ct_out_scale - ct_out_scale.min()) / (ct_out_scale.max() - ct_out_scale.min())

        # Create the Example
        example = tf.train.Example(
            features=tf.train.Features(feature={'input': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(input_net, np.prod(input_net.shape)) )),
                                                'label': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(labels_out, np.prod(labels_out.shape)) )),
                                                'ct': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(ct_out_scale, np.prod(ct_out_scale.shape)) )),
                                                'limits_Z': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(limits_Z, np.prod(limits_Z.shape)) ))
                                                }))


        if study_df.Partition.values[0] == "Train":
            write_train_list[0].write(example.SerializeToString())
            write_train_list[0].flush()
            OUTPUT_NAME = OUTPUT_NAME_TRAIN
        elif study_df.Partition.values[0] == "Validation":
            write_validation_list[0].write(example.SerializeToString())
            write_validation_list[0].flush()
            OUTPUT_NAME = OUTPUT_NAME_VALIDATION
        else:
            raise ValueError('Study not in Train nor Validation list!')



        # Create low resolution Example
        for idx_downSample in range(1, DOWNSAMPLE_SCALES):
            low_res_post_scaling = 2**idx_downSample

            input_net_low_res = ndimage.interpolation.zoom(input=input_net[:,:,:], zoom=[1.0/low_res_post_scaling,1.0/low_res_post_scaling,1.0/low_res_post_scaling], order=order_zoom_use)
            ct_out_scale_low_res = ndimage.interpolation.zoom(input=ct_out_scale[:,:,:], zoom=[1.0/low_res_post_scaling,1.0/low_res_post_scaling,1.0/low_res_post_scaling], order=order_zoom_use)
            labels_out_low_res = np.zeros((int(labels_out.shape[0]/low_res_post_scaling),int(labels_out.shape[1]/low_res_post_scaling),int(labels_out.shape[2]/low_res_post_scaling),labels_out.shape[-1]))            
            for idx_label in range(labels_out.shape[-1]):
                labels_out_low_res[:,:,:,idx_label] = ndimage.interpolation.zoom(input=labels_out[:,:,:,idx_label], zoom=[1.0/low_res_post_scaling,1.0/low_res_post_scaling,1.0/low_res_post_scaling], order=order_zoom_use)
            limits_Z_low_res = np.array([int(limits_Z[0]/low_res_post_scaling),int(limits_Z[1]/low_res_post_scaling)])

            # Clean low resolution labels from interpolation errors
            labels_out_low_res[labels_out_low_res<0.5] = 0
            labels_out_low_res[labels_out_low_res>=0.5] = 1
            # Remove double labeling of pixels whith priority order Bones -> parenchyma -> fuid-fat -> Air
            labels_out[:,:,:,2][(labels_out[:,:,:,0]+labels_out[:,:,:,1]) > 0] = 0
            labels_out[:,:,:,1][labels_out[:,:,:,0] > 0] = 0
            labels_out[:,:,:,3][labels_out[:,:,:,0] > 0] = 0
            labels_out[:,:,:,3][labels_out[:,:,:,1] > 0] = 0
            labels_out[:,:,:,3][labels_out[:,:,:,2] > 0] = 0

            example_low_res = tf.train.Example(
                features=tf.train.Features(feature={'input': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(input_net_low_res, np.prod(input_net_low_res.shape)) )),
                                                    'label': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(labels_out_low_res, np.prod(labels_out_low_res.shape)) )),
                                                    'ct': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(ct_out_scale_low_res, np.prod(ct_out_scale_low_res.shape)) )),
                                                    'limits_Z': tf.train.Feature(float_list=tf.train.FloatList(value= np.reshape(limits_Z_low_res, np.prod(limits_Z_low_res.shape)) ))
                                                    }))



            if study_df.Partition.values[0] == "Train":
                write_train_list[idx_downSample].write(example_low_res.SerializeToString())
                write_train_list[idx_downSample].flush()
            elif study_df.Partition.values[0] == "Validation":
                write_validation_list[idx_downSample].write(example_low_res.SerializeToString())
                write_validation_list[idx_downSample].flush()
            else:
                raise ValueError('Study not in Train nor Validation list!')
                exit()






        # Add to dataframe
        DATASET_frame.loc[len(DATASET_frame)] = [patient, 
                                                 OUTPUT_NAME,
                                                 ct_output.min(),
                                                 ct_output.max(),
                                                 ct_low_lim,
                                                 pet_output.min(),
                                                 pet_output.max(),
                                                 pet_low_lim,
                                                 PET_bone_mean,
                                                 PET_parenchyma_mean,
                                                 PET_fluid_fat_mean,
                                                 PET_air_mean,
                                                Objective_size[0],
                                                Objective_size[1],
                                                Objective_size[2]]        

        print("\tLimit: \n\t\tCT:%0.5f (%d %%)\n\t\tPET:%e (%d %%)"%(ct_low_lim,PORC_HIST_LIM_CT*100,
                                                                        pet_low_lim,PORC_HIST_LIM_PET*100,))



        # Print all figures            
        plt.figure(dpi = 50)
        ax = plt.subplot(1, 2, 1)
        ax.stem(bines_ct[Lim_bines:-1], vals_hst_ct[Lim_bines:])
        ax = plt.subplot(1, 2, 2)
        ax.stem(bines_pet[Lim_bines:-1], vals_hst_pet[Lim_bines:])
        plt.tight_layout()
        plt.show()



        print('Max-Min input = %f ; %f'%(input_net.max(),input_net.min()))
        print('Max-Min ct = %f ; %f'%(ct_out_scale.max(),ct_out_scale.min()))
        print('Limits Min-Max Z = %f ; %f'%(limits_Z[0],limits_Z[1]))
        print('Pet means = %e ; %e ; %e ; %e'%(PET_bone_mean,PET_parenchyma_mean,PET_fluid_fat_mean,PET_air_mean))

        X_ct = int(ct_output.shape[0]/2)
        Y_ct = int(ct_output.shape[1]/2)
        Z_ct = int(ct_output.shape[2]/2)

        X_pet = int(pet_output.shape[0]/2)
        Y_pet = int(pet_output.shape[1]/2)
        Z_pet = int(pet_output.shape[2]/2)

        print("\tCT couch removal test:")
        plt.figure(dpi=100)
        ax = plt.subplot(2,3,1)
        ax.imshow(print_ct_orig)
        ax = plt.subplot(2,3,2)
        ax.imshow(couch_mask)
        ax = plt.subplot(2,3,3)
        ax.imshow(print_ct_final)
        plt.tight_layout()
        plt.show()

        print("\tOriginal volumes:")
        plt.figure(dpi=100)
        ax = plt.subplot(3,2,1)
        ax.imshow(ct_output[:,:,Z_ct])
        ax = plt.subplot(3,2,2)
        ax.imshow(pet_output[:,:,Z_pet])
        ax = plt.subplot(3,2,3)
        ax.imshow(ct_output[:,Y_ct,:])
        ax = plt.subplot(3,2,4)
        ax.imshow(pet_output[:,Y_pet,:])
        ax = plt.subplot(3,2,5)
        ax.imshow(ct_output[X_ct,:,:])
        ax = plt.subplot(3,2,6)
        ax.imshow(pet_output[X_pet,:,:])
        plt.tight_layout()
        plt.show()


        X_ct = int(input_net.shape[0]/2)
        Y_ct = int(input_net.shape[1]/2)
        Z_ct = int(input_net.shape[2]/2)

        X_pet = int(input_net.shape[0]/2)
        Y_pet = int(input_net.shape[1]/2)
        Z_pet = int(input_net.shape[2]/2)


        print("\t Network Input and CT Volume:")
        plt.figure(dpi=100)
        ax = plt.subplot(3,2,1)
        ax.imshow(input_net[:,:,Z_pet])
        ax = plt.subplot(3,2,2)
        ax.imshow(ct_out_scale[:,:,Z_ct])
        ax = plt.subplot(3,2,3)
        ax.imshow(input_net[:,Y_pet,:])
        ax = plt.subplot(3,2,4)
        ax.imshow(ct_out_scale[:,X_ct,:])
        ax = plt.subplot(3,2,5)
        ax.imshow(input_net[X_pet,:,:])
        ax = plt.subplot(3,2,6)
        ax.imshow(ct_out_scale[X_ct,:,:])
        plt.tight_layout()
        plt.show()


        print("\tLabels Bone and Parenchyma (except lung):")
        plt.figure(dpi=100)
        ax = plt.subplot(3,2,1)
        ax.imshow(labels_out[:,:,Z_ct,0])
        ax = plt.subplot(3,2,2)
        ax.imshow(labels_out[:,:,Z_pet,1])
        ax = plt.subplot(3,2,3)
        ax.imshow(labels_out[:,Y_ct,:,0])
        ax = plt.subplot(3,2,4)
        ax.imshow(labels_out[:,Y_pet,:,1])
        ax = plt.subplot(3,2,5)
        ax.imshow(labels_out[X_ct,:,:,0])
        ax = plt.subplot(3,2,6)
        ax.imshow(labels_out[X_pet,:,:,1])
        plt.tight_layout()
        plt.show()

        print("\tLabel Fluid-Fat and Air:")
        plt.figure(dpi=100)
        ax = plt.subplot(3,2,1)
        ax.imshow(labels_out[:,:,Z_pet,2])
        ax = plt.subplot(3,2,2)
        ax.imshow(labels_out[:,:,Z_pet,3])
        ax = plt.subplot(3,2,3)
        ax.imshow(labels_out[:,Y_pet,:,2])
        ax = plt.subplot(3,2,4)
        ax.imshow(labels_out[:,Y_pet,:,3])
        ax = plt.subplot(3,2,5)
        ax.imshow(labels_out[X_pet,:,:,2])
        ax = plt.subplot(3,2,6)
        ax.imshow(labels_out[X_pet,:,:,3])
        plt.tight_layout()
        plt.show()





        del(ct_output)
        del(pet_output)

        print('\tProcess time: %f'%(time.time()-tac))
        study_number = study_number + 1


elap_min = ((time.time()-tic)/60.0)
print('Total elapsed time: %d hs %d min'%((elap_min/60.0),(((elap_min/60.0)-(elap_min//60.0))*60.0)))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11301202965090927, 'fraction_numerical': 0.01946802780413984, 'mean_word_length': 1.9700617048830655, 'pattern_counts': {'":': 4, '<': 37, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 177, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '262902', 'license': 'MIT', 'n_tokens_mistral': 13996, 'n_tokens_neox': 12978, 'n_words': 1860, 'path': 'datasets/Generate_train_and_validation_datasets.ipynb', 'repo_name': 'Umizio/PET_DeepAttCorrection', 'size': '38987'}