• 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)
• Principal Component Analysis
  • Venables, W. N. and B. D. Ripley (2002). Modern Applied Statistics with S, Springer-Verlag
• Shape Tensors
  • http://www.cs.utah.edu/~gk/papers/vissym04/
  • Doube, M.,et al. (2010). BoneJ: Free and extensible bone image analysis in ImageJ. Bone, 47, 1076–9. doi:10.1016/j.bone.2010.08.023
  • Mader, K. , et al. (2013). A quantitative framework for the 3D characterization of the osteocyte lacunar system. Bone, 57(1), 142–154. doi:10.1016/j.bone.2013.06.026

Models / ROC Curves

• Julia Evans - Recalling with Precision
• Stripe's Next Top Model

• True Positive values in the ring that are classified as Foreground
• True Negative values outside the ring that are classified as Background
• False Positive values outside the ring that are classified as Foreground
• False Negative values in the ring that are classified as Background

ROC Curve

• Recall (sensitivity)= $TP/(TP+FN)$
• Precision = $TP/(TP+FP)$

Reciever Operating Characteristic (first developed for WW2 soldiers detecting objects in battlefields using radar). The ideal is the top-right (identify everything and miss nothing)

While finding a ring might be didactic, it is not really a relevant problem and these terms are much more meaningful when applied to medical images where every False Positives and False Negative can be mean life-threatening surgery or the lack thereof. (Data courtesy of Zhentian Wang)

From these images, an expert labeled the calcifications by hand, so we have ground truth data on where they are:

We can perform the same analysis on an image like this one, again using a simple threshold to evalulate how accurately we identify the calcifications

We can get even better results by using a combination of both absorption and scattering to find the calcifications

This can then be represented as a pseudo-ROC curve by showing all of the points and finding the maximum value.

• Motivation (Why and How?)
• Object Characterization
• Volume
• Center and Extents
• Anisotropy

• Shape Tensor
• Principal Component Analysis
• Ellipsoid Representation
• Scale-free metrics
• Anisotropy, Oblateness
• Meshing
  • Marching Cubes
  • Isosurfaces
• Surface Area

We have dramatically simplified our data, but there is still too much.

• We perform an experiment bone to see how big the cells are inside the tissue $$\downarrow$$

2560 x 2560 x 2160 x 32 bit

56GB / sample

• Filtering and Enhancement!
  $$\downarrow$$
• 56GB of less noisy data

• Segmentation

$$\downarrow$$

2560 x 2560 x 2160 x 1 bit

(1.75GB / sample)

• Still an aweful lot of data

Once we have a clearly segmented image, it is often helpful to identify the sub-components of this image. The easist method for identifying these subcomponents is called component labeling which again uses the neighborhood $\mathcal{N}$ as a criterion for connectivity, resulting in pixels which are touching being part of the same object.

In general, the approach works well since usually when different regions are touching, they are related. It runs into issues when you have multiple regions which agglomerate together, for example a continuous pore network (1 object) or a cluster of touching cells

The more general formulation of the problem is for networks (roads, computers, social). Are the points start and finish connected?

We start out with the network and we grow Begin to its connections. In a brushfire-style algorithm

• For each point $(x,y)\in\{\text{Begin},\text{Pixels Grown}\}$
  • For each point $(x^{\prime},y^{\prime})\in\mathcal{N}(x,y)$
  • Set the label to Pixels Grown
• Repeat until no more labels have been changed

Same as for networks but the neighborhood is defined with a kernel (circular, full, line, etc) and labels must be generate for the image

• Assign a unique label to each point in the image $$L(x,y) = y*\text{Im}_{width}+x$$

• For each point $(x,y)$

  • Check each point $(x^{\prime},y^{\prime})\in\mathcal{N}(x,y)$
  • Set $$L(x,y)=\text{min}(L(x,y),L(x^{\prime},y^{\prime}))$$ $$L(x^{\prime},y^{\prime})=L(x,y)$$

• Repeat until no more $L$ values are changed

The image very quickly converges and after 4 iterations the task is complete. For larger more complicated images with thousands of components this task can take longer, but there exist much more efficient algorithms for labeling components which alleviate this issue.

What we would like to to do

• Count the cells
• Say something about the cells
• Compare the cells in this image to another image
• But where do we start?

$$I_{id}(x,y) = \begin{cases} 1, & L(x,y) = id \\ 0, & \text{otherwise} \end{cases}$$

Define a center

$$\bar{x} = \frac{1}{N} \sum_{\vec{v}\in I_{id}} \vec{v}\cdot\vec{i}$$ $$\bar{y} = \frac{1}{N} \sum_{\vec{v}\in I_{id}} \vec{v}\cdot\vec{j}$$ $$\bar{z} = \frac{1}{N} \sum_{\vec{v}\in I_{id}} \vec{v}\cdot\vec{k}$$

If the gray values are kept (or other meaningful ones are used), this can be seen as a weighted center of volume or center of mass (using $I_{gy}$ to distinguish it from the labels)

Define a center

$$\Sigma I_{gy} = \frac{1}{N} \sum_{\vec{v}\in I_{id}} I_{gy}(\vec{v})$$ $$\bar{x} = \frac{1}{\Sigma I_{gy}} \sum_{\vec{v}\in I_{id}} (\vec{v}\cdot\vec{i}) I_{gy}(\vec{v})$$ $$\bar{y} = \frac{1}{\Sigma I_{gy}} \sum_{\vec{v}\in I_{id}} (\vec{v}\cdot\vec{j}) I_{gy}(\vec{v})$$ $$\bar{z} = \frac{1}{\Sigma I_{gy}} \sum_{\vec{v}\in I_{id}} (\vec{v}\cdot\vec{k}) I_{gy}(\vec{v})$$

