Functions

These are all the functions implemented in RETIMAT

File reading

read_vol(file, varargin)

Read metadata, segmentation, OCT, and slo images contained in a from Spectralis OCT (Heidelberg Engineering)

Input arguments

  • file: String containing the path to the .vol file to be read.

  • varargin: Optional parameters from the list:

    • ‘visu’: Visualize the scanning patter along with B-Scans and fundus image (slo).

    • ‘full_header’: Retrieve the original header with all the parameters (by default only a few important parameters are retrieved).

    • ‘get_coordinates’: retrieve fundus and A-Scan X, Y coordinates

    • ‘raw_pixel’: return raw pixel reflectance instead of visualization-adapted values.

Output arguments

  • header: Structure with .vol file header values.

  • segment: Segmenation data stored in the .vol file.

  • bscan: 3D single image with B-Scans.

  • fundus: 2D fundus image.

Notes

Spectralis OCT data can be exported into both E2E and vol format. We recommend using the latter as it provides easier access to the header information.

Originally writen by Radim Kolar, Brno University, Czech Republic Modified by Markus Mayer, Pattern Recognition Lab, University of Erlangen-Nuremberg

Modified by Kris Sheets, Retinal Cell Biology Lab, Neuroscience Center of Excellence, LSU Health Sciences Center, New Orleans

Modified by Andrew Lang, Image Analysis and Communications Lab, Johns Hopkins University - Modifications to increase efficiency - 5/17/2012

Current version modified by David Romero-Bascones, Biomedical Engineering Department, Mondragon Unibertsitatea, 2023, dromero@mondragon.edu

Examples

Read all the information in a .vol file

[header, segment, bscan, fundus] = read_vol('my_oct.vol')

Read only the header of the .vol file (faster)

header = read_vol('my_oct.vol')
read_e2e(file, varargin)

Read metadata, segmentation, OCT and fundus images in .e2e/.E2E files from Spectralis OCT (Heidelberg Engineering)

Input arguments

  • file: String containing the path to the .vol file to be read.

  • varargin: Optional parameters from the list:

    • ‘verbose’: Display header info during read.

    • ‘raw_pixel’: Display raw pixel intensities.

Output arguments

  • header: Structure with .vol file header values.

  • segment: Segmenation data stored in the .vol file.

  • bscan: 3D single image with B-Scans.

  • fundus: 2D fundus image.

Notes

This function was developed based on the previous reverse engineering attempts [1-3] and is not an official file reader. Therefore, some of the information retrieved can be incorrect/incomplete.

Scan focus, scale_x, scale_y and acquisition pattern are yet to be found.

Spectralis OCT data can be exported into both E2E and vol format. We recommend using the latter as it is easier to parse and it only stores a single eye and acquisition.

References

[1] uocte documentation https://bitbucket.org/uocte/uocte/wiki/Topcon%20File%20Format

[2] OCT-Converter, https://github.com/marksgraham/OCT-Converter

[3] LibE2E, https://github.com/neurodial/LibE2E

Examples

Read all the information in a .e2e/.E2E file

file = 'my_oct.e2e';
[header, segment, bscan, fundus] = read_e2e(file)

Read only the header of the .e2e/.E2E file (faster)

file = 'my_oct.e2e';
 header = read_e2e(file)
read_img(file, scan_size, varargin)

Read OCT images from Cirrus .img files

Input arguments

  • file: String with the path to the img file name to be read.

  • scan_size: Array of 3 numeric values defining the size of the

    scanned region. Sizes must be given following [x y z] convention. By default a cube of [6 6 2] mm is considered.

  • varargin: optional string flags from the list:

    • ‘get_coordinates’: retrieve fundus and A-Scan X, Y coordinates. If the scan pattern is unknown coordinates cannot be computed.

Output arguments

  • header: struct with metadata retrieved from filename and volume dimensions.

  • bscan: 3D matrix of size [n_axial x n_ascan x n_bscan] with image data.

Notes

img files do not include any spatial or metadata information. Therefore, to correctly reconstruct the image volume this function relies on two strategies:

  1. Inspection of the filename.

  2. Matching of the voxel count with this known acquisition patterns:

  • Macular Cube: [200 x 200 x 1024] : 40960000

  • Opitic Disc Cube: [512 x 128 x 1024] : 67108864

  • 5 Line Raster: [1024 x 5 x 1024]

  • Hidef: [512 x 2 x 1024] (not supported yet)

Scans differing from the previous protocols might not be read properly.

Info from previous function: hidef scan consists of 2 orthogonal bscans intersecting in the center of the volume. The first is in the y direction (across B-scans), the second is in the x direction (aligned with the center B-scan). Unfortunately, they are not in good alignment with the volumetric data.

Example

[header, bscan] = read_img('myfile.img')
read_bin(file, im_size)

