.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/02-parameter-inference/02-lookup-table.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via Binder. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_autoexamples_02-parameter-inference_02-lookup-table.py: ==================== MP2RAGE lookup table ==================== The scope of this notebook is to map T1 from a two-block MP2RAGE, by interpolating along a curve and by matching the same curve, and to sweep the number of points to show which of the two is limited by it. With a single unknown a dictionary degenerates: the atoms lie on a curve rather than filling a space. Interpolating between the two nearest then costs nothing and takes the grid spacing out of the answer, which is what :class:`~torchsim.LookupTable` does. .. GENERATED FROM PYTHON SOURCE LINES 17-21 .. colab-link:: :needs_gpu: 0 !pip install torchsim brainweb-dl cmap .. GENERATED FROM PYTHON SOURCE LINES 23-27 The problem is stated over a simulator carrying the sequence and filled in by an estimator. :func:`~torchsim.execution` decides where that work runs, and the timings below are taken inside it. .. GENERATED FROM PYTHON SOURCE LINES 28-171 .. code-block:: Python import time import numpy as np import torch import torchsim from torchsim.estimators import DictionaryMatcher, LookupTable from torchsim.simulators import MP2RAGESimulator .. GENERATED FROM PYTHON SOURCE LINES 172-181 Phantom ------- BrainWeb subject 0, slice 90: an axial slice at 1 mm through the lateral ventricles. BrainWeb publishes fuzzy memberships rather than labels, so each voxel holds a fraction of each tissue, and the relaxation times are weighted by those fractions. A third of the voxels are mixtures, so the truth is a continuum and not four values. .. GENERATED FROM PYTHON SOURCE LINES 182-218 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_001.png :alt: BrainWeb subject 0, slice 90 :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 0.9901188654526398, 'BrainWeb subject 0, slice 90') .. GENERATED FROM PYTHON SOURCE LINES 219-227 Protocol -------- One inversion, two spoiled gradient-echo blocks read at two inversion times. Each block samples the centre of k-space of its own shot train, so a voxel contributes two numbers. The train spoils after every readout, so T1 is the only tissue property that moves them. .. GENERATED FROM PYTHON SOURCE LINES 228-239 .. code-block:: Python PROTOCOL = dict( TI=(800.0, 2700.0), flip=(4.0, 5.0), TRspgr=6.7, TRmp2rage=6000.0, nshots=128, ) INVERSION_EFFICIENCY = 0.96 simulator = MP2RAGESimulator(**PROTOCOL, inv_efficiency=INVERSION_EFFICIENCY) .. GENERATED FROM PYTHON SOURCE LINES 240-248 Signal curve ------------ Neither block alone says T1: both carry the proton density and the receive gain. The unified combination divides that scale out, leaving a number between -0.5 and 0.5 that depends on T1 alone. Which combination is monotonic belongs to the sequence, so it is given rather than assumed. .. GENERATED FROM PYTHON SOURCE LINES 249-256 .. code-block:: Python def unified(blocks): """The MP2RAGE unified image: scale-free, and a function of T1 alone.""" return (blocks[..., 0] * blocks[..., 1]) / blocks.square().sum(-1).clamp_min(1e-12) .. GENERATED FROM PYTHON SOURCE LINES 257-261 The curve is not monotonic over every T1, and where it turns back it has no inverse. The table keeps the longest monotonic run and reports what it spans, so the invertible range is a number rather than an assumption. .. GENERATED FROM PYTHON SOURCE LINES 262-290 .. code-block:: Python sweep = torch.arange(50.0, 6000.0, 10.0) curve = unified(simulator.simulate(T1=sweep, M0=1.0)) .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_002.png :alt: the two blocks, the curve a T1 is read off :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 291-297 Measurement ----------- Both blocks, at the true T1 and proton density of every brain voxel, with noise at half a percent of the peak magnetization. .. GENERATED FROM PYTHON SOURCE LINES 298-351 .. code-block:: Python clean = simulator.simulate(T1=truth, M0=density) NOISE_STD = float(0.005 * clean.abs().max()) generator = torch.Generator().manual_seed(42) measured = clean + NOISE_STD * torch.randn(clean.shape, generator=generator) .. GENERATED FROM PYTHON SOURCE LINES 352-360 Two estimators -------------- Both are given the same T1 grid. The match scores the two-block signal against every atom and takes the nearest; the table reduces both blocks to the unified number and interpolates along the curve, so ``combine`` is all it is told. Neither is given the range in advance. .. GENERATED FROM PYTHON SOURCE LINES 361-369 .. code-block:: Python grid = torch.linspace(50.0, 6000.0, 60) table = LookupTable(simulator.bind(M0=1.0), combine=unified).fit(T1=grid, seed=0) maps = table.map(measured) # {"T1": ...}, one value per voxel match = DictionaryMatcher(simulator.bind(M0=1.0)).fit(T1=grid, seed=0) .. GENERATED FROM PYTHON SOURCE LINES 370-373 Sweeping the grid separates the method from the sampling. Times are the best of three passes over the slice. .. GENERATED FROM PYTHON SOURCE LINES 374-397 .. rst-class:: sphx-glr-script-out .. code-block:: none points match table match table error error time time ----------------------------------------------- 30 7.93% 0.88% 1.9ms 0.4ms 60 3.29% 0.75% 2.8ms 0.4ms 120 0.78% 0.71% 2.9ms 0.5ms 250 0.80% 0.71% 5.2ms 0.5ms 500 0.62% 0.71% 10.6ms 0.5ms 1000 0.66% 0.71% 21.8ms 0.6ms 2000 0.71% 0.71% 47.7ms 0.6ms .. GENERATED FROM PYTHON SOURCE LINES 398-406 The table is at its floor from the coarsest grid and does not move again. The match starts an order of magnitude worse and climbs to the same place, paying for it in points: its search is one comparison per atom per voxel, where the table's binary search grows with the logarithm. The floor both reach is the noise. A fine enough grid matches a table exactly; the table's advantage is that it was never told how fine. .. GENERATED FROM PYTHON SOURCE LINES 407-435 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_003.png :alt: what the grid costs, what it costs to pay it :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 436-442 Maps ---- At the point count each needs: the table at sixty, the match at a grid fine enough not to limit it. .. GENERATED FROM PYTHON SOURCE LINES 443-461 .. rst-class:: sphx-glr-script-out .. code-block:: none the table keeps 38 of 60 points -- the monotonic run -- and spans unified -0.50 to 0.49 .. GENERATED FROM PYTHON SOURCE LINES 462-466 Neither estimates M0. Both answer with a T1, and the two blocks it predicts are a shape the measurement is a multiple of, so the multiple is one inner product per voxel. .. GENERATED FROM PYTHON SOURCE LINES 467-513 .. code-block:: Python def proton_density(maps): """The scale the measurement is, of the blocks the answer predicts.""" predicted = simulator.simulate(T1=maps["T1"], M0=1.0) return (predicted * measured).sum(-1) / predicted.square().sum(-1).clamp_min(1e-12) M0_map = proton_density(maps) .. rst-class:: sphx-glr-script-out .. code-block:: none method train map model peak T1 M0 ------------------------------------------------------------------------------ lookup, 60 points 0.00s 0.4ms 0.00 MiB -- 0.75% 0.46% match, 2000 atoms 0.00s 46.9ms 0.05 MiB -- 0.71% 0.45% .. GENERATED FROM PYTHON SOURCE LINES 514-569 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_004.png :alt: truth, lookup, match :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_004.png :class: sphx-glr-multi-img * .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_005.png :alt: Δ lookup, Δ match :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_02-lookup-table_005.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.785 seconds) .. _sphx_glr_download_generated_autoexamples_02-parameter-inference_02-lookup-table.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/firmlab-pisa/torchsim/gh-pages?urlpath=lab/tree/v0.0.8/examples/generated/autoexamples/02-parameter-inference/02-lookup-table.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02-lookup-table.ipynb <02-lookup-table.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02-lookup-table.py <02-lookup-table.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02-lookup-table.zip <02-lookup-table.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_