Volumes of interest¶

In this tutorial we will look at how to represent volumes of interest in datasets.

In [1]:

import matplotlib.pyplot as plt

import hoa_tools.dataset
import hoa_tools.voi

import matplotlib.pyplot as plt import hoa_tools.dataset import hoa_tools.voi

First we'll get an overview dataset, of the brain of donor S-20-29.

In [2]:

overview_dataset = hoa_tools.dataset.get_dataset(
    "S-20-29_brain_complete-organ_25.33um_bm05"
)

overview_dataset = hoa_tools.dataset.get_dataset( "S-20-29_brain_complete-organ_25.33um_bm05" )

Now we'll define a volume of interest (VOI) object. This represents a cuboid shape volume within the full dataset.

In [3]:

voi = hoa_tools.voi.VOI(
    dataset=overview_dataset,
    downsample_level=4,
    lower_corner={"x": 100, "y": 50, "z": 23},
    size={"x": 40, "y": 20, "z": 1},
)

voi = hoa_tools.voi.VOI( dataset=overview_dataset, downsample_level=4, lower_corner={"x": 100, "y": 50, "z": 23}, size={"x": 40, "y": 20, "z": 1}, )

Once we have a VOI object, there's a number of useful things we can do. Lets get the data array for the VOI and plot it:

In [4]:

data_array = voi.get_data_array()
data_array

data_array = voi.get_data_array() data_array

Out[4]:

<xarray.DataArray 'S-20-29_brain_complete-organ_25.33um_bm05' (z: 1, y: 20,
                                                               x: 40)> Size: 2kB
dask.array<getitem, shape=(1, 20, 40), dtype=uint16, chunksize=(1, 20, 28), chunktype=numpy.ndarray>
Coordinates:
  * z        (z) float64 8B 9.321e+03
  * y        (y) float64 160B 2.026e+04 2.067e+04 ... 2.756e+04 2.796e+04
  * x        (x) float64 320B 4.053e+04 4.093e+04 ... 5.593e+04 5.633e+04

xarray.DataArray

'S-20-29_brain_complete-organ_25.33um_bm05'

• z: 1
• y: 20
• x: 40

• dask.array<chunksize=(1, 20, 28), meta=np.ndarray>

  Array Chunk Bytes 1.56 kiB 1.09 kiB Shape (1, 20, 40) (1, 20, 28) Dask graph 2 chunks in 3 graph layers Data type uint16 numpy.ndarray

• Coordinates: (3)
  • z
    (z)
    float64
    9.321e+03

    units :

    μm

    array([9321.44])

  • y
    (y)
    float64
    2.026e+04 2.067e+04 ... 2.796e+04

    units :

    μm

    array([20264.  , 20669.28, 21074.56, 21479.84, 21885.12, 22290.4 , 22695.68,
           23100.96, 23506.24, 23911.52, 24316.8 , 24722.08, 25127.36, 25532.64,
           25937.92, 26343.2 , 26748.48, 27153.76, 27559.04, 27964.32])

  • x
    (x)
    float64
    4.053e+04 4.093e+04 ... 5.633e+04

    units :

    μm

    array([40528.  , 40933.28, 41338.56, 41743.84, 42149.12, 42554.4 , 42959.68,
           43364.96, 43770.24, 44175.52, 44580.8 , 44986.08, 45391.36, 45796.64,
           46201.92, 46607.2 , 47012.48, 47417.76, 47823.04, 48228.32, 48633.6 ,
           49038.88, 49444.16, 49849.44, 50254.72, 50660.  , 51065.28, 51470.56,
           51875.84, 52281.12, 52686.4 , 53091.68, 53496.96, 53902.24, 54307.52,
           54712.8 , 55118.08, 55523.36, 55928.64, 56333.92])

In [5]:

plt.figure()
ax = data_array.plot(cmap="Grays_r")
ax.axes.set_aspect("equal")
ax.axes.set_title("Downsample level = 4")

plt.figure() ax = data_array.plot(cmap="Grays_r") ax.axes.set_aspect("equal") ax.axes.set_title("Downsample level = 4")

Out[5]:

Text(0.5, 1.0, 'Downsample level = 4')