Read grayscale slo (fundus) images from Cirrus .bin files

Input arguments

  • file: path of the .bin file to read.

  • im_size: size of the slo image. If not provided it is assumed to be [664 512];

Output arguments

  • fundus: fundus grayscale image.

Notes

This code was developed by reverse engineering a few .bin files and may not work for every .bin image.

Example

I = read_bin('myfile.bin')
read_fda(file, varargin)

Reads the header, segmentation and images (OCT + fundus) contained in an .fda Topcon file.

Input arguments

  • file: Path of the .fda file to read.

  • varargin: Optional flags from the list:

    • ‘verbose’: If provided, reading info is displayed.

    • ‘get_coordinates’: If provided A-scan coordinates are returned.

Output arguments

  • header Structure with .fda file header values.

  • segment Segmenation data stored in the .fda file.

  • bscan 3D single image with B-Scans.

  • fundus 2D fundus image.

Notes

This code is heavily based on [1,2], which were developed by reverse engineering fda files. Therefore, data read using this function may be incomplete/incorrect.

The reading process is as follows:

  1. Identify data chunks present in file (not all are always present)

  2. Read specific chunks with the relevant information

‘get_coordinates’ flag only works if the scanning pattern is horizontal raster (pending to implement others).

References

[1] uocte documentation https://bitbucket.org/uocte/uocte/wiki/Topcon%20File%20Format

[2] OCT-Converter: https://github.com/marksgraham/OCT-Converter

Example

[header, segment, bscan, fundus] = read_fda(file)
read_xml_iowa(file, varargin)

Read segmentation computed by OCTExplorer (IOWA reference algorithm)

Input arguments

  • file: path to xml output file ending with Surfaces_Iowa.xml

  • varargin: different string flags:

    • ‘get_coordinates’: to return A-Scan coordinates as part of the header.

    • ‘keep_names’: to use layer naming convention from OCTExplorer instead of RETIMAT.

    • ‘keep_nan’: to keep regions with noisy segmentation as numerical values instead of nan.

    • ‘verbose’: to plot additional info on read layers.

Output arguments

  • header: struct with metadata related with the segmentation: version, resolution, dimensions.

  • seg: struct with point segmentation values (in pixels).

Notes

This function has been only tested for raster acquisitions and OCTExplorer version version 3.8.0.

The axes convention used in OCTExplorer differs from the one used here:

  • x: horizontal (temporal - nasal)

  • y: vertical (inferior - superior)

  • z: axial (depth)

Layer naming convention returned by this function differs from the one used by OCTExplorer.

OCTExplorer uses different ways to determine an invalid region:

  • Segmentation = 0

  • Segmentation NA (not sure)

  • Undefined region chunk (full square regions I believe)

References

[1] The Iowa Reference Algorithms (Retinal Image Analysis Lab, Iowa Institute for Biomedical Imaging, Iowa City, IA), https://www.iibi.uiowa.edu/oct-reference

[2] Abramoff MD et al., Retinal Imaging and Image Analysis. IEEE Reviews in Biomedical Engineering, 2010. doi:10.1109/RBME.2010.2084567

[3] Kang L et al., Optimal Surface Segmentation in Volumetric Images – A Graph-Theoretic Approach. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2006. doi:10.1109/TPAMI.2006.19

[4] Garvin MK et al., Automated 3-D Intraretinal Layer Segmentation of Macular Spectral-Domain Optical Coherence Tomography Images. IEEE Trans Med. Imaging, 2009. doi:10.1109/TMI.2009.2016958

Example

[header, seg] = read_xml_iowa('my_file_Surfaces_Iowa.xml')

Visualization

generate_report(bscan, seg, fundus, layers, varargin)

Create a summary figure with fundus, en-face reflectance map, thickness maps and bscans.

Input arguments (mandatory)

  • bscan: Raw b-scans with shape [n_axial, n_ascan, n_bscan]

  • seg: Struct with segmentation data. Each field corresponds to a different boundary. Expected shape of each field [n_bscan, n_ascan]

  • fundus: Fundus image.

  • layers: List of layers for which to plot thickness maps.

Input arguments (optional)

  • n_plot_bscan: How many bscans to plot. Default: maximum between 5 and n_bscan.

  • n_col_max: Maximum number of columns in the figure. Default: 5

  • file_name: File name to save the plot. Default: [] (plot is not saved).

  • visible: ‘on’ to show the figure window and ‘off’ to hide it. Default: ‘on’

  • title: text to us as title. Default: [] (no title).

Output arguments

  • fig: Created figure handle.

Notes

This reports are useful to identify segmentation errors or poor signal quality.

Example

[~, segment, bscan, fundus] = read_vol(file)
generate_report(bscan, segment, fundus, {'TRT','RNFL'})
plot_sectors(Z, Sectors, varargin)

