Source code for ClearMap.ImageProcessing.Differentiation.Hessian

"""
Hessian
=======

Module to compute curvature measures based on Hessian Matrix

Usefull for filtering vasculature data
"""
__author__    = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__   = 'GPLv3 - GNU General Pulic License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__   = 'http://idisco.info'
__download__  = 'http://www.github.com/ChristophKirst/ClearMap2'

import numpy as np
import scipy.ndimage as ndi    

import ClearMap.IO.IO as io

import pyximport
pyximport.install(setup_args={"include_dirs":np.get_include()}, reload_support=True)

from . import HessianCode as code

#import ClearMap.IO.IO as io

__all__ = ['hessian', 'hessian_eigenvalues', 'hessian_eigensystem', 'lambda123', 'tubeness'];


###############################################################################
### Curvature measures
###############################################################################
          
[docs] def hessian(source, sink = None, sigma = None): """Returns the hessian matrix at each location calculatd via finite differences. Arguments --------- source : array Input array. sink : array Output, if None, a new array is allocated. Returns ------- hessian : array: 5d array with the hessian matrix in the first two dimensions. """ return _apply_code(code.hessian, source, sink, sink_shape_per_pixel = (3,3), parameter = None, sigma = sigma);
[docs] def hessian_eigenvalues(source, sink = None, sigma = None): """Hessian eigenvalues of source data Arguments --------- source : array Input array. sink : array Output, if None, a new array is allocated. sigma : float or None If not None, a Gaussian filter with std sigma is applied initialliy. Returns ------- sink : array The three eigenvalues of the Hessian matrix of the source. """ return _apply_code(code.eigenvalues, source, sink, sink_shape_per_pixel = (3,), parameter = None, sigma = sigma);
[docs] def tubeness(source, sink = None, threshold = None, sigma = None): """Tubeness mesure of source data Arguments --------- srouce : array Input array. sink : array Output, if None, a new array is allocated. threshold : float or None If float, the tubeness is thresholded at this level. sigma : float or None If not None, a Gaussian filter with std sigma is applied initialliy. Returns ------- sink : 3-D array Tubness output. Note ---- The tubness is the geometric mean of the two smallest eigenvalues. """ if threshold is None: return _apply_code(code.tubeness, source, sink, sigma = sigma, sink_dtype = float) else: return _apply_code(code.tubeness_threshold, source, sink, parameter = threshold, sigma = sigma, sink_dtype = bool)
[docs] def lambda123(source, sink = None, gamma12 = 1.0, gamma23 = 1.0, alpha = 0.25, sigma = None, threshold = None): """Generalized tubness measure of source data. Arguments --------- source : array Input array. sink : array Output, if None, a new array is allocated. gamma12, gamma23, alpha : float Parameters for the tubeness measure. sigma : float or None If not None, a Gaussian filter with std sigma is applied initialliy. Returns ------- sink : array The tubness measure. Note ---- Reference: Sato et al. Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images, Medical Image Analysis 1998, pp 143--168. """ if threshold is None: parameter = np.array([gamma12, gamma23, alpha], dtype = float); return _apply_code(code.lambda123, source, sink, parameter=parameter, sigma=sigma, sink_dtype=float); else: parameter = np.array([gamma12, gamma23, alpha, threshold], dtype = float); return _apply_code(code.lambda123_threshold, source, sink, parameter=parameter, sigma=sigma, sink_dtype=bool);
[docs] def hessian_eigensystem(source, sink = None, sigma = None, eigenvectors = 3): """Hessian eigensystem of source data Arguments --------- source : array Input array. sink : array Output, if None, a new array is allocated. sigma : float or None If not None, a Gaussian filter with std sigma is applied initialliy. Returns ------- sink : array The three eigenvalues and up to 3 eigenvectors of the Hessian matrix of the source. """ parameter = np.array([eigenvectors], dtype=float); sink_shape_per_pixel = (3 + eigenvectors * 3,) return _apply_code(code.eigensystem, source, sink, sink_shape_per_pixel=sink_shape_per_pixel, parameter=parameter, sigma = sigma).reshape(source.shape + (-1, 3));
#Note: for testing only def hessian_test(source, sink = None, sigma = None): """Returns the hessian matrix at each location calculatd via finite differences. Arguments --------- source : array Input array. sink : array Output, if None, a new array is allocated. Returns ------- hessian : array: 5d array with the hessian matrix in the first two dimensions. """ return _apply_code(code.eigensystem_test, source, sink, sink_shape_per_pixel = (3 + 6,), parameter = None, sigma = sigma); # def eigenvalues_test(source, sink = None, sigma = None): # """Hessian eigenvalues of source data # Arguments # --------- # source : array # Input array. # sink : array # Output, if None, a new array is allocated. # sigma : float or None # If not None, a Gaussian filter with std sigma is applied initialliy. # Returns # ------- # sink : array # The three eigenvalues along the first axis for each source. # """ # return _apply_code(code.eigenvalues_test, source, sink, sink_shape_per_pixel = (3,), parameter = None, sigma = sigma); ############################################################################### ### Helpers ############################################################################### def _apply_code(function, source, sink, sink_dtype = None, sink_shape_per_pixel = None, parameter = None, sigma = None): """Helper to apply the core functions""" if source.ndim != 3: raise ValueError('The tubness measure is implemented for 3d data, found %dd!' % source.ndim); if source is sink: raise ValueError("Cannot perform operation in place!") if sink_shape_per_pixel is None: shape_per_pixel = (1,); else: shape_per_pixel = sink_shape_per_pixel; if sink is None: if sink_dtype is None: sink_dtype = float sink = np.zeros(source.shape + shape_per_pixel, dtype = sink_dtype, order = 'F') else: if shape_per_pixel != (1,): if sink.shape != source.shape + shape_per_pixel: raise ValueError('The sink of shape %r does not have expected shape %r!' % (sink.shape, source.shape + shape_per_pixel)); if sink.shape != source.shape + shape_per_pixel: sink = sink.reshape(source.shape + shape_per_pixel) if sink.dtype == bool: s = sink.view('uint8') else: s = sink; sink_stride = io.element_strides(sink)[-1]; if parameter is None: parameter = np.zeros(0); parameter = np.asarray([parameter], dtype = float).flatten(); if sigma is not None: data = ndi.gaussian_filter(np.asarray(source, dtype=float), sigma=sigma); else: data = np.asarray(source, dtype=float); function(source=data, sink=s, sink_stride=sink_stride, parameter=parameter) if sink_shape_per_pixel is None: sink = sink.reshape(sink.shape[:-1]); return sink ############################################################################### ### Tests ############################################################################### def _test(): import numpy as np import ClearMap.Tests.Files as tst import ClearMap.Visualization.Plot3d as p3d import ClearMap.ImageProcessing.Differentiation as dif from importlib import reload reload(dif); # test hessian import scipy.ndimage as ndi sigma = 2; data = np.random.rand(20,20,1); d1 = ndi.gaussian_filter(data[:,:,0], sigma=sigma, order=(1,1)); import ClearMap.ImageProcessing.Differentiation.Hessian as hes d2 = hes.hessian(data, sigma=sigma) clim = (np.min([d1, d2[:,:,0,0,1]]), np.max([d1, d2[:,:,0,0,1]])) import matplotlib.pyplot as plt plt.figure(1); plt.clf(); plt.subplot(1,3,1) plt.imshow(d1, clim=clim) plt.subplot(1,3,2) plt.imshow(d2[:,:,0,0,1], clim=clim) plt.subplot(1,3,3); plt.imshow(d1 - d2[:,:,0,0,1], clim=clim) plt.plot(d1[:]) plt.plot(d2[:,0,0,0,0]) # test hessian eigenvalues source = np.array(tst.source('v')[:7, :7, :7], dtype=float, order = 'F') hessian = dif.hessian(source, sigma=None) # eigenvalues evs_np = np.sort(np.linalg.eigvals(hessian), axis=-1)[...,::-1] evs = dif.hessian_eigenvalues(source, sigma=None) print(np.allclose(evs, evs_np)) # eigenvectors es = dif.hessian_eigensystem(source, sigma=None) evs = es[...,0,:]; evecs = es[...,1:,:] axis = -1 index = list(np.ix_(*[np.arange(i) for i in evs.shape])) index[axis] = evs.argsort(axis) index = tuple(index) evs = evs[index]; evecs = evecs[index] evs_np, evecs_np = np.linalg.eigh(hessian); evecs_np = evecs_np.transpose((0,1,2,4,3)) # numpy gives eigenvectros in the form evec[:,index] axis = -1 index_np = list(np.ix_(*[np.arange(i) for i in evs_np.shape])) index_np[axis] = evs_np.argsort(axis) index_np = tuple(index_np) evs_np = evs_np[index_np] evecs_np = evecs_np[index_np] print(np.allclose(evs, evs_np)) print(np.allclose(np.abs(evecs / evecs_np), 1)) #wrong = np.array(np.where(np.logical_not(np.isclose(np.abs(evecs / evecs_np), 1)))).T #wrong = wrong[::9][:,:3] sink = np.zeros(source.shape + (3,), order = 'C'); dif.hessian_eigenvalues(source, sink, sigma = 1.0); p3d.plot([source] + list(sink.transpose([3,0,1,2]))); sink = np.zeros(source.shape + (3,), order = 'F'); dif.hessian_eigenvalues(source, sink, sigma = 1.0); p3d.plot([source] + list(sink.transpose([3,0,1,2]))); sink = np.zeros(source.shape); dif.lambda123(source, sink, sigma = 1.0); p3d.plot([source, sink]) sink = dif.lambda123(source, sigma = 1.0, threshold = 0.2); p3d.plot([source, sink]); #test = dif.hessian_test(source, sigma=None) # Hessian: #A B C #B D E #C E F #hessian_test = np.array([[test[...,3], test[...,4], test[...,5]],[test[...,4], test[...,6], test[...,7]], [test[...,5], test[...,7], test[...,8]]]) #hessian_test = Hessian.transpose((2,3,4,0,1)) #np.allclose(hessian_test, hessian) # Python implementations of HessianCode for testing #from scipy import ndimage as ndi # #def hessian(source): # """Returns the hessian matrix at each location calculatd via finite differences. # # Arguments # --------- # source : array # Input array. # sink : array # Output, if None, a new array is allocated. # # Returns # ------- # ndarray: # 5d array with the hessian matrix in the first two dimensions # """ # c = 2*source; # h[0,0] = mm[0:-2,1:-1,1:-1] - c + mm[2:,1:-1,1:-1]; # h[1,1] = mm[1:-1,0:-2,1:-1] - c + mm[1:-1,2:,1:-1]; # h[2,2] = mm[1:-1,1:-1,0:-2] - c + mm[1:-1,1:-1,2:]; # # h[0,1] = h[1,0] = (mm[2:,2:,1:-1] - mm[:-2,2:,1:-1] - mm[2:,:-2,1:-1] + mm[:-2,:-2,1:-1]) / 4; # h[0,2] = h[2,0] = (mm[2:,1:-1,2:] - mm[:-2,1:-1,2:] - mm[2:,1:-1,:-2] + mm[:-2,1:-1,:-2]) / 4; # h[1,2] = h[2,1] = (mm[1:-1,2:,2:] - mm[1:-1,:-2,2:] - mm[1:-1,2:,:-2] + mm[1:-1,:-2,:-2]) / 4; # # return h; # #def eigenvalues3D(m): # """Find the coefficients of the characteristic polynomial: # http://en.wikipedia.org/wiki/Eigenvalue_algorithm # # // The matrix looks like: # /* # A B C # B D E # C E F # */ # """ # # A = m[0,0]; # B = m[0,1]; # C = m[0,2]; # D = m[1,1]; # E = m[1,2]; # F = m[2,2]; # # a = -1; # b = A + D + F; # c = B * B + C * C + E * E - A * D - A * F - D * F; # d = A * D * F - A * E * E - B * B * F + 2 * B * C * E - C * C * D; # # third = 0.333333333333333333333333333333333333; # # #Now use the root-finding formula described here: # #http://en.wikipedia.org/wiki/Cubic_equation#oot-finding_formula # # q = (3*a*c - b*b) / (9*a*a); # r = (9*a*b*c - 27*a*a*d - 2*b*b*b) / (54*a*a*a); # # discriminant = q*q*q + r*r; # #print discriminant.shape # # eigenValues = np.zeros((3,) + m.shape[2:]); # #print eigenValues.shape # # #ids = discriminant > 0; # #if np.any(ids): # #should never happen for symmetric matrix # # return None; # # ids = discriminant < 0; # # rr = r[ids]; # # rootThree = 1.7320508075688772935; # innerSize = np.sqrt( rr*rr - discriminant[ids] ); # # ids2 = rr > 0; # innerAngle = np.zeros(rr.shape); # # innerAngle[ids2] = np.arctan(np.sqrt(-discriminant[ids][ids2]) / rr[ids2] ); # # ids2 = np.logical_not(ids2); # innerAngle[ids2] = ( np.pi - np.arctan( np.sqrt(-discriminant[ids][ids2]) / -rr[ids2] ) ); # # # So now s is the cube root of innerSize * e ^ ( innerAngle * i ) # # and t is the cube root of innerSize * e ^ ( - innerAngle * i ) # # stSize = np.power(innerSize,third); # # sAngle = innerAngle / 3; # #tAngle = - innerAngle / 3; # # sPlusT = 2 * stSize * np.cos(sAngle); # # eigenValues[0][ids] = ( sPlusT - (b[ids] / (3*a)) ); # firstPart = - (sPlusT / 2) - (b[ids] / 3*a); # lastPart = - rootThree * stSize * np.sin(sAngle); # # eigenValues[1][ids] = ( firstPart + lastPart ); # eigenValues[2][ids] = ( firstPart - lastPart ); # # # #The discriminant is zero (or very small), # #so the second two solutions are the same: # ids = discriminant == 0; # rr = r[ids]; # # ids2 = rr >= 0; # # # sPlusT = np.zeros(rr.shape); # sPlusT[ids2] = 2 * np.power(rr[ids2],third); # ids2 = np.logical_not(ids2); # sPlusT[ids2] = -2 * np.power(-rr[ids2],third); # # bOver3A = b[ids] / (3 * a); # # eigenValues[0][ids] = ( sPlusT - bOver3A ); # eigenValues[1][ids] = (-sPlusT/2 - bOver3A ); # eigenValues[2][ids] = eigenValues[1][ids]; # # return eigenValues; ## ## #def eigenvalues(m, sigma = 1.0, sort = True): # if sigma is not None: # smoothed = ndi.gaussian_filter(np.asarray(m, dtype = float), sigma=sigma); # else: # smoothed = m; # # h3 = hessian(smoothed); # e3 = eigenvalues3D(h3); # #e3a = np.abs(e3); # e3a = e3; # # if sort: # index = list(np.ix_(*[np.arange(i) for i in e3a.shape])) # index[0] = e3a.argsort(axis = 0) # return e3[index]; # else: # return e3; # # #def tubeness(m, sigma = 1.0): # if sigma is not None: # smoothed = ndi.gaussian_filter(np.asarray(m, dtype = float), sigma=sigma); # else: # smoothed = m; # # h3 = hessian(smoothed); # e3 = eigenvalues3D(h3); # e3a = np.abs(e3); # # index = list(np.ix_(*[np.arange(i) for i in e3a.shape])) # index[0] = e3a.argsort(axis = 0) # e3s = e3[index]; # # tb = np.zeros(m.shape); # ids = np.logical_and(e3s[1] < 0, e3s[2] < 0); # tb[ids] = np.sqrt(e3s[1][ids] * e3s[2][ids]); # # return tb; # # #def frangi(m, alpha = 0.5, beta = 0.5, gamma = 100, sigma = 1.0): # if sigma is not None: # smoothed = ndi.gaussian_filter(np.asarray(m, dtype = float), sigma=sigma); # else: # smoothed = m; # # h3 = hessian(smoothed); # e3 = eigenvalues3D(h3); # e3a = np.abs(e3); # # index = list(np.ix_(*[np.arange(i) for i in e3a.shape])) # index[0] = e3a.argsort(axis = 0) # e3s = e3[index]; # # ra = np.abs(e3s[1]) / np.abs(e3s[2]); # rb = np.abs(e3s[0]) / np.sqrt(np.abs(e3s[1] * e3s[2])); # s = np.sqrt(np.square(e3s[0]) + np.square(e3s[1]) + np.square(e3s[2])) # # plate = 1 - np.exp(- np.square(ra) / (2 * np.square(alpha))); # blob = np.exp(-np.square(rb) / (2 * np.square(beta))); # background = 1 - np.exp(-np.square(s) / (2 * np.square(gamma))); # # f = plate * blob * background; # # ids = np.logical_and(e3s[1] >= 0, e3s[2] >= 0); # f[ids] = 0; # # return f; #