Transforming between dowsample levels¶

VOIs can be transformed to the different resolution levels of the same dataset. Lets get a higher resolution version of the data we plotted above:

In [6]:

high_res_voi = voi.change_downsample_level(new_downsample_level=2)

high_res_voi = voi.change_downsample_level(new_downsample_level=2)

Because we've increased the resolution (equivalently, decreased the downsampling level) the VOI is now bigger.

In [7]:

high_res_data_array = high_res_voi.get_data_array()

high_res_data_array = high_res_voi.get_data_array()

In [8]:

fig, axs = plt.subplots(nrows=2)
data_array.isel(z=0).plot(cmap="Grays_r", ax=axs[0])
axs[0].set_aspect("equal")
axs[0].set_title("Downsample level = 4")

high_res_data_array.isel(z=0).plot(cmap="Grays_r", ax=axs[1])
axs[1].set_aspect("equal")
axs[1].set_title("Downsample level = 2")

fig.tight_layout()

fig, axs = plt.subplots(nrows=2) data_array.isel(z=0).plot(cmap="Grays_r", ax=axs[0]) axs[0].set_aspect("equal") axs[0].set_title("Downsample level = 4") high_res_data_array.isel(z=0).plot(cmap="Grays_r", ax=axs[1]) axs[1].set_aspect("equal") axs[1].set_title("Downsample level = 2") fig.tight_layout()

Transforming between registered datasets¶

Now we'll step through transforming a VOI between two datasets that have been registered to each other. We'll start by getting one of the children of the overview dataset used above.

In [9]:

children = overview_dataset.get_children()
child = children[8]
child.name

children = overview_dataset.get_children() child = children[8] child.name

Out[9]:

'S-20-29_brain_VOI-04_6.5um_bm05'

In [10]:

child_array = child.data_array(downsample_level=0)
child_array

child_array = child.data_array(downsample_level=0) child_array

Out[10]:

<xarray.DataArray 'S-20-29_brain_VOI-04_6.5um_bm05' (z: 4176, y: 3844, x: 3844)> Size: 123GB
dask.array<array, shape=(4176, 3844, 3844), dtype=uint16, chunksize=(128, 128, 128), chunktype=numpy.ndarray>
Coordinates:
  * z        (z) float64 33kB 0.0 6.5 13.0 ... 2.712e+04 2.713e+04 2.714e+04
  * y        (y) float64 31kB 0.0 6.5 13.0 ... 2.497e+04 2.497e+04 2.498e+04
  * x        (x) float64 31kB 0.0 6.5 13.0 ... 2.497e+04 2.497e+04 2.498e+04

xarray.DataArray

'S-20-29_brain_VOI-04_6.5um_bm05'

• z: 4176
• y: 3844
• x: 3844

• dask.array<chunksize=(128, 128, 128), meta=np.ndarray>

  Array Chunk Bytes 114.94 GiB 4.00 MiB Shape (4176, 3844, 3844) (128, 128, 128) Dask graph 31713 chunks in 2 graph layers Data type uint16 numpy.ndarray

• Coordinates: (3)
  • z
    (z)
    float64
    0.0 6.5 ... 2.713e+04 2.714e+04

    units :

    μm

    array([0.00000e+00, 6.50000e+00, 1.30000e+01, ..., 2.71245e+04, 2.71310e+04,
           2.71375e+04], shape=(4176,))

  • y
    (y)
    float64
    0.0 6.5 ... 2.497e+04 2.498e+04

    units :

    μm

    array([0.00000e+00, 6.50000e+00, 1.30000e+01, ..., 2.49665e+04, 2.49730e+04,
           2.49795e+04], shape=(3844,))

  • x
    (x)
    float64
    0.0 6.5 ... 2.497e+04 2.498e+04

    units :

    μm

    array([0.00000e+00, 6.50000e+00, 1.30000e+01, ..., 2.49665e+04, 2.49730e+04,
           2.49795e+04], shape=(3844,))

Now we'll select a smaller volume of interest from the child array, and plot it

In [11]:

child_voi = hoa_tools.voi.VOI(
    dataset=child,
    downsample_level=0,
    lower_corner={"x": 3434, "y": 2060, "z": 2656},
    size={"x": 256, "y": 256, "z": 128},
)
child_array = child_voi.get_data_array()