Visualize sectorization results

Input arguments (mandatory)

  • Z: Values of each sector.

  • Sectors: Struct defining the sectorization. Either manually defined or obtained after using sectorize_map().

Input arguments (name/value pair)

  • n_point: Number of points in each non-straight line segment. Higher number increases resolution and computation time. Default: 50.

  • alpha: Transparency of patches. Range: [0,1]. Default: 1.

  • edge_color: Color of the edges of each patch. Default: ‘k’ (black).

  • axis_off: Set it to true to hide axis. Default: true.

  • axis_equal: Set it to true to set a aspect ratio of 1. Default: true.

Notes

ETDRS sectorization is defined as a structure with fields:

  • type: ‘wedge’

  • radius: [0.5 1.5 3]

  • n_angle = 4

  • theta_0 = -pi/4

Example

Plot ETDRS sectorization

[header, seg] = read_vol(file, 'verbose', 'get_coordinates');
Thickness = compute_thickness(seg, 'TRT', header.scale_z);
[TRT, Sectors] = sectorize_map(X, Y, TRT, 'mean', 'etdrs');
plot_sectors(TRT, Sectors);
spider_plot(Data, varargin)

Visualize sectorization results

Input arguments (mandatory)

  • Data: matrix with N (observations)x M (angles)

Input arguments (name/value pair)

  • param: structure with parameters for plotting configuration.

    • center: ‘mean’ or ‘median’

    • superior: interval superior value % percentile (e.g. 75)

    • inferior: interval inferior value % percentile (e.g. 25)

    • circleRadius: radius of the radar outer circle. Change to hide outliers

    • plotDots: flag to overlay raw data points or not

Notes

This is a legacy function that should be updated at some point.

Reflectance

compute_attenuation(bscan, scale_z, method, seg)

Compute voxelwise attenuation coefficient based on [1] or [3].

Input arguments

  • bscan: matrix with B-scans of dimensions n_axial x n_ascan x n_bscan. If 2D matrix it will be interpreted as a single B-scan.

  • scale_z: axial (depth) resolution in mm.

  • method: method used to compute the attenuation coefficient.

    • ‘vermeer_2014’ (default)

    • ‘vanderschoot_2012’

  • seg: struct with segmentation data used with the ‘vanderschoot_2012’ method [3]

Output arguments

  • att: attenuation coefficient per voxel.

Notes

The precise measurement of tissue attenuation coefficient is difficult and and still a research topic. This function provides an approximate value based on assumptions on optical tissue properties [1]. You may want to read [2] to check how this method compares to others.

The method in [3] is currently only supported for RNFL and returns a total attenuation value for each a-scan instead of a 3D matrix.

For the methods to work it is important to use raw voxel intensity. For instance reading .vol files with the ‘raw_pixel’ flag.

References

[1] Vermeer, Depth-resolved model-based reconstruction of attenuation coefficients in optical coherence tomography, Biomedical Optics Express 2014, https://doi.org/10.1364/BOE.5.000322

[2] Chang, Review of methods and applications of attenuation coefficient measurements with optical coherence tomography, Journal of Biomedical Optics, 2019, https://doi.org/10.1117/1.JBO.24.9.090901

[3] van der Schoot, The Effect of Glaucoma on the Optical Attenuation Coefficient of the Retinal Nerve Fiber Layer in Spectral Domain Optical Coherence Tomography Images, IOVS, 2012, https://doi.org/10.1167/iovs.11-8436

Example

Load a .vol file and compute the attenuation coefficient

[header,~,bscan] = read_vol('my_file.vol','raw_pixel');
att = compute_attenuation(bscan, header.scale_z);
image_quality(I, metric, varargin)

Compute an image quality metric for each bscan in I.

Input arguments (mandatory)

  • I: 2D or 3D matrix with bscan data. If 3D data the 3rd dimension is assumed to be the bscan index.

  • metric: Metric used to compute image quality. Accepted options

    • ‘mTCI’: maximum tissue contrast index [3].

    • ‘snr’: signal to noise ratio.

    • ‘psnr’: peak to noise ratio [1].

    • ‘cnr’: contrast to noise ratio [2].

Input arguments (optional)

  • scanner: String defining the OCT scanner (when using mTCI). If not provided an arbitrary value will be used.

    • ‘Cirrus’

    • ‘RTVue’

    • ‘Spectralis’

    • ‘3D-OCT-1000’

  • seg: struct with segmentation used in ‘snr’,’psnr’,’cnr’.

Output arguments

  • iq: Image quality metric.

Notes

The definition of SNR, CNR varies across the literature. Check the references for the precise definition implemented here.

Here SNR, PSNR and CNR rely on an accurate segmentation of both the ILM and the BM. If the image is pure nois, highly ocluded, or the segmentation is wrong these metrics might not work.

