We illustrate several approaches of using the fast marching algorithm for the segmentation of vessels in a medical image. In particular, techniques such as dimension lifting and curvature penalization are demonstrated. The segmentation of tubular structures in images is a domain of research in itself, on which the author only has moderate expertise, hence this notebook has no pretention to be exhaustive or state of the art. It should only be regarded as an advanced demo of the HFM library.
Summary of volume Fast Marching Methods, this series of notebooks.
Main summary of the Adaptive Grid Discretizations book of notebooks, including the other volumes.
This Python® notebook is intended as documentation and testing for the HamiltonFastMarching (HFM) library, which also has interfaces to the Matlab® and Mathematica® languages. More information on the HFM library in the manuscript:
Copyright Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
import sys; sys.path.insert(0,"..") # Allow import of agd from parent directory (useless if conda package installed)
#from Miscellaneous import TocTools; print(TocTools.displayTOC('Tubular','FMM'))
from agd import Eikonal
from agd import Metrics
from agd.Plotting import savefig,imread; #savefig.dirName = 'Figures/Tubular'
from agd import AutomaticDifferentiation as ad
import numpy as np; xp=np
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # Plots geodesics
from scipy import ndimage
from itertools import groupby
from copy import deepcopy
import matplotlib.image
# Commodity functions
def DeleteConsecutiveDuplicates(L):
return np.array([x[0] for x in groupby(L, lambda x:tuple(x))])
def MyArgmin(positions,indices,cost,lvl):
return np.array([(x,y,lvl[np.argmin(cost[i,j,:])]) for ((x,y),(i,j)) in zip(positions,indices)])
Uncomment the following line to use the GPU eikonal solver
#Eikonal.dictIn.default_mode = 'gpu_transfer'
Alternatively, with a recent version of cupy
#xp,plt,ndimage,Eikonal = map(ad.cupy_friendly,(xp,plt,ndimage,Eikonal))
Several options are available to import images to be further processed by a Python script. Matplotlib provides the simplest approach, but on the author's computer it fails if the imported image is simultanously non-PNG, and internet-downloaded. In that case, the PIL library provides a workaround.
image0 = imread("Notebooks_FMM/TestImages/subImHRF31_preview.png")
In this notebook, the processed image is assumed to be in grayscale, with intensity values rescaled to the interval [0,1].
image = image0[:,:,0] # This particular image is black and white, so take any channel
image = image/np.max(image) # Rescale image values to [0,1]
image = xp.asarray(image) # Transfer to GPU if applicable
There is a surpising lack of uniformity about the way to display a two-dimensional array of grayscale values on screen. Depending on the software, coordinate axes can be exchanged, reversed, ...
In this notebook, we avoid the (i.m.o terrible) YXZ axes convention, and use an adequate imshow variant to display the results.
image = image.T # Transpose coordinates. Aka indexing = 'ij' rather than 'xy'
def imshow(image,**kwargs): plt.imshow(image.T,origin='lower',**kwargs)
# Use the origin = 'lower' option for compatibility with e.g. numpy.contourf
plt.title('The tubular structures to be segmented')
imshow(image,cmap='gray');
The next cell defines a number of points of interest in the image. The HFM library will be used to compute minimal geodesics between the appropriate seeds and tips, and with respect to the appropriate metric. In practical applications, these seeds and tips may either be automatically generated, or user defined.
seeds0 = xp.array([
[327, 15], #On main vessel
[259, 14], # On small bottom vessel
]);
seeds1 = xp.array([
[350, 16] # On secondary vessel
# [327, 92] # On small right vessel
]);
tips0 = xp.array([
[291,65], #On primary vessel, easy
[203,41], # On a small vessel, attached to the primary vessel
[250,115], #On primary vessel, after crossing
[116,111], #On primary vessel, possible shortcut
[64,24], #On primary vessel, possible shortcut, far
[12,108], #On primary vessel, two possible shortcuts
[115,54] #On small bottom vessel, after crossing
]);
tips1 = xp.array([
[302,77], #On secondary vessel, easy
[244,105], #On secondary vessel, after crossing
[189,188], #On secondary vessel, possible shortcut
[64,189], #On secondary vessel, two possible shortcuts
[20,25], #On secondary vessel, three possible shortcuts
[304,176] #On small right vessel
]);
seeds = np.concatenate((seeds0,seeds1))
tips = np.concatenate((tips0,tips1))
In this notebook, we always define the image origin as the point $(0,0)$, and using the physical grid scale $h=1$. In particular, each pixel occupies a square domain of the form $[a,a+1] \times [b,b+1]$. As a result, the (multi-)indices of the seeds and tips in the data arrays are also their 'physical' coordinates. In case of more general coordinate systems, such as in the other notebooks of this series, the conversion between physical points and array indices can be performed using the member functions IndexFromPoint
and PointFromIndex
of the agd.Eikonal.dictIn
class.
seedIndices0, seedIndices1, seedIndices = seeds0.astype(int), seeds1.astype(int), seeds.astype(int)
tipIndices0, tipIndices1, tipIndices = tips0.astype(int), tips1.astype(int), tips.astype(int)
plt.figure(figsize=(12,6))
plt.scatter(seeds[:,0],seeds[:,1],color=['cyan','red','orange']) # Seeds in red
plt.scatter(tips[:,0],tips[:,1],color='pink') # Tips in blue
imshow(image,cmap='gray');
We attach one distinct color to each seed, that will be used for coloring the backtracked geodesics.
def SeedColor(pt):
seedIndex = np.argmin((seeds[:,0]-pt[0])**2+(seeds[:,1]-pt[1])**2)
return ['cyan','red','orange'][seedIndex]
# Named colors : 'blue', 'orange', 'green', 'red', 'purple', 'brown', 'pink', 'gray', 'olive', 'cyan'
The proposed approach to tubular structure segmentation in images works in two steps:
We use local filters to detect, in a preliminary image filtering, wether a vessel or the edge of an object is likely present at a point of the provided image. The chosen filters are closely related with the Gabor filter, up to an additional anisotropy parameter. More precisely, our convolution kernels depend on several parameters
The following function returns a local system of coordinates $(X,Y)$ in the plane, a Gaussian bump, and two convolution filters that are closely related with the real and imaginary part of a Gabor filter. They can be used to detect, respectively, tube-like features and object boundaries.
def FeatureKernels(radius, theta, anisotropy, darkObjects = False):
# Construct coordinate system
side = int(radius * anisotropy * 1.4)+1
ran = xp.arange(-side,side+1)
X,Y = xp.meshgrid(ran,ran,indexing='ij')
# Rotate and rescale the coordinate system
Xr,Yr = np.cos(theta)*X+np.sin(theta)*Y, -np.sin(theta)*X+np.cos(theta)*Y
Xr,Yr = Xr/radius, Yr/radius
bump = np.exp(-((Xr/anisotropy)**2+Yr**2)/2) # Anisotropic gaussian
bump = bump/bump.sum() # Normalized bump function
#Constant 2 is somewhat arbitrary
edge = bump*np.sin(2*Yr); # Similar to imaginary part of Gabor filter
vessel = bump*np.cos(2*Yr); # Similar to real part of Gabor filter
vessel = vessel - vessel.sum() * bump # Filter should have zero mean.
sign= -1 if darkObjects else 1
return [X, Y, bump, sign*vessel, sign*edge]
X, Y, meanFilter, vesselFilter, edgeFilter = FeatureKernels(3,np.pi/6,2,darkObjects=True)
The first filter, the vessel filter, is efficient at detecting (the centerline of) tubular structures, whose radius and orientation match the parameters.
fig = plt.figure(figsize=[3,3]); plt.title('Filter for vessel detection.');
plt.contourf(X,Y,vesselFilter); plt.axis('equal');
plt.title('Image response to the vessel detection filter.');
imshow(ndimage.convolve(image,vesselFilter),cmap='gray'); #plt.colorbar();
The second filter can be used to detect object boundaries. However, this is irrelevant for the present application, which is the detection of tubular centerline. For instance, two cells below, we obtain a positive and a negative reponse on the side of the vessel, whose radius and orientation match the parameters. In contrast the previous filter yields a positive response in the vessel centerline, which is of direct interest.
fig = plt.figure(figsize=[3,3]); plt.title('Filter for edge detection.');
plt.contourf(X,Y,edgeFilter); plt.axis('equal');
plt.title('Image response to the edge detection filter.');
imshow(ndimage.convolve(image,edgeFilter),cmap='gray');
In this subsection, we compute the response of the image to the Gabor-like filters for various orientation and radius parameters.
Our first step is to choose the kernel parameters. We fix the anisotropy of the kernel, and define a discretization space for the radius and orientation.
kernelAnisotropy = 2
nHalfTheta = 30
Thetas = np.linspace(0,np.pi,nHalfTheta,endpoint=False) # 30 angular orientations in [0,\pi]
radiusMin = 1.2; radiusMax = 5; nRadius = 5;
Radii = np.exp(np.linspace(np.log(radiusMin),np.log(radiusMax),nRadius))
#Radii = np.geomspace(radiusMin,radiusMax,nRadius) # Equivalent, but geomspace is not found on my computer
print(Radii,"\n",Thetas)
[1.2 1.71446426 2.44948974 3.49963551 5. ] [0. 0.10471976 0.20943951 0.31415927 0.41887902 0.52359878 0.62831853 0.73303829 0.83775804 0.9424778 1.04719755 1.15191731 1.25663706 1.36135682 1.46607657 1.57079633 1.67551608 1.78023584 1.88495559 1.98967535 2.0943951 2.19911486 2.30383461 2.40855437 2.51327412 2.61799388 2.72271363 2.82743339 2.93215314 3.0368729 ]
In the next cell, we convolve the input image with the chosen Gabor-like filters. The resulting array is $4$-dimensional, with the following coordinate system: $(r,\theta,x,y)$. A multiplicative factor factor $1/\sqrt r$ is introduced so as to slightly increase the visibility of vessels of small radius $r$.
# Note : takes a few seconds
import time; start_time = time.time()
response = np.array([[
ndimage.convolve(image, FeatureKernels(radius,theta,kernelAnisotropy,darkObjects=True)[3] )/ np.sqrt(radius)
for theta in Thetas] for radius in Radii])
print("--- %s seconds ---" % (time.time() - start_time))
--- 4.5042970180511475 seconds ---
The above computed $4$ dimensional response array is unfortunately already too high dimensional for most practical applications. Indeed, in order to keep execution times low and to user experience good, most PDE-based numerical methods for image processing involve domains of dimension $2$ or $3$ only. In the next cell, we perform dimension crushing by taking the maximum of the filter response w.r.t. the radius, the angle, or both.
responseRadius = response.max(axis=1).transpose(1,2,0) # Radius dependent response
responseTheta = response.max(axis=0).transpose(1,2,0) # Orientation dependent response
response0 = responseRadius.max(axis=2) # Position only dependent response
print('Original image size',response0.shape,
', radius lifted:',responseRadius.shape,
', orientation lifted',responseTheta.shape,'.')
Original image size (375, 200) , radius lifted: (375, 200, 5) , orientation lifted (375, 200, 30) .
The next figure shows the maximum, over all possible orientations and radii, of the image response to the Gabor-like filter. We regard this as a measure of vesselness.
plt.title("Measure of vesselness")
imshow(response0); #plt.colorbar();
The radius-dependent response does, as expected and desired, detect large vessels at large scales and small vessels at small scales.
plt.figure(figsize=(16,8))
for i in range(5):
plt.subplot(2,3,1+i)
plt.title(f"Response with radius {np.round(Radii[i],2)}")
imshow(responseRadius[:,:,i]);
The orientation-dependent response is positive only for the vessels parts whose orientation matches the one of the filter.
plt.figure(figsize=(16,8))
for i in range(6):
plt.subplot(2,3,1+i)
iTheta = int(nHalfTheta*i/6)
plt.title(f"Response with angle {np.round(Thetas[iTheta],2)}")
imshow(responseTheta[:,:,iTheta]);
#fig=plt.figure(figsize=(6,3.5)); plt.axis('off')
#imshow(responseTheta[:,:,5]); plt.title(r'Vessel detection at orientation $\theta=\pi/6$');
#imshow(responseTheta[:,:,0]); plt.title(r'Vessel detection at orientation $\theta=0$');
#savefig(fig,'VesselDetection1.png');
#fig=plt.figure(figsize=(6,3.5)); plt.axis('off')
#imshow(responseRadius[:,:,0]); plt.title(r'Vessel detection at minimal radius');
#imshow(responseRadius[:,:,3]); plt.title(r'Vessel detection at maximal radius');
#savefig(fig,'VesselRadius1.png');
In this section, we demonstrate vessel extraction using isotropic fast marching methods, following an approach introduced by Cohen and Kimmel (1997). The tubular structure centerlines are extracted as minimal paths, joining the prescribed seeds and tips, with respect to a metric that is locally proportionnal to the Euclidean metric: \begin{equation*} \mathrm{len}(\gamma) = \int_0^1 c(\gamma(t)) \| \gamma'(t)\| \, \mathrm d t, \end{equation*} where $c$ is a data driven cost function. Let us immediately acknowledge that the HFM library is only one of the many available software which can address this particularly simple isotropic path cost model, in contrast with the curvature penalized methods discussed in the next section.
We also illustrate the concept of dimension lifting, in which an additional abstract coordinate is introduced. The additional coordinate is independent from the image coordinates, but instead accounts for a specific feature of the vessel to be extracted. More precisely, we illustrate dimension lifting according to the radius of the extracted vessel (Li-Yezzi, 2007), and to its grayscale value (Chen-Yang-Cohen, 2013). Common dimension lifted models usually relay on diagonal rather than isotropic metrics, which means that motion along different coordinate axes may have distinct costs. (In contrast with general anisotropy, such axis aligned anisotropy raises no specific computational difficulty.)
Choice of origin.
There is a shift of one half pixel between the imshow display function and the internal representations of the eikonal solver. In order to address this discrepancy, and for the best overlay of the image and minimal geodesic paths, we shift the origin
by one half pixel in the following code cell.
isoIn=Eikonal.dictIn({
'model':'Isotropic2',
'dims':image.shape,
'gridScale':1,
'exportValues':1,
'origin':[-0.5,-0.5], #Better alignement with imshow
})
In the previous section, we convolved the input image with Gabor-like filters, so as to extract a measure of vesselness. As illustrated in the next cell, vessels appear as positive values in the processed image.
imshow(response0>0.01);
In the next cells, we define a cost function as the negative exponential of the previously computed measure of vesselness, following the original Cohen and Kimmel model. The cost function reads
$$
c(x) = c_{\min} + \exp(-c_{\exp} \Phi_0(x)),
$$
where $c_{\min}$ and $c_{\exp}$ are two constants, and $\Phi_0(x)=$response0(x)
is the vesselness measure maximized over all radii and orientations.
isoMin = 0.1
isoExp = 50
isoIn['cost']= isoMin+np.exp(-isoExp*response0)
imshow(isoIn['cost']); plt.colorbar();
In the next cells, we perform front propagation and geodesic backtracking from the first two seeds, towards the corresponding tips. For these endpoints and parameters, the results are rather satisfying.
isoIn['seeds']=seeds0
isoIn['tips']=tips0
isoOut = isoIn.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.013375 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
plt.title("Vessel segmentation using isotropic fast marching")
for geo in isoOut['geodesics']: plt.plot(*geo,color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
Recall that the geodesics presented above were backtracked after a front was propagated over all the domain, so as to compute a distance map to the seeds.
plt.title("Distance map for vessel segmentation")
imshow(isoOut['values']);
By design, the front propagation is much faster along the vessels, as illustrated by the level sets.
plt.title("Distance map for vessel segmentation")
plt.contourf(isoOut['values'].T,levels=12);
The correct extraction of the tubular structures obtained above is in part attributable to the favorable configuration of the seeds and tips. Turning to the second set of endpoints, we immediately see the shortcomings of the original Cohen and Kimmel model.
isoIn['seeds']=seeds1
isoIn['tips']=tips1
isoOut = isoIn.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.013307 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
Several issues can be identified:
#plt.figure(figsize=(16,8))
plt.title("Leaks and shortcuts in vessel segmentation using isotropic eikonal")
for geo in isoOut['geodesics']: plt.plot(*geo,color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
Radius lifting is a dimension lifting technique introduced by Li and Yezzi (2007). The radius of the tubular structures is extracted simultaneously with its position, using a three dimensional fast marching method. The vessel radius is often an important feature for medical diagonostic, hence having a reliable and mathematically well posed estimation method is a real plus.
In principle, dimension lifting techniques such as this one may also help segmentation in difficult configurations, such as the one displayed above, by disentangling the intersecting vessels in the higher dimensional domain. Unfortunately, no such improvement is noticeable here, see the next subsection for a more successful approach.
nR = Radii.size
radIn=Eikonal.dictIn({
'model':'Diagonal3',
'dims':(*image.shape,nR),
'gridScale':1,
'origin':[-0.5,-0.5,-0.5]
})
In the next cell, we define the diagonal cost of motion for our minimal path model, at any point $(x,y,r)$. We choose to distinguish:
In other words, we use the diagonal metric
$$
\begin{pmatrix}
c(x,r) & &\\
& c(x,r) &\\
& & c_r
\end{pmatrix}
\qquad \text{where} \quad
c(x,r) = c_{\min} + \exp(-c_{\exp} \Phi_{\text{rad}}(x,r)),
$$
and where $c_{\min}$ and $c_{\exp}$ are constants, and $\Phi_{\text{rad}}(x,r)$=responseRadius(x,r)
is the vesselness measure maximized over orientations.
radMin = 0.1
radExp = 50
physCost = radMin+np.exp(-radExp*responseRadius) # Cost of motion in the physical domain
radIn['metric'] = Metrics.Diagonal([physCost,physCost,2*np.ones(physCost.shape)])
The seeds and tips must now be defined as points of a three dimensional domain. For a given physical position $(x,y)$, we choose to lift the seed or tip to the radius $r$ which minimizes the cost function. This is, if we did well, the most likely radius of the vessel at this point.
radIn['seeds'] = MyArgmin(seeds0,seedIndices0,physCost,range(nR))
radIn['tips'] = MyArgmin(tips0,tipIndices0,physCost,range(nR))
radOut = radIn.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.121273 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
The segmentation results remain correct for the first set of seeds and tips, and incorrect for the second one (not shown).
for geo in radOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
However, some new information is extracted from the image: the radius of the backtraced vessels, as illustrated in the following figure. For that purpose, we need to convert radii from the linear scale used with the HFM library to the original exponential scale.
def PhysRadius(r): return Radii[0]*(Radii[1]/Radii[0])**r
print(Radii,'\n',[PhysRadius(r) for r in range(nR)])
[1.2 1.71446426 2.44948974 3.49963551 5. ] [1.2, 1.7144642578192795, 2.4494897427831774, 3.4996355115805815, 4.999999999999997]
Note : Correctly and reliably setting the pix2pt
conversion factor in the next code snippet, from input image pixels to displayed image points, would involve a complex expression involving dpi, inches, images pixel size, ... For simplicity, this constant was adjusted by hand.
pix2pt=1.4
fig = plt.figure(figsize=(8,4))
plt.title('Tubular structure extraction using a radius lifted model')
for geo in radOut['geodesics']:
pts=geo[:,::10] #subsampling to make the vessel visible
plt.scatter(pts[0],pts[1],
s=pix2pt**2*PhysRadius(pts[2])**2,
marker='o',color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
savefig(fig,'VesselRadius.png');
The extracted geodesics are genuine three dimensional paths.
ax = plt.axes(projection='3d')
for geo in radOut['geodesics']: ax.plot(*geo)
A second dimension lifting technique is discussed in this subsection, based on the grayscale value of the extracted vessel, following Cohen and Chen (2013). Indeed, distinct vessels often have distinct gray levels, and this feature varies rather slowly along the vessel. (In particular, veins and arteries have different colors.) Thus grayscale intensity can be used to disentangle, in an abstract three dimensional space, the vessels intersecting each other in the original image.
In order to obtain locally consistent grayscale values, we slightly blur the input image.
grayRef = ndimage.convolve(image, FeatureKernels(2,0,1)[2] )
imshow(grayRef,cmap='gray');
We build in the next cell a scale consisting of $20$ gray intensities in the range $[0.2,0.7]$. Note that we could have instead discretized whole brightness interval $[0,1]$, but limiting ourselves to the sub-range actually encountered for the vessels saves a bit of computation time.
grayMin = 0.2
grayMax = 0.7
grayN = 40
grayScale = (grayMax-grayMin)/grayN
grayLevels = Eikonal.CenteredLinspace(grayMin,grayMax,grayN)
The cost function, for motions $(\dot x, \dot y, 0)$ along the physical dimensions, is defined as follows: $$ c(x,y,g) = c_0(x,y) \left(1+ \frac{(g-G(x,y))^2}{2 \sigma^2}\right). $$ We denoted by $g$ the coordinate in the extra 'grayscale' dimension, and by $G$ the slightly blurred image used as a grayscale reference. We denoted by $c_0$ the cost function used in the experiment with the classical Cohen and Kimmel model, and by $\sigma>0$ a standard deviation in the grayscale dimension that is to be determined. By construction, this cost strongly increases when the current grayscale value differs from the image reference by more than $\sigma$.
baseMin = 0.1
baseExp = 50
baseCost = baseMin+np.exp(-baseExp*response0)
graySigma = 3*grayScale
grayCost = np.stack([baseCost*(1+(grayRef-lvl)**2 / (2*graySigma**2) ) for lvl in grayLevels],axis=-1)
Some slices of the cost function, corresponding to several coordinates in the grayscale dimension, are shown in the following cell.
plt.figure(figsize=(16,8))
for i in range(6):
plt.subplot(2,3,1+i)
igray = 2+6*i
plt.title(f"Cost function for gray level {igray}")
imshow(1/grayCost[:,:,igray]);
Motions along the grayscale dimension $(0,0,\dot g)$ are penalized w.r.t. a cost inversely proportionnal with the grayscale intensity: \begin{equation*} c(g) = \frac \alpha g, \end{equation*} where $\alpha>0$ is a suitable constant. This definition amounts to a uniform penalization in the logarithmic coordinates $g\mapsto \alpha \ln g$.
grayLevelCost = np.broadcast_to((100/grayLevels),(*image.shape,grayN))
We next proceed to define the test case.
grayIn=Eikonal.dictIn({
'model':'Diagonal3',
'exportValues':1,
# Domain dimensions
'dims':grayCost.shape,
'gridScales':[1,1,grayScale],
'origin':[0,0,grayMin],
'metric':Metrics.Diagonal([grayCost,grayCost,grayLevelCost]),
'seeds':MyArgmin(seeds1,seedIndices1,grayCost,grayLevels),
'tips': MyArgmin(tips1,tipIndices1,grayCost,grayLevels)
})
grayIn['stopWhenAllAccepted'] = grayIn['tips']
grayOut = grayIn.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 1.2242 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
As shown in the next cell, the backtracked minimal path is now correct for four tips out of six, shown as pink disks. In contrast, only a single backtracking was successful in the previous experiments, corresponding to the tip closest to the seed.
Note that the planar projections of several of the backtracked paths cross each other. This phenomenon is a signature of dimension lifting techniques.
plt.scatter(tips1[0:3,0],tips1[0:3,1],color='pink');
plt.scatter(tips1[5,0],tips1[5,1],color='pink');
for geo in grayOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
We can help the tubular structure extraction by exchanging the roles of the seed and of the third tip.
grayIn['tips'][2,:],grayIn['seeds'][0,:]=deepcopy(grayIn['seeds'][0,:]),deepcopy(grayIn['tips'][2,:])
grayOut = grayIn.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.765834 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
All extracted minimal geodesic now run along the secondary gray vessel, as desired, with the exception of the faint side vessel on the right.
for geo in grayOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
Regarding the first set of seeds and tips, the vessel extraction remains successful, except for a small error in the bottom, arguably excusable given the type of models considered.
grayIn['seeds']= MyArgmin(seeds0,seedIndices0,grayCost,grayLevels)
grayIn['tips'] = MyArgmin(tips0,tipIndices0,grayCost,grayLevels)
grayOut = grayIn.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 1.03565 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985 Tip {12,108,0.20625} yields 7 restarts, geodesicVolumeBound increased to 24.134
for geo in grayOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
This section is devoted to tubular segmentation using minimal paths with respect to curvature penalizing energies, of the form \begin{equation*} \int_0^L c(\gamma(t), \gamma'(t)) \mathcal C(\gamma''(t)) \mathrm d t, \end{equation*} where $\gamma$ is a path in the image domain, of euclidean length $L$ and parametrized at unit speed. The heuristic intent of curvature penalization is to add some "inertia" to the geodesics, so that they do not switch from one tubular structure to another at crossings, thus avoiding the problem of "shortcuts". The HFM library rephrases the second order energy defined above into a first order path length with respect to a degenerate metric, and computes a global optimum using a generalized fast marching algorithm.
The cost function involved in curvature penalized models depends on both the current position and tangent orientation of the path. We denote it by $c(x,\theta)$, abusively identifying the unit tangent vector with the corresponding angle. The next cell defines such a cost function, as a negative exponential of the input image response to the Gabor-like convolution filters: $$ c(x,\theta) = c_{\min} + \exp(-c_{\exp} \Phi(x,\theta)), $$ where $c_{\min}$ and $c_{\exp}$ are positive constants, and $\Phi$ is the vesselness map, maximized over radii.
orientedCost = 0.1+np.exp(-100*responseTheta)
The seeds and tips, originally defined as points of the image domain, must be lifted to three dimensional configuration space of positions and orientations.
We do this by selecting the orientation angle minimizing the cost function, at any given position. Alternatively, one could resort to the seeds_Unoriented
and tips_Unoriented
fields when calling the HFM library.
orientedSeeds = MyArgmin(seeds,seedIndices,orientedCost,Thetas)
orientedTips = MyArgmin(tips,tipIndices,orientedCost,Thetas)
imshow(-image,cmap='Greys');
plt.quiver(seeds[:,0],seeds[:,1],np.cos(orientedSeeds[:,2]),np.sin(orientedSeeds[:,2]),color='red')
plt.quiver(tips[:,0],tips[:,1],np.cos(orientedTips[:,2]),np.sin(orientedTips[:,2]),color='blue');
Note that the above definition of the oriented cost and lifted seeds and tips assumes an angular variable $\theta \in [0,\pi]$ defined up to a multiple of $\pi$. Some models require the full specification $\theta \in [0,2\pi]$ of an orientation and a direction, but not the first one that we consider.
The Reeds-Shepp "car" model would be a suitable abstraction for a wheelchair or a Segway®-like vehicle. In particular, such vehicles can rotate into place for a bounded cost, and can move both forward and backwards. From a mathematical standpoint, the Reeds-Shepp model is one of the most studied curvature penalized vehicle, thanks to its appealing sub-Riemannian structure. The function $\mathcal C$ used to penalize the path curvature $\kappa = \|\gamma''(t)\|$ is defined by \begin{equation*} \mathcal C(\kappa) := \sqrt{1+ \xi^2 \kappa^2}, \end{equation*} where $\xi$ is a parameter defined by the user, which has the dimension of a radius of curvature.
A variant, referred to as the Reeds-Shepp forward model, implemented in the HFM library but not further discussed here, allows to remove the reverse gear.
reedsIn = Eikonal.dictIn({
'model':'ReedsShepp2',
'projective':1, # Angular coordinates are restricted to the (periodic) interval [0,pi]
'xi':40,
'dims':(*image.shape,nHalfTheta),
'gridScale':1, # Physical grid scale (not angular)
'cost':orientedCost,
})
For the first set of seeds and tips, the Reeds-Shepp model extracts the desired vessels.
reedsIn['seeds']= MyArgmin(seeds0,seedIndices0,reedsIn['cost'],Thetas)
reedsIn['tips'] = MyArgmin(tips0, tipIndices0, reedsIn['cost'],Thetas)
reedsIn['stopWhenAllAccepted']=reedsIn['tips']
reedsOut = reedsIn.Run()
Field verbosity defaults to 1 Field eps defaults to 0.1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.731887 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
for geo in reedsOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
For the second, harder, collection of seeds and tips, the Reeds-Shepp model extracts the desired vessels when backtracking from four of the given tips, similarly to the grayscale lifted method, and much better than the isotropic or radius lifted approaches (one successful tip only). Observe that some of (the planar projections of) the backtracked geodesics cross each other, which we recall is a signature of dimension lifting techniques.
reedsIn['seeds']= MyArgmin(seeds1,seedIndices1,reedsIn['cost'],Thetas)
reedsIn['tips'] = MyArgmin(tips1,tipIndices1,reedsIn['cost'],Thetas)
reedsIn['stopWhenAllAccepted']=reedsIn['tips']
reedsOut = reedsIn.Run()
Field verbosity defaults to 1 Field eps defaults to 0.1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 1.09793 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
plt.scatter(tips1[0:3,0],tips1[0:3,1],color='pink');
plt.scatter(tips1[5,0],tips1[5,1],color='pink');
for geo in reedsOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
Similarly to the grayscale lifted case, we can help the backtracking by exchanging the role of the seed and of the third tip.
reedsIn['tips'][2,:],reedsIn['seeds'][0,:]=deepcopy(reedsIn['seeds'][0,:]),deepcopy(reedsIn['tips'][2,:])
reedsOut = reedsIn.Run()
Field verbosity defaults to 1 Field eps defaults to 0.1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 1.1748 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
fig=plt.figure(figsize=(6,3.5)); plt.axis('off'); plt.title('Tubular structure extraction using the Reeds-Shepp model')
for geo in reedsOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
savefig(fig,'ReedsSheppBacktracking.png')
Note the cusp, in other words the regression point, in the trajectory leading to the faint vessel on the right. In this case, and contrary to often, the cusp appears to be meaningful and well placed. Cusps can be regarded as a signature of sub-Riemannian models.
This model distinguishes itself in two ways from the previously mentioned Reeds-Shepp car:
where again $\xi$ is a user provided parameter. This cost has the physical interpretation of the bending energy of an elastic bar, hence the name.
The Euler-Mumford elastica model is a bit costly in terms of computation time, but it often yields the most satisfying results.
eulerIn=Eikonal.dictIn({
'model':'Elastica2',
'xi':60,
'dims':(*image.shape,2*nHalfTheta),
'gridScale':1,
})
The angular coordinate $\theta$ belongs to the full interval $[0,2 \pi]$, rather than $[0,\pi]$ as in our previous experiment. In more precise mathematical terms, the Euler-Mumford elastica model is posed on the configuration space $\mathbb R^2 \times \mathbb S^1$, in contrast with the projective Reeds-Shepp model discussed above which is defined on $\mathbb R^2 \times \mathbb P^1$. As a result, the array defining the cost function must be duplicated.
eulerIn['cost']=np.concatenate((orientedCost,orientedCost),axis=2)
For our first trial, we lift the seeds and tips to the orientation domain using the same routine as for the projective Reeds-Shepp model. As shown below, this naive attempt is a complete failure.
eulerIn['seeds']= MyArgmin(seeds0,seedIndices0,orientedCost,Thetas)
eulerIn['tips'] = MyArgmin(tips0, tipIndices0, orientedCost,Thetas)
eulerIn['stopWhenAllAccepted']=eulerIn['tips']
%%time
# Warning : takes up to a minute.
eulerOut =eulerIn.Run()
Field verbosity defaults to 1 Field eps defaults to 0.1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 10.4481 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985 CPU times: user 10 ms, sys: 29.3 ms, total: 39.4 ms Wall time: 10.8 s
plt.quiver(eulerIn['tips'][:,0],eulerIn['tips'][:,1],np.cos(eulerIn['tips'][:,2]),np.sin(eulerIn['tips'][:,2]),color='red')
for geo in eulerOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
The geodesics could not be successfully backtracked because their tips were not correctly oriented. There are two possible remedies to this situation:
The same reasoning applies to the seed points, but they are here, by chance, correctly oriented.
for i in [1,3,4]: eulerIn['tips'][i,2]+=np.pi # Solution 1 : reverse the direction of ill oriented tips.
eulerIn['tips_Unoriented'] = tips0 # Solution 2 : forget about the tips orientation
eulerIn['stopWhenAllAccepted']=eulerIn['tips']
Note that the computation time is several times is shorter. That is because the tips (correctly oriented, or unoriented) are reached earlier, hence the front propagation can be aborted sooner.
%%time
# Warning : takes up to a minute.
eulerOut = eulerIn.Run()
Field verbosity defaults to 1 Field eps defaults to 0.1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 1.29631 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985 CPU times: user 9.66 ms, sys: 27.3 ms, total: 36.9 ms Wall time: 1.69 s
# With correctly oriented tips
plt.quiver(eulerIn['tips'][:,0],eulerIn['tips'][:,1],np.cos(eulerIn['tips'][:,2]),np.sin(eulerIn['tips'][:,2]),color='red')
for geo in eulerOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
# With unoriented tips
for geo in eulerOut['geodesics_Unoriented']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
We next address the second, harder, set of seeds and tips. Sadly, and a bit surprisingly, the results are not as good as with the Reeds-Shepp model in this instance.
eulerIn.pop('tips_Unoriented', None);
eulerIn.pop('seeds_Unoriented',None);
eulerIn['seeds'] = MyArgmin(seeds1,seedIndices1,orientedCost,Thetas)
eulerIn['tips'] = MyArgmin(tips1,tipIndices1,orientedCost,Thetas)
eulerIn['tips'][4,2] += np.pi # Solution 1 : Reverse the uncorrectly oriented tip
eulerIn['stopWhenAllAccepted']=eulerIn['tips']
eulerOut = eulerIn.Run()
Field verbosity defaults to 1 Field eps defaults to 0.1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 2.38049 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
plt.quiver(eulerIn['tips'][:,0],eulerIn['tips'][:,1],np.cos(eulerIn['tips'][:,2]),np.sin(eulerIn['tips'][:,2]),color='red')
for geo in eulerOut['geodesics']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
eulerIn.pop('tips',None);
eulerIn.pop('seeds',None);
eulerIn.pop('stopWhenAllAccepted',None);
Trying, to exchange the roles of the seed and third tip. We set all of them unoriented for simplicity.
eulerIn['tips_Unoriented']=deepcopy(tips1)
eulerIn['seeds_Unoriented']=deepcopy(seeds1)
eulerIn['tips_Unoriented'][2,:],eulerIn['seeds_Unoriented'][0,:] \
=deepcopy(eulerIn['seeds_Unoriented'][0,:]),deepcopy(eulerIn['tips_Unoriented'][2,:])
eulerOut = eulerIn.Run()
Field verbosity defaults to 1 Field eps defaults to 0.1 Field order defaults to 1 Fast marching solver completed in 16.8078 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
for geo in eulerOut['geodesics_Unoriented']: plt.plot(geo[0],geo[1],color=SeedColor(geo[:,-1]))
imshow(-image,cmap='Greys');
# Use this cell to inspect the tips and reverse their orientation if necessary
#i=4
#plt.quiver(eulerIn['tips'][:,0],eulerIn['tips'][:,1],np.cos(eulerIn['tips'][:,2]),np.sin(eulerIn['tips'][:,2]),color='red')
#plt.quiver(eulerIn['tips'][i,0],eulerIn['tips'][i,1],np.cos(eulerIn['tips'][i,2]),np.sin(eulerIn['tips'][i,2]),color='blue')
#plt.imshow(-image,origin='lower',cmap='Greys');
As discussed in the introduction, many other models can be thought of, and implemented using the HFM library, for tubular structure segmentation. For instance, this notebook does not discuss at all Riemannian models, which are more flexible than isotropic or diagonal models, while being only marginally more expensive.