child_voi = hoa_tools.voi.VOI( dataset=child, downsample_level=0, lower_corner={"x": 3434, "y": 2060, "z": 2656}, size={"x": 256, "y": 256, "z": 128}, ) child_array = child_voi.get_data_array()

In [12]:

plt.figure()
ax = child_array.isel(z=0).plot(cmap="Grays_r")
ax.axes.set_aspect("equal")
ax.axes.set_title("Zoom dataset slice")

plt.figure() ax = child_array.isel(z=0).plot(cmap="Grays_r") ax.axes.set_aspect("equal") ax.axes.set_title("Zoom dataset slice")

Out[12]:

Text(0.5, 1.0, 'Zoom dataset slice')

Now lets transform the VOI to the overview dataset

In [13]:

overview_voi = child_voi.transform_to(overview_dataset)
overview_voi

overview_voi = child_voi.transform_to(overview_dataset) overview_voi

Out[13]:

VOI(dataset=Dataset(name='S-20-29_brain_complete-organ_25.33um_bm05', dataset_type='overview', voi='complete-organ', data=Data(shape=[5750, 5750, 6627], voxel_size_um=25.33, gcs_url='n5://gs://ucl-hip-ct-35a68e99feaae8932b1d44da0358940b/S-20-29/brain/25.33um_complete-organ_bm05/'), sample=Sample(organ='brain', organ_context=None, post_mortem_interval_hours=None, fixation_method=None, fixation_medium='formalin', organ_infilled=None, stabilisation_medium='crushed agar', degassing_method='vacuum', scan_solvent='ethanol', scan_solvent_concentration=70.0, scan_temperature='room temperature'), donor=Donor(id='S-20-29', age=Age(root=80), sex='M', weight=Weight(root=90.0), height=Height(root=170.0), cause_of_death='Acute right heart failure with extensive Pulmonary artery thrombembolism, diffuse alveolar damage, shock kidneys bilateral with acute kidney failure, and acute liver failure COVID-19 disease', medical_history='Kidney failure, arterial hypertension, Z.n. Spondylodesis LWK 2 and sacral vertebra 2, SARS-CoV2 infection, symptoms of memory loss and predementia, Dandy-Walker syndrome variant, post-mortem diagnostic.', diabetes='Type 2', hypertension=None, smoker=None), scan=Scan(date=datetime.date(2021, 3, 4), beamline='BM05', energy=Energy(root=80.0), current_start=None, filling_mode=None, n_frames=NFrames(root=9990), n_ref=NRef(root=3), n_dark=NDark(root=0), latency_time=LatencyTime(root=0.0025), exposure_time=ExposureTime(root=0.0025), subframe_time=None, n_subframes=NSubframes(root=5), scan_type='zseries', scan_range=360.0, n_scans=None, acquisition='quarter', z_step=None, scan_time=None, filters=['Mo 0.24mm', 'SiO2 8mm block'], scintillator=None, optic='DZoom', optic_magnification=None, distance_source_sample=None, distance_sample_detector=DistanceSampleDetector(root=3500.0), psho=None, sensor_name='PCO edge 4.2 CLHS', sensor_mode='FFM', sensor_roi_x_size=SensorRoiXSize(root=3100), sensor_roi_y_size=SensorRoiYSize(root=6629), sensor_binning=None, pixel_size=25.33, xray_magnification=None, technique='Hierarchical Phase-Contrast Tomography (HiP-CT)', experiment_type='Tomography'), proposal=Proposal(proposal_number='md1252', title="Multiscale Quantification of Covid-19's impact on lung vasculature from whole lobe to alveolar/microvascular scales", proposers=['Peter Lee', 'Paul Tafforeau', 'Danny Jonigk', 'Maximilian Ackermann', 'Mark P. Kuhnel', 'Elodie Boller', 'Willi Wagner']), registration=None, citation=Citation(title='Overview at 25.33um of the brain of donor S-20-29, scanned at ESRF on beamline BM05.', contributors=[Contributor(first_name='The', last_name='Human Organ Atlas Collaboration', orcid=None, roles=[]), Contributor(first_name='Claire', last_name='Walsh', orcid=Orcid(root='0000-0003-3769-3392'), roles=['Conceptualization', 'Funding acquisition']), Contributor(first_name='Danny', last_name='Jonigk', orcid=Orcid(root='0000-0002-5251-2281'), roles=['Funding acquisition']), Contributor(first_name='David', last_name='Stansby', orcid=Orcid(root='0000-0002-1365-1908'), roles=['Data curation', 'Software']), Contributor(first_name='Elodie', last_name='Boller', orcid=None, roles=['Funding acquisition']), Contributor(first_name='Guillaume', last_name='Gaisne', orcid=Orcid(root='0000-0002-3401-7930'), roles=['Data curation', 'Software']), Contributor(first_name='Hector', last_name='Dejea I Velardo', orcid=Orcid(root='0000-0003-2584-9812'), roles=['Data curation', 'Investigation']), Contributor(first_name='Joanna', last_name='Purzycka', orcid=Orcid(root='0000-0001-8742-5279'), roles=['Data curation']), Contributor(first_name='Joseph', last_name='Brunet', orcid=Orcid(root='0000-0002-8424-9510'), roles=['Data curation', 'Investigation', 'Software']), Contributor(first_name='Mark', last_name='P. Kuhnel', orcid=Orcid(root='0000-0003-3558-2576'), roles=['Funding acquisition']), Contributor(first_name='Maximilian', last_name='Ackermann', orcid=Orcid(root='0000-0001-9996-2477'), roles=['Funding acquisition']), Contributor(first_name='Paul', last_name='Tafforeau', orcid=Orcid(root='0000-0002-5962-1683'), roles=['Conceptualization', 'Data curation', 'Funding acquisition']), Contributor(first_name='Peter', last_name='Lee', orcid=Orcid(root='0000-0002-3898-8881'), roles=['Conceptualization', 'Funding acquisition']), Contributor(first_name='Theresa', last_name='Urban', orcid=Orcid(root='0000-0002-0352-8180'), roles=['Data curation', 'Investigation', 'Software']), Contributor(first_name='Willi', last_name='Wagner', orcid=None, roles=['Funding acquisition'])], doi='10.15151/ESRF-DC-1773961413', author_list=['The Human Organ Atlas Collaboration', 'Claire Walsh', 'Danny Jonigk', 'David Stansby', 'Elodie Boller', 'Guillaume Gaisne', 'Hector Dejea I Velardo', 'Joanna Purzycka', 'Joseph Brunet', 'Mark P. Kuhnel', 'Maximilian Ackermann', 'Paul Tafforeau', 'Peter Lee', 'Theresa Urban', 'Willi Wagner']), metadata_version='1.1'), downsample_level=0, lower_corner=ArrayCoordinate(x=2975, y=1689, z=4316), size=ArrayCoordinate(x=69, y=69, z=33))

