def load_reconstruction(folder, slice_start=0, slice_end=-1):
    """Load a volume from folder, also returns the corresponding partition.

    Parameters
    ----------
    folder : str
        Path to the folder where the DICOM files are stored.
    slice_start : int
        Index of the first slice to use. Used for subsampling.
    slice_end : int
        Index of the final slice to use.

    Returns
    -------
    partition : `odl.RectPartition`
        Partition describing the geometric positioning of the voxels.
    data : `numpy.ndarray`
        Volumetric data. Scaled such that data = 1 for water (0 HU).

    Notes
    -----
    DICOM data is highly non trivial. Typically, each slice has been computed
    with a slice tickness (e.g. 3mm) but the slice spacing might be
    different from that.

    Further, the coordinates in DICOM is typically the *middle* of the pixel,
    not the corners as in ODL.

    This function should handle all of these peculiarities and give a volume
    with the correct coordinate system attached.
    """
    file_names = sorted([f for f in os.listdir(folder) if f.endswith(".IMA")])

    if len(file_names) == 0:
        raise ValueError('No DICOM files found in {}'.format(folder))

    volumes = []
    datasets = []

    file_names = file_names[slice_start:slice_end]

    for file_name in tqdm.tqdm(file_names, 'loading volume data'):
        # read the file
        dataset = dicom.read_file(folder + '/' + file_name)

        # Get parameters
        pixel_size = np.array(dataset.PixelSpacing)
        pixel_thickness = float(dataset.SliceThickness)
        rows = dataset.Rows
        cols = dataset.Columns

        # Get data array and convert to correct coordinates
        data_array = np.array(np.frombuffer(dataset.PixelData, 'H'),
                              dtype='float32')
        data_array = data_array.reshape([cols, rows], order='C')
        data_array = np.rot90(data_array, -1)

        # Convert from storage type to densities
        # TODO: Optimize these computations
        hu_values = (dataset.RescaleSlope * data_array +
                     dataset.RescaleIntercept)
        densities = (hu_values + 1000) / 1000

        # Store results
        volumes.append(densities)
        datasets.append(dataset)

    voxel_size = np.array(list(pixel_size) + [pixel_thickness])
    shape = np.array([rows, cols, len(volumes)])

    # Compute geometry parameters
    mid_pt = (np.array(dataset.ReconstructionTargetCenterPatient) -
              np.array(dataset.DataCollectionCenterPatient))
    reconstruction_size = (voxel_size * shape)
    min_pt = mid_pt - reconstruction_size / 2
    max_pt = mid_pt + reconstruction_size / 2

    # axis 1 has reversed convention
    min_pt[1], max_pt[1] = -max_pt[1], -min_pt[1]

    if len(datasets) > 1:
        slice_distance = np.abs(
            float(datasets[1].DataCollectionCenterPatient[2]) -
            float(datasets[0].DataCollectionCenterPatient[2]))
    else:
        # If we only have one slice, we must approximate the distance.
        slice_distance = pixel_thickness

    # The middle of the minimum/maximum slice can be computed from the
    # DICOM attribute "DataCollectionCenterPatient". Since ODL uses corner
    # points (e.g. edge of volume) we need to add half a voxel thickness to
    # both sides.
    min_pt[2] = -np.array(datasets[0].DataCollectionCenterPatient)[2]
    min_pt[2] -= 0.5 * slice_distance
    max_pt[2] = -np.array(datasets[-1].DataCollectionCenterPatient)[2]
    max_pt[2] += 0.5 * slice_distance

    partition = odl.uniform_partition(min_pt, max_pt, shape)

    volume = np.transpose(np.array(volumes), (1, 2, 0))

    return partition, volume
def load_reconstruction ( folder , slice_start = 0 , slice_end = - 1 ) : file_names = sorted ( [ f for f in os . listdir ( folder ) if f . endswith ( ".IMA" ) ] ) if len ( file_names ) == 0 : raise ValueError ( 'No DICOM files found in {}' . format ( folder ) ) volumes = [ ] datasets = [ ] file_names = file_names [ slice_start : slice_end ] for file_name in tqdm . tqdm ( file_names , 'loading volume data' ) : dataset = dicom . read_file ( folder + '/' + file_name ) pixel_size = np . array ( dataset . PixelSpacing ) pixel_thickness = float ( dataset . SliceThickness ) rows = dataset . Rows cols = dataset . Columns data_array = np . array ( np . frombuffer ( dataset . PixelData , 'H' ) , dtype = 'float32' ) data_array = data_array . reshape ( [ cols , rows ] , order = 'C' ) data_array = np . rot90 ( data_array , - 1 ) hu_values = ( dataset . RescaleSlope * data_array + dataset . RescaleIntercept ) densities = ( hu_values + 1000 ) / 1000 volumes . append ( densities ) datasets . append ( dataset ) voxel_size = np . array ( list ( pixel_size ) + [ pixel_thickness ] ) shape = np . array ( [ rows , cols , len ( volumes ) ] ) mid_pt = ( np . array ( dataset . ReconstructionTargetCenterPatient ) - np . array ( dataset . DataCollectionCenterPatient ) ) reconstruction_size = ( voxel_size * shape ) min_pt = mid_pt - reconstruction_size / 2 max_pt = mid_pt + reconstruction_size / 2 min_pt [ 1 ] , max_pt [ 1 ] = - max_pt [ 1 ] , - min_pt [ 1 ] if len ( datasets ) > 1 : slice_distance = np . abs ( float ( datasets [ 1 ] . DataCollectionCenterPatient [ 2 ] ) - float ( datasets [ 0 ] . DataCollectionCenterPatient [ 2 ] ) ) else : slice_distance = pixel_thickness min_pt [ 2 ] = - np . array ( datasets [ 0 ] . DataCollectionCenterPatient ) [ 2 ] min_pt [ 2 ] -= 0.5 * slice_distance max_pt [ 2 ] = - np . array ( datasets [ - 1 ] . DataCollectionCenterPatient ) [ 2 ] max_pt [ 2 ] += 0.5 * slice_distance partition = odl . uniform_partition ( min_pt , max_pt , shape ) volume = np . transpose ( np . array ( volumes ) , ( 1 , 2 , 0 ) ) return partition , volume
Load a volume from folder, also returns the corresponding partition.

    Parameters
    ----------
    folder : str
        Path to the folder where the DICOM files are stored.
    slice_start : int
        Index of the first slice to use. Used for subsampling.
    slice_end : int
        Index of the final slice to use.

    Returns
    -------
    partition : `odl.RectPartition`
        Partition describing the geometric positioning of the voxels.
    data : `numpy.ndarray`
        Volumetric data. Scaled such that data = 1 for water (0 HU).

    Notes
    -----
    DICOM data is highly non trivial. Typically, each slice has been computed
    with a slice tickness (e.g. 3mm) but the slice spacing might be
    different from that.

    Further, the coordinates in DICOM is typically the *middle* of the pixel,
    not the corners as in ODL.

    This function should handle all of these peculiarities and give a volume
    with the correct coordinate system attached.