11.7. Estimating Spectrum from Generated 1D Phantom data

The most crucial tool in this software is to estimate the spectrum from the multidimensional MR data. For an unknown spectrum \(\boldsymbol{f}_i\) for \(i\)-th voxel, we solve following nonegatively constrained optimisation problem with spatial smoothness regulariser,

where \(\boldsymbol{\mathrm{d}}_i\) is the vectorised multidimensional MR data acquired from the voxel \(i\) and \(\boldsymbol{\mathrm{f}}\) is the vector formed by concatenating the unknown spectra from all:math:L voxels. \(\boldsymbol{\mathrm{K}}\) is the dictionary matrix that maps the spectral amplitude to the MR data. We have shown here how to create the dictionary tensor. We just flatten the spectral dimension of the dictionary tensor to get the above mentioned matrix form. The novelty of this method is the use of the spatial smoothness regulariser term along with more general non-negativity constraint. This spatial smoothness regulariser minimises the the difference between the neighboring voxels’ spectra. This paper elaborately explains the concecpt behind the regulariser. The smoothness is controlled by the regularization parameter \(\lambda\). If you choose \(\lambda\) too large, you may get solutions that look over-smoothed. If you choose lambda too small, you may get solutions that are too noisy to be useful. This is an engineering trade-off, and it is up to the user to identify what they like. This spatially regularized reconstruction is often hundreds of times better than unregularized reconstruction from an estimation theoretic point of view (cite), although is computationally slower. We provide two algorithms for regularized reconstruction: LADMM (which is often the fastest cite), ADMM (which is generally slower but is still provided for reference cite), and the conventional unregularized NNLS solution (which is very fast but frequently produces worse results). These can be chosen by appropriately setting parameters in the configuration file.

Citation

If you use ADMM algorithm for your spectral estimation, please cite:

Kim et al., Diffusion-relaxation correlation spectroscopic imaging: A multidimensional approach for probing microstructure, 2017. DOI: 10.1002/mrm.26629

Liu et al., An Efficient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping With Applications to Multicomponent Diffusion and Relaxation MRI,2025. DOI: 10.1109/TCI.2025.3609974

If you use LADMM algorithm for your spectral estimation, please cite:

Liu et al., An Efficient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping With Applications to Multicomponent Diffusion and Relaxation MRI,2025. DOI: 10.1109/TCI.2025.3609974

11.7.1. Goal

This tutorial will guide you to estimate the spectrum from the 1D phantom data generated in the last tutorial. Specifically we will be,

Creating the configuration file (as discussed here) to specify algorithm parameters for spectrum estimation.
Choosing optimal penalty parameter for ADMM/LADMM algorithm
Finally estimating spectrum using relevant files.

11.7.2. Get started

Before going to the process of estimating spectra, lets check if we have the folder with the following files (we will be working in the Documents directory):
```
/Documents/Phantom1D/
├── Phantom_data.mat
├── Phantom_spectrum_info.mat
├── Phantom_mask.mat
└── Phantom_mask_beta_calc.mat
```
Once we have these files, we have to create configfile to provide details of the algorithm. We have created the following .ini file using the text editor to use LADMM algorithm and saved it in the ‘Documents’ folder with file name 'Phantom1D_ladmm.ini'.
```
dc_comp = 1
lambda = 0.1

[solver]
name='LADMM'
num_iter=3000
beta=0.1
save_inter=100
check_tol=500
tol=1e-3

[low_rank]
flag = 1
rank = 21
```
This file specifies the parameters of the spectrum estimation tool:

Top-level (global)
- dc_comp: Assigning to 1 we are enabling the estimation algorithm to get rid of any constant offset present in the data. If the data does not seem to have a constant offset value, we can assign ‘dc_comp’ value to 0. Only required for ADMM or LADMM algorithm.
- lambda parameter is the value of spatial regulariser. Higher value of this parameter will result in smoother spectrum. Only required for 'ADMM' or 'LADMM' algorithm. 'NNLS' does not support usage of the spatial regulariser.
[solver]
- name: Solver to run, For this tutorial we chose to use 'LADMM' algorithm. Change the name to choose from other two available algorithms - 'ADMM' , 'NNLS'. 'NNLS' solves the first optimization problem mentioned above, whereas 'ADMM' or 'LADMM' solves the second optimzation problem.
- num_iter: Maximum iterations before stopping (hard cap). Set to 3000 here. If the tolerance criteria as mentioned below is not met, the algorithm will stop after 2000 iteration. Higher number of iteration will make sure the final spectra is a converged solution.
- beta: LADMM penalty parameter, chosen 0.1 arbitrarily. Follow the next step to fine tune this value for faster convergence.
- save_inter: Save intermediate outputs 100 iterations in this case. Output types vary depending on the function used.
- tol: Convergence tolerance on the stopping metric. If the tolerance value of 1e-3 is achieved before num_iter, the spectrum estimation process will stop.
- check_tol: Check the stopping criterion after 500 iteration for skipping the initial transient behavior of cost function. If the tol value is reached before check_tol iterations, the algorithm won’t stop, because the tolerance value is checked after check_tol iterations.
[low_rank]
- flag: Enable/disable low-rank approximation of the matrix, K (1 = on, 0 = off).
- rank: Target rank for the matrix K. Lower the rank, faster the LADMM algorithm. The user should choose the rank by looking at the distribution of the singular values of the tensor, K.
See more details about each entry of this configuration file here.
Although we have provided a beta value in the config file, we are not sure if that is going to give us fastest convergence for the above dataset. So now we will create another text file containing beta values to compare convergence speed. Open a text editor and save the file with name ‘betafile_phantom1D.txt’ in ‘Documents’ folder. Note that the variable name should always be ‘beta_vals’ in the text file.
```
beta_vals= [0.001 0.01 0.1 1 10];
```
Now the we can run the tool that will calculate the convergence speed for the above beta values. Note that our Phantom generation code generates two spatial masks, one for the whole MR volume and another for the penalty parameter selection. We use the second spatial mask named 'Phantom1D/Phantom_mask_beta_calc.mat' in this case. You can also create your own patches of masks following the mask file format. We recommend having multiple 3x3 patches for each slice and choose the region such that you include spatial region having most partial voluming.
SH
plot_beta_sweep.sh -imgfile 'Phantom1D/Phantom_data.mat' -betafile 'betafile_phantom1D.txt' -spect_infofile 'Phantom1D/Phantom_spectrm_info.mat' -configfile 'Phantom1D_ladmm.ini' -spatmaskfile 'Phantom1D/Phantom_mask_beta_calc.mat' -outprefix 'Phantom1D/Phantom1D_data_ladmm_spect' -file_types 'png'
MATLAB
plot_beta_sweep('imgfile','Phantom1D/Phantom_data.mat','betafile','betafile_phantom1D.txt', 'spect_infofile','Phantom1D/Phantom_spectrm_info.mat','configfile','Phantom1D_ladmm.ini','spatmaskfile','Phantom1D/Phantom_mask_beta_calc.mat','outprefix','Phantom1D/Phantom1D_data_ladmm_spect','file_types','png')
This will save the image ‘Phantom1D/Phantom1D_data_ladmm_spect_costVSiter.png’ as shown below. As you can see that the convergence rate is almost same for all the beta values except 0.001 and 0.01. So we don’t need to change the config file, 'Phantom1D_ladmm.ini'.

Note

This plot of cost function vs iteration number for multiple beta can give us an approximate number of maximum iterations required for achieving a converged solution. As you can see from the figure above that, for beta= 0.5, 1, and 5 the cost function is not changing much after 2000 iterations (so, converged!), whereas, the cost function still haven’t converged for beta = 0.05. So if you choose beta value of 0.05 and use num_iter value to 3000, you will not ever get a converged solution. Most of the times, these unconverged solution vary significantly from the converged solution. So its better to start with a higher value of num_iter to generate the above plot. This way we can get an idea about the maximum iteration number, approximately needed for the whole volume spectrum estimation.

Once we have all the required files set, we can run the spectral estimation solver,

SH

estimate_spectra.sh -imgfile  'Phantom1D/Phantom_data.mat' -spect_infofile 'Phantom1D/Phantom_spectrm_info.mat' -configfile 'Phantom1D_ladmm.ini' -spatmaskfile 'Phantom1D/Phantom_mask.mat' -outprefix 'Phantom1D/Phantom1D_data_ladmm_spect'

MATLAB

estimate_spectra('imgfile','Phantom1D/Phantom_data.mat', 'spect_infofile', 'Phantom1D/Phantom_spectrm_info.mat','configfile','Phantom1D_ladmm.ini','spatmaskfile','Phantom1D/Phantom_mask.mat','outprefix','Phantom1D/Phantom1D_data_ladmm_spect')

Note

Here we do not show the effect of the spatial regularizer, \(\lambda\), in estimating the spectrum. You may look into the tutorial for the spectrum estimation of 2D phantom data for more details. Similar argument follows here too.

11.7.3. Output

This command will save the estimated spectrum in the file named ‘Phantom1D/Phantom1D_data_ladmm_spect.mat’.

Phantom1D_data_ladmm_spect.mat — spectroscopic image file (see here)
- spectral_image: double matrix, size [300,28,38], the 70x70 spectrum for each of the 28x38 volxels in one file.
- axes: 1×5 struct with fields: type (char), sample (1×N vector), name (char), unit (char), spacing ('log' or 'linear')
  - axis 1 (relaxation): type, 'spectral', name 'T_2', unit ‘\(ms\)’, 300 samples, spacing 'log'.
  - axis 2 : type, 'spatial', name 'x', unit ‘\(mm\)’, empty samples, spacing 'linear'.
  - axis 3 : type, 'spatial', name 'y', unit ‘\(mm\)’, empty samples, spacing 'linear'.
- spectral_dim: 300, storing the spectral dimension for axes 1 .
- spatial_dim`: ``[28, 38], storing the spatial dimension for axes 2 and 3 respectively.
- transform: 4×4 affine (voxel → world, Identity matrix)
- resolution: [Δx, Δy, Δz] in mm.

Note that some variables in this file are directly transferred from the 'Phantom_spectrm_info.mat' file. In the next tutorial we will be visualising the estimated spectra.