A more accurate noise estimation procedure may require acquiring a NOISE profiling or reference image [4].

References

[1] Shirasawa, Objective Determination of Optimal Number of Spectral- Domain Optical Coherence Tomographic Images of Retina to Average, PLOS ONE, 2014. http://dx.doi.org/10.1371/journal.pone.0110550

[2] Rico-Jimenez J J, Real-time OCT image denoising using a self-fusion neural network,” Biomed. Opt. Express, 2022. https://doi.org/10.1364/BOE.451029

[3] Huang, Signal Quality Assessment of Retinal Optical Coherence Tomography Images, IOVS, 2012. Https://dx.doi.org/10.1167%2Fiovs.11-8755

[4] Sahu, Statistical modeling and Gaussianization procedure based de-speckling algorithm for retinal OCT images, Journal of Ambient Intelligence and Humanized Computing, 2018.

Example

Compute mTCI for a whole volume

[~, ~, bscan] = read_vol(file);
mTCI = image_quality(bscan, 'mTCI', 'Spectralis');
normalize_reflectance(bscan, seg, method)

Normalize the intensity (reflectance) of one or more bscans based on the reflectance of the RPE layer.

Input arguments

  • bscan: Matrix with B-scans of dimensions n_axial x n_ascan x N_bscan. If 2D matrix then it will be interpreted as a single B-scan.

  • seg: Struct with the segmentation of the ILM, IZ_RPE and BM.

  • method: Method used to normalize.

    • ‘ascan’: to normalize each A-scan separately

    • ‘bscan’: to normalize each B-scan as a whole (default).

Output arguments

  • bscan_norm: Normalized bscan images.

Notes

This function relies on an accurate segmentation of the ILM, IZ_RPE and BM.

The normalization sets the low and hight reflectance references to be the vitreous and the RPE, respectively. Can be thought as the % of deviation from those parametesr. Values outside [0,100] are expected.

The performance might be somewhat slow due to the necessity of looping through all the A-scans. There might be computational improvements using poly2mask but that needs to be adapted to handle NaNs properly (usual in segmentation). It might be also cleaner to use here reflectance_map to avoid duplicated code.

References

[1] Sharafeldeen, Precise higher‐order refectivity and morphology models for early diagnosis of diabetic retinopathy using OCT images, Scientific Reports, 2021. https://doi.org/10.1038/s41598-021-83735-7

Example

[~,seg,bscan] = read_vol('my_file.vol');
bscan_norm = normalize_reflectance(bscan, seg);
reflectance_map(bscan, method, metric, seg, varargin)

Create a 2D map (en-face image) of each A-Scan reflectance

Input arguments (mandatory)

  • bscan: 3D Volume with b-scans images.

  • method: Method to compute reflectance.

    • ‘raw’ (default)

    • ‘normalized’

    • ‘attenuation’

  • metric: Metric to be used.

    • ‘mean’

    • ‘total’

    • ‘layer_index’

  • seg: Struct with boundary segmentation data (in voxel units measured from the top of each B-Scan). The dimensions must match the provided volume (bscan). If not provided, all the voxels in each A-scan are used.

Input arguments (optional)

  • scale_z: axial (depth) resolution of the image. Necessary if the metric is ‘total reflectance’.

  • top: Name of the upper boundary of the layer to be analyzed. It must correspond to a field in seg.

  • bottom: Name of the bottom boundary of the layer to be analyzed. It must correspond to a field in seg.

If only one optional argument is provided that is assumed to be scale_z. When two optionals are provided then those are considered as top/bottom boundaries. When the three are provided they are assumed in the order: scale_z, top and bottom.

Output arguments

  • R: 2D matrix with reflectance values.

Notes

Reflectance may be derived from raw voxel intensities, normalized intensities or attenuation coefficient.

When using ‘normalized’, reflectance is corrected based on the vitreous and the RPE layer. See normalize_reflectance.m for details.

Attenuation coefficient computation relies on several assumptions on the optical properties of the tissue. See compute_attenuation.m for details.

‘layer_index’ metric is described in [1]. It includes a sort of normalization at a bscan level (the 99 percentile is used for normalization).

References

[1] Varga et al., “Investigating Tissue Optical Properties and Texture Descriptors of the Retina in Patients with Multiple Sclerosis”, PLos One, 2015

Example

[header, seg, bscan] = read_vol('my_file.vol');
R = reflectance_map(bscan, 'raw', 'total', 'ILM', 'RPE')
stack_bscans(bscan, seg, layers)

Stack the bscan portion showing the same layer into the a single image as proposed in [1].

Input arguments

  • bscan: Matrix with B-scans of dimensions n_axial x n_ascan x n_bscan.

  • seg: Struct with segmentation data.

  • layers: String or a cell array of strings with the name of the layer to be stacked. Note that, depending on the layer of

    interest, the segmentation must follow a certain naming convention (see the beggining of the function).

