Books

• Jean Claude, Morphometry with R
• Online through ETHZ
• Buy it
• John C. Russ, “The Image Processing Handbook”,(Boca Raton, CRC Press)
• Available online within domain ethz.ch (or proxy.ethz.ch / public VPN)

Papers / Sites

• Thickness
  • [1] Hildebrand, T., & Ruegsegger, P. (1997). A new method for the model-independent assessment of thickness in three-dimensional images. Journal of Microscopy, 185(1), 67–75. doi:10.1046/j.1365-2818.1997.1340694.x
• Curvature
  • http://mathworld.wolfram.com/MeanCurvature.html
  • [2] “Computation of Surface Curvature from Range Images Using Geometrically Intrinsic Weights”*, T. Kurita and P. Boulanger, 1992.

• Motivation (Why and How?)
• What are Distance Maps?
• Skeletons
  • Tortuosity
• What are thickness maps?
  • Thickness with Skeletons

• Watershed Segmentation
  • Connected Objects
• Curvature
  • Characteristic Shapes

To simplify our data, but an ellipse model is too simple for many shapes

So while bounding box and ellipse-based models are useful for many object and cells, they do a very poor job with the sample below.

Why

• We assume an entity consists of connected pixels (wrong)
• We assume the objects are well modeled by an ellipse (also wrong)

What to do?

• Is it 3 connected objects which should all be analzed seperately?
• If we could divide it, we could then analyze each spart as an ellipse
• Is it one network of objects and we want to know about the constrictions?
• Is it a cell or organelle with docking sites for cell?
• Neither extents nor anisotropy are very meaningful, we need a more specific metric which can characterize

A map (or image) of distances. Each point in the map is the distance that point is from a given feature of interest (surface of an object, ROI, center of object, etc)

Definition

If we start with an image as a collection of points divided into two categories

• \( Im(x,y)= \) {Foreground, Background}
• We can define a distance map operator (\( dist \)) that transforms the image into a distance map

\[ dist(\vec{x}) = \textrm{min}(||\vec{x}-\vec{y}|| \forall \vec{y} \in \textrm{Background}) \]

We will use Euclidean distance \( ||\vec{x}-\vec{y}|| \) for this class but there are other metrics which make sense when dealing with other types of data like Manhattan/City-block or weighted metrics.

Using this rule a distance map can be made for the euclidean metric

Similarly the Manhattan or city block distance metric can be used where the distance is defined as \[ \sum_{i} |\vec{x}-\vec{y}|_i \]

The distance map is one of the crictical points where the resolution of the imaging system is important.

• We measure distances computationally in pixels or voxels
• but for them to have a meaning physically they must be converted
  • Isotropic imaging (1 \( \mu \) m x 1\( \mu \) m x 1 \( \mu \) m) is fine

Anisotropic

Ideally

• as part of filtering, resample and convert to an isotropic scale.

If that is not possible

• custom distance map algorithms which use the side-lengths of the voxels to calculate distance rather than assuming 1x1x1

We can create 2 distance maps

1. Foreground \( \rightarrow \) Background

   • Information about the objects size and interior

2. Background \( \rightarrow \) Foreground

   • Information about the distance / space between objects

One of the most useful components of the distance map is that it is relatively insensitive to small changes in connectivity.

• Component Labeling would find radically different results for these two images
  • One has 4 small circles
  • One has 1 big blob

Foreground

Background

For some structures like cellular materials and trabecular bone, we want a more detailed analysis than just thickness. We want to know

• which structures are connected
• how they are connected
• express the network in a simple manner
  • quantify tortuosity
  • branching

The first step is to take the distance transform the structure \[ I_d(x,y) = \textrm{dist}(I(x,y)) \] We can see in this image there are already local maxima that form a sort of backbone which closely maps to what we are interested in.

By using the laplacian filter as an approximate for the derivative operator which finds the values which high local gradients.

\[ \nabla I_{d}(x,y) = (\frac{\delta^2}{\delta x^2}+\frac{\delta^2}{\delta y^2})I_d \approx \underbrace{\begin{bmatrix} -1 & -1 & -1 \\ -1 & 8 & -1 \\ -1 & -1 & -1 \end{bmatrix}}_{\textrm{Laplacian Kernel}} \otimes I_d(x,y) \]

We can locate the local maxima of the structure by setting a minimum surface distance \[ I_d(x,y)>MIN_{DIST} \] and combining it with a minimum slope value \[ \nabla I_d(x,y) > MIN_{SLOPE} \]

Thresholds

Harking back to our first lectures, this can be seen as a 2D threshold of the entire dataset.

• We first make the dataset into a tuple

\[ \textrm{cImg}(x,y) = \langle \underbrace{I_d(x,y)}_1, \underbrace{\nabla I_d(x,y)}_2 \rangle \]

\[ \textrm{skelImage}(x,y) = \] \[ \begin{cases} 1, & \textrm{cImg}_1(x,y)\geq MIN-DIST \\ & \& \textrm{ cImg}_2(x,y)\geq MIN-SLOPE \\ 0, & \textrm{otherwise} \end{cases} \]

The structure is a overgrown

1. Stricter 'thresholds'
   • Sensitive to 'threshold' selection
   • Resolution dependent
2. Dedicated pruning algorithms
   • Ideally model-based
   • Minimum branch length
   • Minimum branch width

With the skeleton which is ideally one voxel thick, we can characterize the junctions in the system by looking at the neighborhood of each voxel.

