SE2_RetinaTracking (the Wolfram Language for Students - Personal Use Only : www.wolfram.com)

---Preamble

Path to the fast marching code binaries

Load the LieAnalysis package

---Functions

HFM FileIO (The old approach using disk writing)

Data IO (hard-disk <> RAM)

HFM WLL (The new approach using a direct Wolfram Library Link)

The library link functions

Back tracking via LieAnalysis`

---Retinal Vessel Tracking

0. Detect images/test data

1. Import image

Load the stored seed points

The format of the seed point data is

with is the tip belonging to the seed.

In the stored format the coordinates are specified as follows: seed = (k,i,j), with i and j respectively the row and column index of the image ( (1,1) is the upper left pixel in the image ), and k is the discrete index for orientation. The seed points were saved at an angular resolution of sθ=2π/64, and index k=1 corresponds to θ=0. The corresponding angles for each k are given by θ=(k-1)*sθ.

Downsample data

We down sample the data by a factor sx for faster computations.
Then the discrete angle-indices are converted to continuous angles.

Note: usually we use the ordering (x,y,θ), so we would like to have (i,j,k) ordering when working with the discrete data, so we re order the coordinates as well.

2. Compute Orientation Score

Here we compute the orientation score transform of input image im. The angular resolution is specified by sθ=2π/Nθ:

Note: In the orientation scores objects (ObjPositionOrientationData), the data array is stored as a 2D multichannel image (hence the Image[ data ] in the Affix statement. The data of the orientation score can be accessed via os[“Data”], this returns a 3D array (not an image object..).

3. Vessel enhancement via multi-scale LIDOS

Here we apply vessel enhancement by taking the second order derivatives in the perpendicular vessel direction . We do this at multiple Gaussian scales, and sum the responses.

Note: the factor is a scale normalization factor, which ensures equal responses for blobs/lines at different scales. Here we included an additional parameter γ which for γ=1 does not give any bias, but for γ=2 gives a slight bias for larger scales in the reconstruction (summation).
Note: second order derivative of orientation score U is computed with the function LeftInvariantDerivative and derivative specification 22.

4. From vesselness not cost/speed function (to be used in the metric/dual metric)

The vesselness function above was only computed from 0 to π, here we construct the full vesselness function V(x,y,θ) stored from 0 to 2π.

The cost function is then defined as

with λ a parameter that balances the influence of the cost relative to no data-adaptivity, p is a contrast parameter.

In the end, in the code we only need to specify a “speed function” (dual cost), which is given by speed = 1/cost.

5. Vessel tracking (the experimental set up)

Select a seed point out of the Length[gList] number of seeds.

Here we specify tips normal arrival directions (g1dListForward), and tips with reversed arrival directions (g1dListBackward)).

Show setup

Note: in order to plot Graphics as an overlay on an image in Mathematica you have to take the Graphics coordinate system into account. Mathematica uses the convention:
(x,y), with x the horizontal axis, running from left to right, and y along the vertical axis, running from top to bottom.
See the mapping {#[[2]]-0.5,Nx-#[[1]]+0.5}&/@ in the code below to go from array coordinates (i,j) to graphics coordinates.