Output arguments

  • stacked: Struct with a stacked b-scan for each layer in layers.

Notes

The stacking process is sensitive to NaN values in the segmentation. This function performs a rough extrapolation of NaN values that may not work if the number of NaN values is high.

References

[1] Tazarjani, Retinal OCT Texture Analysis for Differentiating Healthy Controls from Multiple Sclerosis (MS) with/without Optic Neuritis, BioMed Research International, 2021. https://doi.org/10.1155/2021/5579018

Example

Generate stacked b-scans of RNFL and GCIPL layers

[~, seg, bscan] = read_vol(file);
layers = {'RNFL', 'GCIPL'};
stacked = stack_bscans(bscan, seg, layers)

Structural analysis

compute_thickness(seg, layers, scale_z)

Compute thickness the specified layers specified by using the segmentation data stored in ‘seg’ struct.

Input arguments

  • seg: struct with the segmentation of the boundaries. Each field must correspond to a specific boundary.

  • layers: string or a cell array of strings indicating the layers for which to compute thickness.

  • scale_z: (optional). Axial (depth) resolution of the images in mm. If provided, it is used to transform thicknes values from pixel to um units.

Output arguments

  • Thickness Struct with thickness values for each layer.

Notes

The naming convention of the boundaries and layers might differ both in the literature and in different segmentation algorithms. Here we mostly follow the APOSTEL 2.0 recommendations ([1]) as specified below:

  • BM: Bruch’s membrane

  • ELM: external limiting membrane

  • EZ: ellipsoid zone

  • GCL: ganglion cell layer

  • GCIPL: ganglion cell and inner plexiform layer (composite)

  • ILM: inner limiting membrane

  • INL: inner nuclear layer

  • IPL: inner plexiform layer

  • IRL: inner retinal layers (composite)

  • IZ: interdigitation zone

  • MZ: myoid zone

  • ONL: outer nuclear layer

  • ONPL: outer nuclear - plexiform layer (composite)

  • OPL: outer plexiform layer

  • OSP: outer segment of the photoreceptors

  • RNFL: retinal nerve fiber layer

  • RPE: retinal pigment epithelium

For boundaries and layer composites not specified in [1] we have used underscores to denote top/bottom boundaries. To consider the following naming case:

  • EZ_OSP: boundary between EZ and OSP

  • EZOSP: composite layer with EZ + OSP

References

[1] Aytulun et al., “APOSTEL 2.0 Recommendations for Reporting Quantitative Optical Coherence Tomography Studies”, Neurology, 2021 https://doi.org/10.1212/WNL.0000000000012125

Example

Compute thickness for the NFL and GCL layers

[header, seg, ~, ~] = read_vol(file,'verbose', 'coordinates');
Thickness = compute_thickness(seg, {'NFL','GCL'}, header.scale_z);
fit_pit_model(theta, rho, Z, pit_model, varargin)

Fit a mathematical model to the foveal pit

Input arguments (mandatory)

  • theta: Matrix with theta coordinates (polar). It expects data in format n_directions x n_points/direction.

  • rho: Matrix with rho value of polar coordinates.It expects data in format n_directions x n_points/direction.

  • Z Thickness map.

  • pit_model String defining the mathematical model to be used:

    • ‘Breher’

    • ‘Ding’

    • ‘Dubis’

    • ‘Liu’

    • ‘Scheibe’

    • ‘Yadav’

Input arguments (name-value pairs)

  • max_iter: Maximum number of iterations for each fitting. Default: 1000

  • tol_x: Tolerance of the errors during fitting. Default: 1e-6

  • tol_fun: Tolerance of the function during fitting. Default: 1e-6

  • x0’: Initial iteration coefficient values. By default specific values related to the model are used.

  • ‘x_low: Inferior limits of model coefficients.

  • x_sup: Superior limits of model coefficients.

Output arguments

  • Z_fit: Matrix with fitted values

  • Fit_coeff: Struct with estimated model coefficients

References

[1] Romero-Bascones et al., Foveal Pit Morphology Characterization: A Quantitative Analysis of the Key Methodological Steps, Entropy, 2021 https://doi.org/10.3390/e23060699

[2] Breher K. et al., Direct Modeling of Foveal Pit Morphology from Distortion-Corrected OCT Images, Biomedical Optics Express, 2019.

[3] Ding Y. et al., Application of an OCT Data-Based Mathematical Model of the Foveal Pit in Parkinson Disease, Journal fo Neural Transmission, 2014.

[4] Dubis A.M. et al., “Reconstructing Foveal Pit Morphology from Optical Coherence Tomography Imaging”, British Journal of Ophthalmology, 2009.