As with nearly every operation, the neighborhood we define is important and we can see that if we use a box neighborhood vs a cross neighborhood (4 vs 8 adjacent) we count diagonal connections differently

Once we have classified the skeleton into terminals, segments, and junctions we can analyze the various components and assess metrics like connectivity, branching, and many other.

The easist way is to further process the image by

• component labeling the different segments and joints
• determining the average length of each segment
• number of incoming segments at a given joint

Most of the other metrics can be simply counted

Cross Neighborhood

Box Neighborhood

Thickness is a metric for assessing the size and structure of objects in a very generic manner. For every point \( \vec{x} \) in the image you find the largest sphere which:

• contains that point
• is entirely contained within the object

Taken from [1]

• The image shows a typical object
• The sphere centered at point \( p \) with a radius \( r \) is the largest that fits

Application

• Ideal for spherical and cylindrical objects since it shows their radius
• Also relevant for flow and diffusion since it can highlight bottlenecks in structure

Calculating the distance map by drawing a sphere at each point is very time consuming (\( O(n^3) \)).

• The skeleton (last section) is very closely related to the thickness.
• We found the local maxima in the image using the Laplace
• We can thus grow the Spheres from these points instead of all
• Start by instead of thresholding transforming the image to the distance at each point

\[ \textrm{thSkelImage}(x,y) = \] \[ \begin{cases} \textrm{cImg}_1(x,y) , & \textrm{cImg}_1(x,y)\geq MIN-DIST \\ & \& \textrm{ cImg}_2(x,y)\geq MIN-SLOPE \\ 0, & \textrm{otherwise} \end{cases} \]

   user  system elapsed 
  4.366   0.081   4.464 

   user  system elapsed 
 13.735   0.413  14.224 

It depends

• Small structures are lost
• They might not have been very important or noisy anyways
• Higher values are very similar

   user  system elapsed 
  1.813   0.027   1.841 

While the examples and demonstrations so far have been shown in 2D, the same exact technique can be applied to 3D data as well. For example for this liquid foam structure

• The thickness can be calculated of the background (air) voxels in the same manner.

• With care, this can be used as a proxy for bubble size distribution in systems where all the bubbles are connected to it is difficult to identify single ones.

Watershed is a method for segmenting objects without using component labeling.

• It utilizes the shape of structures to find objects
• From the distance map we can make out substructures with our eyes
• But how to we find them?!

So we start by distributing pixels all across the image

If any of the neighbors have a higher distance (more downhill) then move to that position

We can then take the points from these basins and regrow them back to their original size and these represent the watershed regions of our image

With the notion of surface normal defined (\( \vec{N} \)), we can define many curvatures at point \( \vec{x} \) on the surface.

• This is because there are infinitely many planes which contain both point \( \vec{x} \) and \( \vec{N} \)
• More generally we can define an angle \( \theta \) about which a single plane containing both can be freely rotated
• We can then define two principal curvatures by taking the maximum and minimum of this curve

\[ \kappa_1 = \textrm{max}(\kappa(\theta)) \] \[ \kappa_2 = \textrm{min}(\kappa(\theta)) \]

Mean Curvature

The mean of the two principal curvatures \[ H = \frac{1}{2}(\kappa_1+\kappa_2) \]

Gaussian Curvature

The mean of the two principal curvatures \[ K = \kappa_1\kappa_2 \]

• positive for spheres (or spherical inclusions)
  • curvatures agree in sign
• negative for saddles (hyperboloid surfaces)
  • curvatures disgree in sign
• 0 for planes

Examining a complex structure with no meaningful ellipsoidal or watershed model. The images themselves show the type of substructures and shapes which are present in the sample.

• gaussian curvature: the comparative amount of surface at, above, and below 0
  • from spherical particles into annealed mesh of balls

There are hundreds of other techniques which can be applied to these complicated structures, but they go beyond the scope of this course. Many of them are model-based which means they work well but only for particular types of samples or images. Of the more general techniques several which are easily testable inside of FIJI are

• Directional Analysis = Looking at the orientation of different components using Fourier analysis (Analyze \( \rightarrow \) Directionality)
• Tubeness / Surfaceness (Plugins \( \rightarrow \) Analyze \( \rightarrow \)) characterize binary images and the shape at each point similar to curvature but with a different underlying model

• Fractal Dimensionality = A metric for assessing the structure as you scale up and down by examining various spatial relationships
  • Ma, D., Stoica, A. D., & Wang, X.-L. (2009). Power-law scaling and fractal nature of medium-range order in metallic glasses. Nature Materials, 8(1), 30–4. doi:10.1038/nmat2340
• Two (or more) point correlation functions = Used in theoretical material science and physics to describe random materials and can be used to characterize distances, orientations, and organization in complex samples
  • Jiao, Y., Stillinger, F., & Torquato, S. (2007). Modeling heterogeneous materials via two-point correlation functions: Basic principles. Physical Review E, 76(3). doi:10.1103/PhysRevE.76.031110
  • Andrey, P., Kiêu, K., Kress, C., Lehmann, G., Tirichine, L., Liu, Z., … Debey, P. (2010). Statistical analysis of 3D images detects regular spatial distributions of centromeres and chromocenters in animal and plant nuclei. PLoS Computational Biology, 6(7), e1000853. doi:10.1371/journal.pcbi.1000853
  • Haghpanahi, M., & Miramini, S. (2008). Extraction of morphological parameters of tissue engineering scaffolds using two-point correlation function, 463–466. Retrieved from http://portal.acm.org/citation.cfm?id=1713360.1713456