Finally, lets plot the overview and child datasets together. We can see that the same structures are visible in both, so the registration is approximately accurate, and all the transformations worked.

In the next tutorial section we'll see how to do a more accurate registration of just these small sub-volumes of data in order to do quantitative comparisons.

In [14]:

overview_array = overview_voi.get_data_array()
overview_array

overview_array = overview_voi.get_data_array() overview_array

Out[14]:

<xarray.DataArray 'S-20-29_brain_complete-organ_25.33um_bm05' (z: 33, y: 69,
                                                               x: 69)> Size: 314kB
dask.array<getitem, shape=(33, 69, 69), dtype=uint16, chunksize=(33, 69, 69), chunktype=numpy.ndarray>
Coordinates:
  * z        (z) float64 264B 1.093e+05 1.093e+05 ... 1.101e+05 1.101e+05
  * y        (y) float64 552B 4.278e+04 4.281e+04 ... 4.448e+04 4.45e+04
  * x        (x) float64 552B 7.536e+04 7.538e+04 ... 7.705e+04 7.708e+04

xarray.DataArray

'S-20-29_brain_complete-organ_25.33um_bm05'

• z: 33
• y: 69
• x: 69

• dask.array<chunksize=(33, 69, 69), meta=np.ndarray>

  Array Chunk Bytes 306.86 kiB 306.86 kiB Shape (33, 69, 69) (33, 69, 69) Dask graph 1 chunks in 3 graph layers Data type uint16 numpy.ndarray