[5] Liu L. et al., Sloped Piecemeal Gaussian Model for Characterising Foveal Pit Shape, Ophthalmic Physiological Optics, 2016.

[6] Scheibe P. et al., Parametric Model for the 3D Reconstruction of Individual Fovea Shape from OCT Data, Experimental Eye Research, 2014.

[7] Yadav S. et al., CuBe: Parametric Modeling of 3D Foveal Shape Using Cubic Bézier, Biomedical Optics Express, 2017.

get_morph_params(rho, Z, parameters, average)

Compute morphological parameters of the foveal pit

Input arguments

  • rho: Matrix with rho coordinates (polar coordinates). Each row is a different angular direction.

  • Z: Matrix with thickness profile values. Each row is a different angular direction.

  • parameters: String or cell array of strings with the name of the parameters to computed. By default all available parameters are returned.

  • average: If true then radial parameters such as max slope are averaged across angular directions. Default is true.

Output arguments

  • X: Struct with computed features.

Notes

This parameters are usually computed for the TRT profile but can be computed as well for layers with a convex thickness profile such as GCL, IPL or INL

References

[1] Scheibe P. et al., Parametric Model for the 3D Reconstruction of Individual Fovea Shape from OCT Data, Experimental Eye Research, 2014.

[2] Yadav S. et al., CuBe: Parametric Modeling of 3D Foveal Shape Using Cubic Bézier, Biomedical Optics Express, 2017.

smooth_pit(rho, Z, span)

Smooth the foveal pit based on LOESS smoothing

Input arguments

  • rho: Matrix with radial coordinates of thickness values.

  • Z: Matrix with foveal pit values.

  • span: Percentage of samples included in each window.

    Output arguments:

    ‘Z_smooth’ Smoothed version of the thickness map.

Notes

The function expects the data in radial format: [n_angle x n_points]. To use it for entire B-Scans convert data to radial format first.

References

[1] Romero-Bascones et al., Foveal Pit Morphology Characterization: A Quantitative Analysis of the Key Methodological Steps, Entropy, 2021 https://doi.org/10.3390/e23060699

Example

[header, seg]  = read_vol('my_file.vol','coordinates');
Thickness = compute_thickness(seg, 'TRT');
[X, Y, TRT] = resample_map(header.X_oct, header.Y_oct, Thickness.TRT, ...
'star', 'n_point', 100, 'max_d', 2.4, 'n_angle', n_angle);     
[theta, rho] = cart2pol(X, Y);
TRT_smooth = smooth_pit(rho, TRT, 20);

Spatial analysis

find_fovea(X, Y, TRT, method, max_d)

Find foveal center based on total retinal thickness (TRT) map.

Input arguments

  • X: Matrix with X coordinates of each A-Scan.

  • Y: Matrix with Y coordinates of each A-Scan.

  • TRT: Matrix with total retinal thickness values.

  • method: Method used to find the foveal center. Options [‘flood’ (default), ‘min’, ‘resample_min’, ‘smooth_min’]

  • max_d: Maximum alignment error. Default: 0.85

Output arguments

  • x_fovea: X coordinate of the foveal center.

  • y_fovea: Y coordinate of the foveal center.

Notes