Exents or caliper lenghts are the size of the object in a given direction. Since the coordinates of our image our $x$ and $y$ the extents are calculated in these directions

Define extents as the minimum and maximum values along the projection of the shape in each direction $$\text{Ext}_x = \left\{ \forall \vec{v}\in I_{id}: max(\vec{v}\cdot\vec{i})-min(\vec{v}\cdot\vec{i}) \right\}$$ $$\text{Ext}_y = \left\{ \forall \vec{v}\in I_{id}: max(\vec{v}\cdot\vec{j})-min(\vec{v}\cdot\vec{j}) \right\}$$ $$\text{Ext}_z = \left\{ \forall \vec{}\in I_{id}: max(\vec{v}\cdot\vec{k})-min(\vec{v}\cdot\vec{k}) \right\}$$

• Lots of information about each object now
• But, I don't think a biologist has ever asked “How long is a cell in the $x$ direction? how about $y$?”

By definition (New Oxford American): varying in magnitude according to the direction of measurement.

• It allows us to define metrics in respect to one another and thereby characterize shape.
• Is it tall and skinny, short and fat, or perfectly round

Due to its very vague definition, it can be mathematically characterized in many different very much unequal ways (in all cases 0 represents a sphere)

$$Aiso1 = \frac{\text{Longest Side}}{\text{Shortest Side}} - 1$$

$$Aiso2 = \frac{\text{Longest Side}-\text{Shortest Side}}{\text{Longest Side}}$$

$$Aiso3 = \frac{\text{Longest Side}}{\text{Average Side Length}} - 1$$

$$Aiso4 = \frac{\text{Longest Side}-\text{Shortest Side}}{\text{Average Side Length}}$$

$$\cdots \rightarrow \text{ ad nauseum}$$

Let's take some sample objects

Anisotropy vs $x$ extents

Objects with uniformally distributed, independent $x$ and $y$ extents

While many of the topics covered in Linear Algebra and Statistics courses might not seem very applicable to real problems at first glance, at least a few of them come in handy for dealing distributions of pixels (they will only be briefly covered, for more detailed review look at some of the suggested material)

Principal Component Analysis

Similar to K-Means insofar as we start with a series of points in a vector space and want to condense the information. With PCA instead of searching for distinct groups, we try to find a linear combination of components which best explain the variance in the system.

As an example we will use a very simple example of corn and chicken prices vs time

The first principal component condenses the correlated information in both the chicken and corn prices (perhaps the underlying cost of fuel) since it explains the most variance in the final table of corn and chicken prices.

The second principal component is then related to the unique information seperating chicken from corn prices but neither indices directly themselves (maybe the cost of antibiotics)

How do these statistical analyses help us?

Going back to a single cell, we have the a distribution of $x$ and $y$ values.

• are not however completely independent
• greatest variance does not normally lie in either x nor y alone.

A principal component analysis of the voxel positions, will calculate two new principal components (the components themselves are the relationships between the input variables and the scores are the final values.)

• An optimal rotation of the coordinate system

While the eigenvalues and eigenvectors are in their own right useful

• Not obvious how to visually represent these tensor objects
• Ellipsoidal (Ellipse in 2D) representation alleviates this issue

Ellipsoidal Representation

1. Center of Volume is calculated normally
2. Eigenvectors represent the unit vectors for the semiaxes of the ellipsoid
3. $\sqrt{\text{Eigenvalues}}$ is proportional to the length of the semiaxis ($\mathcal{l}=\sqrt{5\lambda_i}$), derivation similar to moment of inertia tensor for ellipsoids.

• The bounding box matches boundaries better
• The ellipse captures the dimensions and shape better
• Extents can be calculated from either method
• Anisotropy can also be calculated from either method

We see that there seems to be a general, albeit weak, correlation between the two measures. The most concerning portion is however the left side where the extents or bounding box method reports 0 anisotropy and the elliptical method reports substancial amounts of it.

The models we have done are all applicable to both 2D and 3D images. The primary difference is when looking at 3D images there is an extra dimension to consider.

It can also be shown as a series of 2D slices.

Anisotropy applies to these samples as well but it only gives us information about the shortest and the longest dimensions. Is that enough?

Each of these images has the exact same anisotropy because the shortest semiaxis length remains

If we take a slice through the image we can image the final shape looking like a very small thin rod in the case where the second dimension is equal to 3 (short in two directions and long in the other) and a pankcake (short in one direction and long in the other two).

We can thus introduce a new metric for assessing the second-degree anisotropy in the object and this we shall somewhat more formally call Oblateness.

$$\textrm{Ob} = 2\frac{\lambda_{2}-\lambda_{1}}{\lambda_{3}-\lambda_{1}}-1$$

The value like in anisotropy is bound between -1 and 1.

• -1 indicates the middle semiaxis is the same length of the shortest semiaxis $\lambda_2=\lambda_1$ and thus the structure is rod-like $\rightarrow$ prolate
• +1 indicates the middle axis is the same length as the longest semiaxis $\lambda_2=\lambda_3$ and thus the structure is pancake-like $\rightarrow$ oblate

We see that the dilation and erosion affects are strongly related to the surface area of an object: the more surface area the larger the affect of a single dilation or erosion step.

Constructing a mesh for an image provides very different information than the image data itself. Most crucially this comes when looking at physical processes like deformation.

While the images are helpful for visualizing we rarely have models for quantifying how difficult it is to turn a pixel off

If the image is turned into a mesh we now have a list of vertices and edges. For these vertices and edges we can define forces. For example when looking at stress-strain relationships in mechanics using Hooke's Model $$\vec{F}=k (\vec{x}_0-\vec{x})$$ the force needed to stretch one of these edges is proportional to how far it is stretched.

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