• Coordinates: (3)
  • z
    (z)
    float64
    1.093e+05 1.093e+05 ... 1.101e+05

    units :

    μm

    array([109324.28, 109349.61, 109374.94, 109400.27, 109425.6 , 109450.93,
           109476.26, 109501.59, 109526.92, 109552.25, 109577.58, 109602.91,
           109628.24, 109653.57, 109678.9 , 109704.23, 109729.56, 109754.89,
           109780.22, 109805.55, 109830.88, 109856.21, 109881.54, 109906.87,
           109932.2 , 109957.53, 109982.86, 110008.19, 110033.52, 110058.85,
           110084.18, 110109.51, 110134.84])

  • y
    (y)
    float64
    4.278e+04 4.281e+04 ... 4.45e+04

    units :

    μm

    array([42782.37, 42807.7 , 42833.03, 42858.36, 42883.69, 42909.02, 42934.35,
           42959.68, 42985.01, 43010.34, 43035.67, 43061.  , 43086.33, 43111.66,
           43136.99, 43162.32, 43187.65, 43212.98, 43238.31, 43263.64, 43288.97,
           43314.3 , 43339.63, 43364.96, 43390.29, 43415.62, 43440.95, 43466.28,
           43491.61, 43516.94, 43542.27, 43567.6 , 43592.93, 43618.26, 43643.59,
           43668.92, 43694.25, 43719.58, 43744.91, 43770.24, 43795.57, 43820.9 ,
           43846.23, 43871.56, 43896.89, 43922.22, 43947.55, 43972.88, 43998.21,
           44023.54, 44048.87, 44074.2 , 44099.53, 44124.86, 44150.19, 44175.52,
           44200.85, 44226.18, 44251.51, 44276.84, 44302.17, 44327.5 , 44352.83,
           44378.16, 44403.49, 44428.82, 44454.15, 44479.48, 44504.81])

  • x
    (x)
    float64
    7.536e+04 7.538e+04 ... 7.708e+04

    units :

    μm

    array([75356.75, 75382.08, 75407.41, 75432.74, 75458.07, 75483.4 , 75508.73,
           75534.06, 75559.39, 75584.72, 75610.05, 75635.38, 75660.71, 75686.04,
           75711.37, 75736.7 , 75762.03, 75787.36, 75812.69, 75838.02, 75863.35,
           75888.68, 75914.01, 75939.34, 75964.67, 75990.  , 76015.33, 76040.66,
           76065.99, 76091.32, 76116.65, 76141.98, 76167.31, 76192.64, 76217.97,
           76243.3 , 76268.63, 76293.96, 76319.29, 76344.62, 76369.95, 76395.28,
           76420.61, 76445.94, 76471.27, 76496.6 , 76521.93, 76547.26, 76572.59,
           76597.92, 76623.25, 76648.58, 76673.91, 76699.24, 76724.57, 76749.9 ,
           76775.23, 76800.56, 76825.89, 76851.22, 76876.55, 76901.88, 76927.21,
           76952.54, 76977.87, 77003.2 , 77028.53, 77053.86, 77079.19])

In [15]:

fig, axs = plt.subplots(nrows=2, figsize=(6, 9))

ax = axs[0]
child_array.isel(z=0).plot(cmap="Grays_r", ax=ax)
ax.set_aspect("equal")
ax.set_title("Zoom dataset")

ax = axs[1]
overview_array.isel(z=0).plot(cmap="Grays_r", ax=ax)
ax.set_aspect("equal")
ax.set_title("Overview dataset")
fig.tight_layout()

fig, axs = plt.subplots(nrows=2, figsize=(6, 9)) ax = axs[0] child_array.isel(z=0).plot(cmap="Grays_r", ax=ax) ax.set_aspect("equal") ax.set_title("Zoom dataset") ax = axs[1] overview_array.isel(z=0).plot(cmap="Grays_r", ax=ax) ax.set_aspect("equal") ax.set_title("Overview dataset") fig.tight_layout()