The smooth+min method is based on the foveafinder.m function of AURA tools. If you use it, please provide appropriate credit to the original work (https://www.nitrc.org/projects/aura_tools/).

References

[1] Romero-Bascones et al., Foveal Pit Morphology Characterization: A Quantitative Analysis of the Key Methodological Steps, Entropy, 2021 https://doi.org/10.3390/e23060699

Example

file = '../data/raster.vol';
[header, seg, ~, ~] = read_vol(file, 'coordinates');
Thickness = compute_thickness(seg, 'TRT', header.scale_z);
[x_fovea, y_fovea] = find_fovea(header.X, header.Y, Thickness.TRT)
flip_coordinates(X, eye, eye_ref)

Flip X coordinates if eye is equal to eye_ref

Input arguments

  • X: input x coordinates (in any shape)

  • eye: string with eye to be tested

  • eye_ref: reference to decide when to flip (usually OS)

Output arguments

  • X: flipped coordinates

get_ascan_coordinates(header)

Compute A-scan coordinates (X, Y) based on the scanning protocol

Input arguments

  • header: File header obtained when reading an OCT file.

Output arguments

  • X: 2D matrix with X coordinates pointing temporal to nasal.

  • Y: 2D matrix with Y coordinates pointing inferior to superior.

resample_map(X, Y, Z, grid_type, varargin)

Resample a 2D map into a new grid by interpolation/extrapolation

Input arguments (mandatory)

  • X: Original X coordinates.

  • Y: Original Y coordinates.

  • Z: Original Z values (e.g., thickness values).

  • grid_type: String specifying the resampling grid type. Depending on the grid type it must be followed by specific extra arguments. Options: [‘regular’, ‘star’]

Input arguments (name/value pair)

  • max_d: Maximum X and Y values (‘regular’) or radius (‘star’)

  • n_point: If grid_type == ‘regular’, the number of points in X and Y directions. If ‘star’, the number points from 0 to max_d

  • n_angle: Number of directions to use (when the grid_type is ‘star’)

  • theta_0: Initial angle for star resampling

  • interp_method: Interpolation method [‘linear’ (default),’cubic’]

  • extrapolate If true then points with NaN outside the original range are extrapolated. Default: false

Output arguments

  • X1 New X coordinates.

  • Y1 New Y coordinates.

  • Z1 New Z values.

Notes

If the region to cover is bigger than the original data region it might not be possible to extrapolate values accurately.

Example

Resample original A-scans into a regular grid

[header, seg, ~, ~] = read_vol(myfile.vol, 'coordinates');
Thickness = compute_thickness(seg, 'TRT', header.scale_z);

[X, Y, TRT] = resample_map(header.X_oct, header.Y_oct, Thickness.TRT, ...
'regular', 'n_point', 100, 'max_d', 2.5);
sectorize_map(X, Y, Z, metric, sector_type, varargin)

Sectorize a 2D map into several sectors (e.g., ETDRS average)

Input arguments

  • X: X coordinates of map points

  • Y: Y coordinates of map points

  • Z: Z coordinates of map points

  • metric metric to be used to sectorize data. Options: [‘mean’, ‘std’, skewness’, ‘kurtosis’]

  • sector_type: sectorization definition. Two options:

    • String defining the sectorization type. It must be followed by appropriate arguments (see examples).

    • Struct with sectorization info created beforehand.

  • varargin: extra arguments to define the sectorization.

Output arguments

  • Zs: Sectorized values for each sector.

  • Sectors: Struct with spatial data defining the sectorization. Useful when the Sectors input variable is a string.

Notes

All angles are expected in radians. To convert from degrees to radians use deg2rad() function.

Examples

ETDRS sectorization

[header, seg, ~, ~] = read_vol(file,'verbose', 'get_coordinates');
Thickness = compute_thickness(seg, 'TRT', header.scale_z);
[X, Y, TRT] = resample_map(header.X_oct, header.Y_oct, Thickness.TRT, ...
     'regular', 'n_point', 100, 'max_d', 2.5);
[Z, G] = sectorize_map(X, Y, TRT, 'mean', 'etdrs');

2 ring sectorization

[header, seg, ~, ~] = read_vol(file,'verbose', 'get_coordinates');
Thickness = compute_thickness(seg, 'TRT', header.scale_z);
[X, Y, TRT] = resample_map(header.X_oct, header.Y_oct, Thickness.TRT, ...
     'regular', 'n_point', 100, 'max_d', 2.5);
[Z, G] = sectorize_map(X, Y, TRT, 'mean', 'ring', [0.5 1.5 3]);
sectorize_peripapillar(X, Y, Z, sect_type, varargin)

Sectorize thickness values of a peripapillar circular scan

Input arguments

  • X: X coordinates of map points

  • Y: Y coordinates of map points

  • Z: Z coordinates of map points

  • sector_type: String defining the sectorization type:

    • ‘average’

    • ‘4_quadrant’

    • ‘qustom’: must be followed by the number of angles and the initial angle (see examples).

  • varargin: Extra arguments to define the sectorization.

Output arguments

  • Zs: Sectorized values for each sector.

Notes

X and Y coordinates must follow the temporal-nasal and inferior-superior conventions.

‘4_quadrant’ sectorization returns values in the order: ‘nasal’, ‘superior’, ‘temporal’ and ‘inferior’.

Examples

4 quadrants sectorization

[header, seg, ~, ~] = read_vol(file,'verbose', 'get_coordinates');
Thickness = compute_thickness(seg, 'TRT', header.scale_z);
Z = sectorize_peripapillar(X, Y, TRT, '4_quadrant');

Qustom number of quadrants sectorization

[header, seg, ~, ~] = read_vol(file,'verbose', 'get_coordinates');
Thickness = compute_thickness(seg, 'TRT', header.scale_z);
Z = sectorize_peripapillar(X, Y, TRT, 'qustom', 8, -pi/8);

Texture analysis

fractal_dimension(I, method, visu)

Compute Fractal Dimension of a grayscale image

Input arguments

  • I: Input image as a 2D matrix. Color images are converted to gray scale values.

  • method: String indicating the method to use for fractal dimension computation.

    • ‘DBC’ [1]

    • ‘IR_DBC’ (default) [2]

  • visu: If true a plot with the computed values is shown. Default = false

Output arguments

  • FD: Fractal dimension.

Notes

This implementation does not incorporate a procedure to compute the fractal dimension of neither color nor binary images. The ‘DBC’ method has several assumptions on image size so ‘IR_DBC’ method is prefered.

References

[1] Sarkar N. and Chaudhuri B. “An Efficient Differential Box-Counting Approach to Compute Fractal Dimension of Image”, IEEE Transactions On Systems, Man and Cybernetics, 1991

[2] Long M. and Peng F. “A Box-Counting Method with Adaptable Box Height for Measuring the Fractal Feature of Images”, Radioengineering, 2013

Example

I = imread('cameraman.tif');
FD = fractal_dimension(I);
GLCM_features(GLCM, varargin)

Compute multiple features from a GLCM matrix

Input arguments (mandatory)

  • GLCM: 2D GLCM matrix obtained by graycomatrix()

Input arguments (optional)

  • features: Cell array of strings with the features to be computed. By default the full set of features are computed.

Output arguments

  • X: Structure with computed features.

Notes

Some features assume GLCM matrix to be symmetric

The naming convention of some parameters is messy in the literature. The following issues are important:

  • Homogeneity: equal to inverse difference and inverse difference moment.

  • Dissimilarity: equal to the difference average

  • Cluster tendency is equal to sum variance

  • Intertia is equal to contrast

Additional features not yet implemented: maximal correlation coefficient

References

[1] Haralick R. Shanmugam K. and Dinstein I. “Textural Features for Image Classification”, IEEE Transactions on Systems, Man and Cybernetics, 1973 http://haralick.org/journals/TexturalFeaturesHaralickShanmugamDinstein.pdf#page10

[2] Soh L. “Texture Analysis of SAR Sea Ice Imagery Using Gray Level Co-Ocurrence Matrices”, IEEE Transations on Gegoscience and Remote Sensing, 1999. https://doi.org/10.1109/36.752194

[3] Zwanenburg et al. The Image Biomarker Standardization Initiative: Standardized Quantitative Radiomics for High-Throughput Image-based Phenotyping, 2020, Radiology, https://pubs.rsna.org/doi/full/10.1148/radiol.2020191145 Implementation: https://pyradiomics.readthedocs.io/en/latest/features.html?highlight=imc#radiomics.glcm.RadiomicsGLCM.getImc1FeatureValue

Example

I = imread('cameraman.tif');
M = graycomatrix(I,'NumLevels',9,'G',[])
X = GLCM_features(M)
lacunarity(I, method, r, visu)

Compute lacunarity of a gray-scale/binary image

Lacunarity is a measure of how an image covers the space. Images/masks with bigger gaps will show a higher lacunarity values. Instead of a single value lacunarity is often computed for a set of scales defined by the size of the box used to compute it. It is usually considered as a complementary metric to fractal dimension.

Input arguments

  • I: Input grayscale/binary image.

  • method: Method to compute lacunarity.

    • ‘window’: (default) fast and simple. See [2] Section 2.2.

    • ‘local_bin’: See [3] Section 3.1.

    • ‘box_3d’: Slow. See [3] Section 3.2.

  • r: 1D array with window size values. By default 10 values from 1 to the minimum between rows and columns (min{M, N}) are used.

  • visu: If true a plot with the computed values is shown. Default = false

Output arguments

  • L: 1D array with lacunarity values.

  • r 1D array with box size values.

Notes

Color images will be converted to grayscale. There are multiple definitions of lacunarity, with different calculation methods. This function covers 3 methods with different approaches. The scale of each method is not directly comparable.

There is also a “normalized Lacunarity” proposed in [2] that may be interesting to implement in the future.

References

[1] Allain C. and Clitre M. “Charaterizing the lacunarity of random and deterministic fractal sets”, Physical Review, 1991. https://doi.org/10.1103/PhysRevA.44.3552

[2] Roy A. and Perfect E. “Lacunarity Analyses of Multifractal and Natural Grayscale Patterns”, Fractals, 2014 https://doi.org/10.1142/S0218348X14400039

[3] Backes A.R. “A new approach to estimate lacunarity of texture images”, Pattern Recognition Letters, 2013

Example

Compute lacunarity

I = imread('cameraman.tif');
[L, s] = lacunarity(I)

Compute Lacunarity choosing a method and custom r values

I = imread('cameraman.tif');
r = 2:10;
[L, ~] = lacunarity(I, 'local_bin', r);
LBP_features(I, n_neighbor, features)

Compute features using Local Binary Pattern

Input arguments

  • I: Input grayscale image to be analyzed.

  • n_neighbor: Number of neighbors. Default is 8.

  • features: String or a cell array of strings defining the nmerical features to be computed. By default all the features are returned. If ‘none’ then no features are computed and only the LBP histogram is returned.

Output arguments

  • X: Struct with features computed from LBP histogram.

  • H: LBP histogram.

Example

I = imread('cameraman.tif');
X = LBP_features(I, 8)