diff --git a/.gitignore b/.gitignore
index 55d415e946c5a7bf10daad0c221d82000c9259cd..e0940c2cf00a46463738a744bbe37c0691bb50d6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,11 +1,16 @@
.idea/
.mypy_cache/
+__pycache__/
__material__/
-__pycache__/
+__runtime__/
__version__/
data/
-runtime/
mprofile_*.dat
+
+brick/processing/frangi_c/obj/
+brick/processing/frangi_c/src/original-formatting/
+brick/processing/frangi_c/vesselness.so
+brick/vesselness.so
diff --git a/connection.py b/brick/component/connection.py
similarity index 93%
rename from connection.py
rename to brick/component/connection.py
index 8942fd51fdc845905aaae5a3b8882ea3311ec374..dc0efe6b093ea9e277b8ea660943a6dad22fb675 100644
--- a/connection.py
+++ b/brick/component/connection.py
@@ -29,10 +29,10 @@
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.
-import dijkstra_1_to_n as dk_
-from extension import extension_t
-from soma import soma_t
-from type import array_t, site_h, site_path_h
+import brick.processing.dijkstra_1_to_n as dk_
+from brick.component.extension import extension_t
+from brick.component.soma import soma_t
+from brick.general.type import array_t, site_h, site_path_h
import itertools as it_
from typing import Callable, Sequence, Tuple
diff --git a/extension.py b/brick/component/extension.py
similarity index 83%
rename from extension.py
rename to brick/component/extension.py
index a84f00bdfa62efbe0753f30abf1e0f61a7b3ae4f..0f8a80b7db0f10efda3b7af970e4640390aad19a 100644
--- a/extension.py
+++ b/brick/component/extension.py
@@ -31,10 +31,10 @@
from __future__ import annotations
-import frangi_3d as fg_
-from glial_cmp import glial_cmp_t
-import map_labeling as ml_
-from type import array_t, site_h
+# import brick.processing.frangi_py.frangi_3d as fg_
+import brick.processing.map_labeling as ml_
+from brick.component.glial_cmp import glial_cmp_t
+from brick.general.type import array_t, site_h
from typing import Optional, Sequence, Tuple
@@ -48,6 +48,35 @@ from scipy import ndimage as im_
min_area_c = 100
+
+import ctypes as ct_
+import os.path as ph_ # Pathlib not necessary
+
+folder, script = ph_.split(__file__)
+if folder.__len__() == 0:
+ folder = "./"
+c_extension = ct_.CDLL(
+ ph_.join(folder, "..", "vesselness.so")
+)
+
+vesselness_C = c_extension.vesselness
+vesselness_C.argtypes = (
+ ct_.c_void_p,
+ ct_.c_int,
+ ct_.c_int,
+ ct_.c_int,
+ ct_.c_void_p,
+ ct_.c_float,
+ ct_.c_float,
+ ct_.c_float,
+ ct_.c_float,
+ ct_.c_float,
+ ct_.c_float,
+)
+vesselness_C.restype = None
+
+
+
class extension_t(glial_cmp_t):
#
# soma_uid: connected to a soma somewhere upstream
@@ -156,17 +185,33 @@ class extension_t(glial_cmp_t):
image, size=2, mode="constant", cval=0.0, origin=0
)
- enhanced_img, scale_map = fg_.FrangiEnhancement(
- preprocessed_img,
- scale_range=(0.1, 3),
- scale_step=1,
- alpha=0.8,
- beta=0.5,
- frangi_c=500,
- bright_on_dark=True,
- in_parallel=in_parallel,
+ preprocessed_img = preprocessed_img.astype(np_.float32)
+ enhanced_img = np_.empty(preprocessed_img.shape, dtype=np_.float32)
+ scale_map = np_.empty(preprocessed_img.shape, dtype=np_.float32)
+
+ vesselness_C(
+ preprocessed_img.ctypes.data,
+ *preprocessed_img.shape,
+ enhanced_img.ctypes.data,
+ 0.1,
+ 3.2,
+ 1.0,
+ 0.8,
+ 0.5,
+ 500.0,
)
+ # enhanced_img, scale_map = fg_.FrangiEnhancement(
+ # preprocessed_img,
+ # scale_range=(0.1, 3),
+ # scale_step=1,
+ # alpha=0.8,
+ # beta=0.5,
+ # frangi_c=500,
+ # bright_on_dark=True,
+ # in_parallel=in_parallel,
+ # )
+
# np_.savez_compressed(
# "./runtime/frangi.npz", enhanced_img=enhanced_img, scale_map=scale_map
# )
diff --git a/glial_cmp.py b/brick/component/glial_cmp.py
similarity index 97%
rename from glial_cmp.py
rename to brick/component/glial_cmp.py
index 10c150387cfcba8f7fa6e9e26d227865ff848b17..a4584e32fd131112361e1b8c99027ed09e05e3b3 100644
--- a/glial_cmp.py
+++ b/brick/component/glial_cmp.py
@@ -31,7 +31,7 @@
from __future__ import annotations
-from type import array_t, np_array_picker_h, py_array_picker_h, site_path_h
+from brick.general.type import array_t, np_array_picker_h, py_array_picker_h, site_path_h
from typing import Dict, List, Optional, Tuple
diff --git a/soma.py b/brick/component/soma.py
similarity index 98%
rename from soma.py
rename to brick/component/soma.py
index 85c2cdf3c3ecee584078a5080c2f8dbeef1ce33d..7bdedbaf9b056bf73a573eec7edf95dd3d7e9978 100644
--- a/soma.py
+++ b/brick/component/soma.py
@@ -31,9 +31,9 @@
from __future__ import annotations
-from glial_cmp import glial_cmp_t
-import map_labeling as ml_
-from type import array_t, py_array_picker_h, site_h
+from brick.component.glial_cmp import glial_cmp_t
+import brick.processing.map_labeling as ml_
+from brick.general.type import array_t, py_array_picker_h, site_h
from collections import namedtuple as namedtuple_t
import sys as sy_
diff --git a/features_ext.py b/brick/feature/features_ext.py
similarity index 100%
rename from features_ext.py
rename to brick/feature/features_ext.py
diff --git a/features_soma.py b/brick/feature/features_soma.py
similarity index 100%
rename from features_soma.py
rename to brick/feature/features_soma.py
diff --git a/feedback.py b/brick/general/feedback.py
similarity index 96%
rename from feedback.py
rename to brick/general/feedback.py
index f077e16f93075863d521cd3ab98efecb6aae69ec..c2a341beb1bdcc746f24cd6ef85fc2f3bbe522e5 100644
--- a/feedback.py
+++ b/brick/general/feedback.py
@@ -1,4 +1,3 @@
-from extension import extension_t
# Copyright CNRS/Inria/UNS
# Contributor(s): Eric Debreuve (since 2019)
#
@@ -30,8 +29,9 @@ from extension import extension_t
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.
-import plot as ot_
-from soma import soma_t
+import brick.processing.plot as ot_
+from brick.component.extension import extension_t
+from brick.component.soma import soma_t
from typing import Sequence, Tuple
diff --git a/type.py b/brick/general/type.py
similarity index 100%
rename from type.py
rename to brick/general/type.py
diff --git a/brick/processing/dijkstra_1_to_n.py b/brick/processing/dijkstra_1_to_n.py
new file mode 120000
index 0000000000000000000000000000000000000000..a9cb6fef877e89f67f892c39597ec194b2d42fe1
--- /dev/null
+++ b/brick/processing/dijkstra_1_to_n.py
@@ -0,0 +1 @@
+../../../../../../brick/def/dijkstra-img/dijkstra_1_to_n.py
\ No newline at end of file
diff --git a/brick/processing/frangi_c/makefile b/brick/processing/frangi_c/makefile
new file mode 100644
index 0000000000000000000000000000000000000000..2f9fb57da2213b22f5278a165102678c0d8da364
--- /dev/null
+++ b/brick/processing/frangi_c/makefile
@@ -0,0 +1,43 @@
+SRC_PATH = ./src
+SRCS = $(SRC_PATH)/detector.cpp $(SRC_PATH)/processing.cpp $(SRC_PATH)/kernel3D.cpp $(SRC_PATH)/main.cpp
+TARGET = vesselness
+
+OBJ_PATH = ./obj
+
+CC = g++ -std=c++0x
+RM = rm -f
+
+INCLUDES += -I $(SRC_PATH)
+
+# -g adds debugging information to the executable file
+# -Wall turns on most, but not all, compiler warnings
+# -pg adds profiling information to the executable file
+# -m64 compile for 64-bit
+# -Wno-comment disable warnings about multiline comments
+# CFLAGS = -Wall -Wno-reorder -m64 -fopenmp -Wno-comment -O2 -Wl,--no-as-needed
+CFLAGS = -Wall -O2 -fopenmp
+CFLAGS +=`pkg-config --cflags /usr/lib/pkgconfig/opencv4.pc`
+
+# LFLAGS += -Wl,--no-as-needed -pthread
+LFLAGS += -O2 -pthread -fopenmp -lm
+LFLAGS += `pkg-config --libs /usr/lib/pkgconfig/opencv4.pc`
+
+OBJS = $(SRCS:$(SRC_PATH)/%.cpp=$(OBJ_PATH)/%.o)
+
+all: $(OBJS)
+ mkdir -p $(OBJ_PATH)
+ $(CC) $(OBJS) -o $(TARGET) $(LFLAGS) $(PY_LFLAGS) $(COVERAGE) $(LIBS)
+
+py: PY_CFLAGS = -fpic
+py: PY_LFLAGS = -shared
+py: all
+ mv $(TARGET) $(TARGET).so
+
+cov: COVERAGE = -fprofile-arcs -ftest-coverage
+cov: all
+
+$(OBJ_PATH)/%.o: $(SRC_PATH)/%.cpp
+ $(CC) -c $< -o $@ $(CFLAGS) $(PY_CFLAGS) $(COVERAGE) $(INCLUDES)
+
+clean:
+ $(RM) $(OBJ_PATH)/*.o $(TARGET)
diff --git a/brick/processing/frangi_c/src/__coverage__.txt b/brick/processing/frangi_c/src/__coverage__.txt
new file mode 100644
index 0000000000000000000000000000000000000000..968961f44f3355c3b6f84129d2d636804d82ea85
--- /dev/null
+++ b/brick/processing/frangi_c/src/__coverage__.txt
@@ -0,0 +1,7 @@
+gcc -Wall -fprofile-arcs -ftest-coverage cov.c
+
+Then run executable
+
+gcov cov.c
+
+Then check: cov.c.gcov
diff --git a/brick/processing/frangi_c/src/__formatting__.txt b/brick/processing/frangi_c/src/__formatting__.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4a0ec6418bbd083730d5e247fae88b3ac2da3211
--- /dev/null
+++ b/brick/processing/frangi_c/src/__formatting__.txt
@@ -0,0 +1 @@
+astyle --style=stroustrup --recursive src/*.cpp,*.h
diff --git a/brick/processing/frangi_c/src/data3D.h b/brick/processing/frangi_c/src/data3D.h
new file mode 100644
index 0000000000000000000000000000000000000000..9eb64e5bff691d3ad0aa8bac768d9d0f7dd82b6c
--- /dev/null
+++ b/brick/processing/frangi_c/src/data3D.h
@@ -0,0 +1,338 @@
+#pragma once
+
+#include "type_info.h"
+
+#include <iostream>
+#include <opencv2/core/core.hpp>
+#include <opencv2/highgui/highgui.hpp>
+#include <typeinfo>
+#include <string.h>
+#include "smart_assert.h"
+
+template<typename T>
+class Data3D {
+public:
+ Data3D( const cv::Vec3i& n_size = cv::Vec3i(0,0,0));
+ Data3D( const int width, const int height, const int depth );
+ Data3D( const cv::Vec3i& n_size, const T& value );
+ template <class T2>
+ Data3D( const Data3D<T2>& src );
+ const Data3D<T>& operator=( const Data3D<T>& src );
+
+ virtual ~Data3D(void) { }
+
+ inline void reset( void )
+ {
+ memset( _mat.data, 0, _size_total * sizeof(T) );
+ }
+ inline void reset( const cv::Vec3i& n_size )
+ {
+ _size = n_size;
+ long size_slice = _size[0] * _size[1];
+ _size_total = size_slice * _size[2];
+ _mat = cv::Mat_<T>( _size[2], size_slice );
+ memset( _mat.data, 0, _size_total * sizeof(T) );
+ }
+ inline void reset( const int width, const int height, const int depth )
+ {
+ _size = cv::Vec3i(width, height, depth);
+ long size_slice = _size[0] * _size[1];
+ _size_total = size_slice * _size[2];
+ _mat = cv::Mat_<T>( _size[2], size_slice );
+ memset( _mat.data, 0, _size_total * sizeof(T) );
+ }
+ void reset( const cv::Vec3i& n_size, const T& value )
+ {
+ resize( n_size );
+ for( cv::MatIterator_<T> it=_mat.begin(); it<_mat.end(); it++ )
+ {
+ (*it) = value;
+ }
+ }
+
+ void resize( const cv::Vec3i& n_size )
+ {
+ _size = n_size;
+ long size_slice = _size[0] * _size[1];
+ _size_total = size_slice * _size[2];
+ _mat = cv::Mat_<T>( _size[2], size_slice );
+ }
+
+ inline const int& SX(void) const
+ {
+ return _size[0];
+ }
+ inline const int& SY(void) const
+ {
+ return _size[1];
+ }
+ inline const int& SZ(void) const
+ {
+ return _size[2];
+ }
+ inline const int& get_width(void) const
+ {
+ return _size[0];
+ }
+ inline const int& get_height(void) const
+ {
+ return _size[1];
+ }
+ inline const int& get_depth(void) const
+ {
+ return _size[2];
+ }
+ inline const int& get_size_x(void) const
+ {
+ return _size[0];
+ }
+ inline const int& get_size_y(void) const
+ {
+ return _size[1];
+ }
+ inline const int& get_size_z(void) const
+ {
+ return _size[2];
+ }
+ inline const int& get_size(int i) const
+ {
+ return _size[i];
+ }
+ inline const long& get_size_total(void) const
+ {
+ return _size_total;
+ }
+ inline const cv::Vec3i& get_size(void) const
+ {
+ return _size;
+ }
+ inline const bool is_empty( void ) const
+ {
+ return get_size_total()==0;
+ }
+
+ virtual const inline T& at( const int& i ) const
+ {
+ return _mat( i );
+ }
+ virtual inline T& at( const int& i )
+ {
+ return _mat( i );
+ }
+ virtual inline const T& at( const int& x, const int& y, const int& z ) const
+ {
+ return _mat( z, y * _size[0] + x );
+ }
+ virtual inline T& at( const int& x, const int& y, const int& z )
+ {
+ return _mat( z, y * _size[0] + x );
+ }
+ virtual inline const T& at( const cv::Vec3i& pos ) const
+ {
+ return at( pos[0], pos[1], pos[2] );
+ }
+ virtual inline T& at( const cv::Vec3i& pos )
+ {
+ return at( pos[0], pos[1], pos[2] );
+ }
+
+ inline cv::Mat_<T>& getMat()
+ {
+ return _mat;
+ }
+ inline const cv::Mat_<T>& getMat() const
+ {
+ return _mat;
+ }
+ inline cv::Mat_<T> getMat( int slice ) const
+ {
+ return _mat.row(slice).reshape( 0, get_height() ).clone();
+ }
+
+ void SetFromArray( const T* array );
+ void CopyToArray( T* array );
+
+ template<typename T1, typename T2>
+ friend const Data3D<T1>& operator*=( Data3D<T1>& left, const T2& right );
+ template<typename T1, typename T2, typename T3>
+ friend bool subtract3D( const Data3D<T1>& src1, const Data3D<T2>& src2, Data3D<T3>& dst );
+ template<typename T1, typename T2, typename T3>
+ friend bool multiply3D( const Data3D<T1>& src1, const Data3D<T2>& src2, Data3D<T3>& dst );
+
+ cv::Vec<T, 2> get_min_max_value() const
+ {
+ cv::Vec<T, 2> min_max;
+ cv::Point minLoc, maxLoc;
+ cv::minMaxLoc( _mat, NULL, NULL, &minLoc, &maxLoc);
+ min_max[0] = _mat( minLoc );
+ min_max[1] = _mat( maxLoc );
+ return min_max;
+ }
+
+public:
+ /* copy data of dimension i to a new Image3D structure
+ e.g. this is useful when trying to visualize vesselness which is
+ a multidimensional data structure
+ Yuchen: p.s. I have no idea how to move the funciton body out the
+ class definition nicely. Let me figure it later maybe. */
+ template<typename T2>
+ void copyDimTo( Data3D<T2>& dst, int dim ) const
+ {
+ dst.reset( this->get_size() );
+ int x, y, z;
+ for( z=0; z<get_size_z(); z++ ) for( y=0; y<get_size_y(); y++ ) for( x=0; x<get_size_x(); x++ ) {
+ dst.at(x,y,z) = this->at(x,y,z)[dim];
+ }
+ }
+
+
+ // return true if a index is valide for the data
+ inline bool isValid( const int& x, const int& y, const int& z ) const
+ {
+ return ( x>=0 && x<get_size_x() &&
+ y>=0 && y<get_size_y() &&
+ z>=0 && z<get_size_z() );
+ }
+ inline bool isValid( const cv::Vec3i& v ) const
+ {
+ return isValid( v[0], v[1], v[2] );
+ }
+
+ inline const T* getData(void) const
+ {
+ return (T*) getMat().data;
+ }
+
+ inline T* getData(void)
+ {
+ return (T*) getMat().data;
+ }
+
+protected:
+ cv::Vec3i _size;
+ long _size_total;
+ cv::Mat_<T> _mat;
+};
+
+
+// Default Constructor
+template <typename T>
+Data3D<T>::Data3D( const cv::Vec3i& n_size )
+{
+ reset(n_size);
+}
+template <typename T>
+Data3D<T>::Data3D( const int width, const int height, const int depth )
+{
+ reset(width, height, depth);
+}
+
+// Constructor with size and value
+template <typename T>
+Data3D<T>::Data3D( const cv::Vec3i& n_size, const T& value )
+{
+ reset(n_size, value);
+}
+
+
+// Copy Constructor
+template<typename T>
+template<typename T2>
+Data3D<T>::Data3D( const Data3D<T2>& src )
+{
+ // resize
+ this->resize( src.get_size() );
+ // copy the data over
+ cv::MatIterator_<T> this_it;
+ cv::MatConstIterator_<T2> src_it;
+ for( this_it = this->getMat().begin(), src_it = src.getMat().begin();
+ this_it < this->getMat().end(), src_it < src.getMat().end();
+ this_it++, src_it++ ) {
+ // Yuchen: The convertion from T2 to T may be unsafe
+ *(this_it) = T( *(src_it) );
+ }
+}
+
+
+template<typename T>
+const Data3D<T>& Data3D<T>::operator=( const Data3D<T>& src )
+{
+ // resize
+ this->resize( src.get_size() );
+
+ // copy the data over
+ cv::MatIterator_<T> this_it;
+ cv::MatConstIterator_<T> src_it;
+ for( this_it = this->getMat().begin(), src_it = src.getMat().begin();
+ this_it < this->getMat().end(), src_it < src.getMat().end();
+ this_it++, src_it++ ) {
+ *(this_it) = *(src_it);
+ }
+
+ return *this;
+}
+
+
+template <typename T>
+void Data3D<T>::SetFromArray( const T* array )
+{
+ memcpy(_mat.data, array, _size_total * sizeof(T));
+}
+
+template <typename T>
+void Data3D<T>::CopyToArray( T* array )
+{
+ memcpy(array, _mat.data, _size_total * sizeof(T));
+}
+
+//////////////////////////////////////////////////////////////////////
+// Operator Overloading and other friend functions
+
+template<typename T, typename T2>
+const Data3D<T>& operator*=( Data3D<T>& left, const T2& right )
+{
+ left._mat*=right;
+ return left;
+}
+
+template<typename T1, typename T2, typename T3>
+bool subtract3D( const Data3D<T1>& src1, const Data3D<T2>& src2, Data3D<T3>& dst )
+{
+ smart_return( src1.get_size()==src2.get_size(), "Sizes should match.", false);
+
+ if( dst.get_size() != src1.get_size() ) {
+ dst.reset( src1.get_size() );
+ }
+
+ int x, y, z;
+ for( z=0; z<src1.get_size_z(); z++ ) {
+ for( y=0; y<src1.get_size_y(); y++ ) {
+ for( x=0; x<src1.get_size_x(); x++ ) {
+ dst.at(x,y,z) = T3(src1.at(x,y,z) - src2.at(x,y,z));
+ }
+ }
+ }
+ return true;
+}
+
+
+template<typename T1, typename T2, typename T3>
+bool multiply3D( const Data3D<T1>& src1, const Data3D<T2>& src2, Data3D<T3>& dst )
+{
+ smart_return_false( src1.get_size()==src2.get_size(),
+ "Source sizes are supposed to be cv::Matched.");
+
+ if( dst.get_size() != src1.get_size() ) {
+ dst.reset( src1.get_size() );
+ }
+
+ int x, y, z;
+ for( z=0; z<src1.get_size_z(); z++ ) {
+ for( y=0; y<src1.get_size_y(); y++ ) {
+ for( x=0; x<src1.get_size_x(); x++ ) {
+ dst.at(x,y,z) = src1.at(x,y,z) * src2.at(x,y,z);
+ }
+ }
+ }
+ return true;
+}
diff --git a/brick/processing/frangi_c/src/detector.cpp b/brick/processing/frangi_c/src/detector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d68f8de970fca36af79a3c47be5ebb35f5b4a24
--- /dev/null
+++ b/brick/processing/frangi_c/src/detector.cpp
@@ -0,0 +1,181 @@
+#include "detector.h"
+#include "data3D.h"
+#include "processing_code.h"
+#include "types.h"
+#include "eigen.cpp" // Otherwise the template mechanism does not work. TODO: solve this
+
+using namespace std;
+
+Vesselness hessien_thread_func( const Data3D<float>& src,
+ const int& x, const int& y, const int& z,
+ const float& alpha, const float& beta, const float& gamma )
+{
+
+ ////////////////////////////////////////////////////////////////////
+ // The following are being computed in this function
+ // 1) derivative of images;
+ // 2) Hessian matrix;
+ // 3) Eigenvalue decomposition;
+ // 4) vesselness measure.
+
+ // 1) derivative of the image
+ float im_dx2 = -2.0f * src.at(x,y,z) + src.at(x-1,y,z) + src.at(x+1,y,z);
+ float im_dy2 = -2.0f * src.at(x,y,z) + src.at(x,y-1,z) + src.at(x,y+1,z);
+ float im_dz2 = -2.0f * src.at(x,y,z) + src.at(x,y,z-1) + src.at(x,y,z+1);
+ // 1) derivative of the image (alternative approach, the one above is more accurate)
+ //float im_dx2 = -0.5f * src.at(x,y,z) + 0.25f * src.at(x-2,y,z) + 0.25f * src.at(x+2,y,z);
+ //float im_dy2 = -0.5f * src.at(x,y,z) + 0.25f * src.at(x,y-2,z) + 0.25f * src.at(x,y+2,z);
+ //float im_dz2 = -0.5f * src.at(x,y,z) + 0.25f * src.at(x,y,z-2) + 0.25f * src.at(x,y,z+2);
+
+ float im_dxdy = (
+ + src.at(x-1, y-1, z)
+ + src.at(x+1, y+1, z)
+ - src.at(x-1, y+1, z)
+ - src.at(x+1, y-1, z) ) * 0.25f;
+ float im_dxdz = (
+ + src.at(x-1, y, z-1)
+ + src.at(x+1, y, z+1)
+ - src.at(x+1, y, z-1)
+ - src.at(x-1, y, z+1) ) * 0.25f;
+ float im_dydz = (
+ + src.at(x, y-1, z-1)
+ + src.at(x, y+1, z+1)
+ - src.at(x, y+1, z-1)
+ - src.at(x, y-1, z+1) ) * 0.25f;
+
+ // 3) Eigenvalue decomposition
+ // Reference: http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices
+ // Given a real symmetric 3x3 matrix Hessian, compute the eigenvalues
+ const float Hessian[6] = {im_dx2, im_dxdy, im_dxdz, im_dy2, im_dydz, im_dz2};
+
+ float eigenvalues[3];
+ float eigenvectors[3][3];
+ eigen_decomp( Hessian, eigenvalues, eigenvectors );
+
+ // order eigenvalues so that |lambda1| < |lambda2| < |lambda3|
+ int i=0, j=1, k=2;
+// if( abs(eigenvalues[i]) > abs(eigenvalues[j]) ) std::swap( i, j );
+// if( abs(eigenvalues[i]) > abs(eigenvalues[k]) ) std::swap( i, k );
+// if( abs(eigenvalues[j]) > abs(eigenvalues[k]) ) std::swap( j, k );
+
+ // 4) vesselness measure
+ Vesselness vn;
+ // vesselness value
+ if( eigenvalues[j] > 0 || eigenvalues[k] > 0 ) {
+ vn.rsp = 0.0f;
+ }
+ else {
+ float lmd1 = abs( eigenvalues[i] );
+ float lmd2 = abs( eigenvalues[j] );
+ float lmd3 = abs( eigenvalues[k] );
+
+ float A = (lmd3>1e-5) ? lmd2/lmd3 : 0;
+ float B = (lmd2*lmd3>1e-5) ? lmd1 / sqrt( lmd2*lmd3 ) : 0;
+ float S = sqrt( lmd1*lmd1 + lmd2*lmd2 + lmd3*lmd3 );
+ vn.rsp = ( 1.0f-exp(-A*A/alpha) )* exp( B*B/beta ) * ( 1-exp(-S*S/gamma) );
+ }
+
+ // orientation of vesselness is corresponding to the eigenvector of the
+ // smallest eigenvalue
+ for( int d=0; d<3; d++ ) {
+ vn.dir[d] = eigenvectors[i][d];
+ }
+
+ return vn;
+}
+
+
+
+bool VesselDetector::hessien( const Data3D<float>& src, Data3D<Vesselness>& dst,
+ int ksize, float sigma,
+ float alpha, float beta, float gamma )
+{
+ if( ksize!=0 && sigma<1e-3 ) { // sigma is not set
+ sigma = 0.15f * ksize + 0.35f;
+ }
+ else if ( ksize==0 && sigma>1e-3 ) { // size is not set
+ ksize = int( 6*sigma+1 );
+ // make sure size is an odd number
+ if ( ksize%2==0 ) ksize++;
+ }
+ else {
+ std::cerr << "At lease size or sigma has to be set." << std::endl;
+ return false;
+ }
+
+ Data3D<float> im_blur;
+ bool flag = ImageProcessing::GaussianBlur3D( src, im_blur, ksize, sigma );
+ smart_return( flag, "Gaussian Blur Failed.", false );
+
+ im_blur *= std::pow(sigma, 1.2); //Normalizing for different scale
+
+ dst.reset( im_blur.get_size(), Vesselness(0.0f) );
+
+ #pragma omp parallel
+ {
+ #pragma omp for
+ for( int z = 1; z < src.get_size_z()-1; z++ ) {
+ for( int y = 1; y < src.get_size_y()-1; y++ ) {
+ for( int x = 1; x < src.get_size_x()-1; x++ ) {
+ dst.at(x, y, z) = hessien_thread_func( im_blur, x, y, z, alpha, beta, gamma);
+ }
+ }
+ }
+ }
+
+ return true;
+}
+
+
+
+int VesselDetector::compute_vesselness(
+ const Data3D<float>& src,
+ Data3D<Vesselness_Sig>& dst,
+ float sigma_from, float sigma_to, float sigma_step,
+ float alpha, float beta, float gamma )
+{
+ std::cout << "Computing Vesselness, it will take a while... " << std::endl;
+ std::cout << "Vesselness will be computed from sigma = " << sigma_from << " to sigma = " << sigma_to << std::endl;
+
+ dst.reset( src.get_size() ); // reszie data, and it will also be clear to zero
+
+ // Error for input parameters
+ smart_return( sigma_from < sigma_to,
+ "sigma_from should be smaller than sigma_to ", 0 );
+ smart_return( sigma_step > 0,
+ "sigma_step should be greater than 0 ", 0 );
+
+ int x,y,z;
+ float max_sigma = sigma_from;
+ float min_sigma = sigma_to;
+
+ Data3D<Vesselness> temp;
+ for( float sigma = sigma_from; sigma < sigma_to; sigma += sigma_step ) {
+ cout << '\r' << "Vesselness for sigma = " << sigma << " " << "\b\b\b\b";
+ cout.flush();
+
+ VesselDetector::hessien( src, temp, 0, sigma, alpha, beta, gamma );
+
+ const int margin = 1;
+ for( z=margin; z<src.get_size_z()-margin; z++ ) {
+ for( y=margin; y<src.get_size_y()-margin; y++ ) {
+ for( x=margin; x<src.get_size_x()-margin; x++ ) {
+ // Update vesselness if the new response is greater
+ if( dst.at(x, y, z).rsp < temp.at(x, y, z).rsp ) {
+ dst.at(x, y, z).rsp = temp.at(x, y, z).rsp;
+ dst.at(x, y, z).dir = temp.at(x, y, z).dir;
+ dst.at(x, y, z).sigma = sigma;
+ max_sigma = std::max( sigma, max_sigma );
+ min_sigma = std::min( sigma, min_sigma );
+ }
+ }
+ }
+ }
+ }
+
+ std::cout << std::endl << "The range for sigma is [";
+ std::cout << min_sigma << ", " << max_sigma << "]" << std::endl;
+ std::cout << "Done. " << std::endl << std::endl;
+
+ return 0;
+}
diff --git a/brick/processing/frangi_c/src/detector.h b/brick/processing/frangi_c/src/detector.h
new file mode 100644
index 0000000000000000000000000000000000000000..e08e6eb644f6c9f5f6e71fb32a76e222f36461c7
--- /dev/null
+++ b/brick/processing/frangi_c/src/detector.h
@@ -0,0 +1,31 @@
+#ifndef VESSELDETECTOR_H
+#define VESSELDETECTOR_H
+
+#include "types.h"
+
+template<typename T> class Data3D;
+
+class Vesselness;
+class Vesselness_Sig;
+class Vesselness_Nor;
+class Vesselness_All;
+
+namespace VesselDetector {
+bool hessien(
+ const Data3D<float>& src, Data3D<Vesselness>& dst,
+ int ksize, float sigma,
+ float alpha, float beta, float gamma );
+
+int compute_vesselness(
+ const Data3D<float>& src,
+ Data3D<Vesselness_Sig>& dst,
+ float sigma_from, float sigma_to, float sigma_step,
+ float alpha = 1.0e-1f,
+ float beta = 5.0e0f,
+ float gamma = 3.5e5f );
+};
+
+namespace VD = VesselDetector;
+namespace VF = VesselDetector;
+
+#endif
diff --git a/brick/processing/frangi_c/src/eigen.cpp b/brick/processing/frangi_c/src/eigen.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..71ab97db79e43db9dd132037b0e045558f42d8ed
--- /dev/null
+++ b/brick/processing/frangi_c/src/eigen.cpp
@@ -0,0 +1,576 @@
+#include "eigen.h"
+
+/* Compute the eigen value decomposition of a 3 by 3 symetric matrix
+Input:
+A = [ A11 A12 A13 = [ A[0] A[1] A[2]
+A21 A22 A23 0 A[3] A[4]
+A31 A32 A33 ]; 0 0 A[5] ];
+Output:
+eigenvalues, eigenvectors
+
+Reference: http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices
+*/
+#include <cmath>
+#include <assert.h>
+#include <algorithm>
+#include <iostream>
+// #include <opencv2/core/types_c.h>
+// #include <opencv2/core/core_c.h>
+
+
+#if !defined(SmallEpsilon) && !defined(BigEpsilon)
+#define SmallEpsilon 1.0e-18
+#define BigEpsilon 1.0e-4
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#include "math.h"
+#ifdef MAX
+#undef MAX
+#endif
+#define MAX(a, b) ((a)>(b)?(a):(b))
+#define n 3
+
+// address/from/source: https://github.com/fethallah/neurons/blob/master/frangi_filter_version2a/eig3volume.c
+
+/* Eigen decomposition code for symmetric 3x3 matrices, copied from the public
+ * domain Java Matrix library JAMA. */
+static float hypot2(float x, float y) { return (float) sqrt(x*x+y*y); }
+
+__inline float absd(float val){ if(val>0){ return val;} else { return -val;} };
+
+/* Symmetric Householder reduction to tridiagonal form. */
+static void tred2(float V[n][n], float d[n], float e[n]) {
+
+/* This is derived from the Algol procedures tred2 by */
+/* Bowdler, Martin, Reinsch, and Wilkinson, Handbook for */
+/* Auto. Comp., Vol.ii-Linear Algebra, and the corresponding */
+/* Fortran subroutine in EISPACK. */
+
+// std::cout << "Entering tred2" << std::endl;
+
+ int i, j, k;
+ float scale;
+ float f, g, h;
+ float hh;
+ for (j = 0; j < n; j++) {d[j] = V[n-1][j]; }
+
+ /* Householder reduction to tridiagonal form. */
+
+ for (i = n-1; i > 0; i--) {
+ /* Scale to avoid under/overflow. */
+ scale = 0.0;
+ h = 0.0;
+ for (k = 0; k < i; k++) { scale = scale + fabs(d[k]); }
+ if (scale == 0.0) {
+ e[i] = d[i-1];
+ for (j = 0; j < i; j++) { d[j] = V[i-1][j]; V[i][j] = 0.0; V[j][i] = 0.0; }
+ } else {
+
+ /* Generate Householder vector. */
+
+ for (k = 0; k < i; k++) { d[k] /= scale; h += d[k] * d[k]; }
+ f = d[i-1];
+ g = sqrt(h);
+ if (f > 0) { g = -g; }
+ e[i] = scale * g;
+ h = h - f * g;
+ d[i-1] = f - g;
+ for (j = 0; j < i; j++) { e[j] = 0.0; }
+
+ /* Apply similarity transformation to remaining columns. */
+
+ for (j = 0; j < i; j++) {
+ f = d[j];
+ V[j][i] = f;
+ g = e[j] + V[j][j] * f;
+ for (k = j+1; k <= i-1; k++) { g += V[k][j] * d[k]; e[k] += V[k][j] * f; }
+ e[j] = g;
+ }
+ f = 0.0;
+ for (j = 0; j < i; j++) { e[j] /= h; f += e[j] * d[j]; }
+ hh = f / (h + h);
+ for (j = 0; j < i; j++) { e[j] -= hh * d[j]; }
+ for (j = 0; j < i; j++) {
+ f = d[j]; g = e[j];
+ for (k = j; k <= i-1; k++) { V[k][j] -= (f * e[k] + g * d[k]); }
+ d[j] = V[i-1][j];
+ V[i][j] = 0.0;
+ }
+ }
+ d[i] = h;
+ }
+
+ /* Accumulate transformations. */
+
+ for (i = 0; i < n-1; i++) {
+ V[n-1][i] = V[i][i];
+ V[i][i] = 1.0;
+ h = d[i+1];
+ if (h != 0.0) {
+ for (k = 0; k <= i; k++) { d[k] = V[k][i+1] / h;}
+ for (j = 0; j <= i; j++) {
+ g = 0.0;
+ for (k = 0; k <= i; k++) { g += V[k][i+1] * V[k][j]; }
+ for (k = 0; k <= i; k++) { V[k][j] -= g * d[k]; }
+ }
+ }
+ for (k = 0; k <= i; k++) { V[k][i+1] = 0.0;}
+ }
+ for (j = 0; j < n; j++) { d[j] = V[n-1][j]; V[n-1][j] = 0.0; }
+ V[n-1][n-1] = 1.0;
+ e[0] = 0.0;
+
+// std::cout << "Leaving tred2" << std::endl;
+}
+
+/* Symmetric tridiagonal QL algorithm. */
+static void tql2(float V[n][n], float d[n], float e[n]) {
+
+/* This is derived from the Algol procedures tql2, by */
+/* Bowdler, Martin, Reinsch, and Wilkinson, Handbook for */
+/* Auto. Comp., Vol.ii-Linear Algebra, and the corresponding */
+/* Fortran subroutine in EISPACK. */
+
+// std::cout << "Entering tql2" << std::endl;
+
+ int i, j, k, l, m;
+ float f;
+ float tst1;
+ float eps;
+ int iter;
+ float g, p, r;
+ float dl1, h, c, c2, c3, el1, s, s2;
+
+ for (i = 1; i < n; i++) { e[i-1] = e[i]; }
+ e[n-1] = 0.0;
+
+ f = 0.0;
+ tst1 = 0.0;
+ eps = pow(2.0, -52.0);
+ for (l = 0; l < n; l++) {
+
+ /* Find small subdiagonal element */
+
+ tst1 = MAX(tst1, fabs(d[l]) + fabs(e[l]));
+ m = l;
+ while (m < n) {
+ if (fabs(e[m]) <= eps*tst1) { break; }
+ m++;
+ }
+
+ /* If m == l, d[l] is an eigenvalue, */
+ /* otherwise, iterate. */
+
+ if (m > l) {
+ iter = 0;
+ do {
+ iter = iter + 1; /* (Could check iteration count here.) */
+ /* Compute implicit shift */
+ g = d[l];
+ p = (d[l+1] - g) / (2.0 * e[l]);
+ r = hypot2(p, 1.0);
+ if (p < 0) { r = -r; }
+ d[l] = e[l] / (p + r);
+ d[l+1] = e[l] * (p + r);
+ dl1 = d[l+1];
+ h = g - d[l];
+ for (i = l+2; i < n; i++) { d[i] -= h; }
+ f = f + h;
+ /* Implicit QL transformation. */
+ p = d[m]; c = 1.0; c2 = c; c3 = c;
+ el1 = e[l+1]; s = 0.0; s2 = 0.0;
+ for (i = m-1; i >= l; i--) {
+ c3 = c2;
+ c2 = c;
+ s2 = s;
+ g = c * e[i];
+ h = c * p;
+ r = hypot2(p, e[i]);
+ e[i+1] = s * r;
+ s = e[i] / r;
+ c = p / r;
+ p = c * d[i] - s * g;
+ d[i+1] = h + s * (c * g + s * d[i]);
+ /* Accumulate transformation. */
+ for (k = 0; k < n; k++) {
+ h = V[k][i+1];
+ V[k][i+1] = s * V[k][i] + c * h;
+ V[k][i] = c * V[k][i] - s * h;
+ }
+ }
+ p = -s * s2 * c3 * el1 * e[l] / dl1;
+ e[l] = s * p;
+ d[l] = c * p;
+
+ /* Check for convergence. */
+ } while (fabs(e[l]) > eps*tst1);
+ }
+ d[l] = d[l] + f;
+ e[l] = 0.0;
+ }
+
+ /* Sort eigenvalues and corresponding vectors. */
+ for (i = 0; i < n-1; i++) {
+ k = i;
+ p = d[i];
+ for (j = i+1; j < n; j++) {
+ if (d[j] < p) {
+ k = j;
+ p = d[j];
+ }
+ }
+ if (k != i) {
+ d[k] = d[i];
+ d[i] = p;
+ for (j = 0; j < n; j++) {
+ p = V[j][i];
+ V[j][i] = V[j][k];
+ V[j][k] = p;
+ }
+ }
+ }
+
+// std::cout << "Leaving tql2" << std::endl;
+}
+
+void eigen_decomposition(float A[n][n], float V[n][n], float d[n]) {
+ float e[n];
+ float da[3];
+ float dt, dat;
+ float vet[3];
+ int i, j;
+
+// std::cout << "Entering egien_decompositiom" << std::endl;
+
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < n; j++) {
+ V[i][j] = A[i][j];
+ }
+ }
+ tred2(V, d, e);
+ tql2(V, d, e);
+
+ /* Sort the eigen values and vectors by abs eigen value */
+ da[0]=absd(d[0]); da[1]=absd(d[1]); da[2]=absd(d[2]);
+ if((da[0]>=da[1])&&(da[0]>da[2]))
+ {
+ dt=d[2]; dat=da[2]; vet[0]=V[0][2]; vet[1]=V[1][2]; vet[2]=V[2][2];
+ d[2]=d[0]; da[2]=da[0]; V[0][2] = V[0][0]; V[1][2] = V[1][0]; V[2][2] = V[2][0];
+ d[0]=dt; da[0]=dat; V[0][0] = vet[0]; V[1][0] = vet[1]; V[2][0] = vet[2];
+ }
+ else if((da[1]>=da[0])&&(da[1]>da[2]))
+ {
+ dt=d[2]; dat=da[2]; vet[0]=V[0][2]; vet[1]=V[1][2]; vet[2]=V[2][2];
+ d[2]=d[1]; da[2]=da[1]; V[0][2] = V[0][1]; V[1][2] = V[1][1]; V[2][2] = V[2][1];
+ d[1]=dt; da[1]=dat; V[0][1] = vet[0]; V[1][1] = vet[1]; V[2][1] = vet[2];
+ }
+ if(da[0]>da[1])
+ {
+ dt=d[1]; dat=da[1]; vet[0]=V[0][1]; vet[1]=V[1][1]; vet[2]=V[2][1];
+ d[1]=d[0]; da[1]=da[0]; V[0][1] = V[0][0]; V[1][1] = V[1][0]; V[2][1] = V[2][0];
+ d[0]=dt; da[0]=dat; V[0][0] = vet[0]; V[1][0] = vet[1]; V[2][0] = vet[2];
+ }
+
+// std::cout << "Leaving egien_decompositiom" << std::endl;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*
+template<typename T=float>
+inline void normalize( T v[3] )
+{
+ T length = 0;
+ for( int i=0; i<3; i++ ) {
+ length += v[i] * v[i];
+ }
+// assert( length > SmallEpsilon && "Cannot normalized zero vector" );
+ if ( length <= SmallEpsilon )
+ return;
+
+ length = sqrt( length );
+ for( int i=0; i<3; i++ ) {
+ v[i] /= length;
+ }
+}
+
+template<typename T=float>
+inline void cross_product( const T v1[3], const T v2[3], T v3[3] )
+{
+ v3[0] = v1[1] * v2[2] - v1[2] * v2[1];
+ v3[1] = v1[2] * v2[0] - v1[0] * v2[2];
+ v3[2] = v1[0] * v2[1] - v1[1] * v2[0];
+ normalize( v3 );
+}
+
+template<typename T=float>
+void normal_vectors( const T v1[3], T v2[3], T v3[3] )
+{
+ // an arbitrary combination of the following two vectors
+ // is a normal to v1
+ // alpha * (-v1[1], v1[0], 0) + beta * (0, -v1[2], v[1])
+
+ if( abs(v1[0])>BigEpsilon || abs(v1[1])>BigEpsilon ) {
+ v2[0] = -v1[1];
+ v2[1] = v1[0];
+ v2[2] = 0;
+ }
+ else {
+ v2[0] = 0;
+ v2[1] = -v1[2];
+ v2[2] = v1[1];
+ }
+
+ normalize( v2 );
+
+ cross_product( v1, v2, v3 );
+}
+*/
+
+template<typename T=float>
+void eigen_decomp( const T A[6], T eigenvalues[3], T eigenvectors[3][3] )
+{
+ float Asym[n][n];
+ Asym[0][0] = A[0]; Asym[0][1] = A[1]; Asym[0][2] = A[2];
+ Asym[1][0] = A[1]; Asym[1][1] = A[3]; Asym[1][2] = A[4];
+ Asym[2][0] = A[2]; Asym[2][1] = A[4]; Asym[2][2] = A[5];
+
+ /* Function below sorts the eigen values and vectors by abs eigen value */
+ eigen_decomposition(Asym, eigenvectors, eigenvalues);
+
+ return;
+
+// static CvMat mat, *eigenVec, *eigenVal;
+// static float Ain1D[9];
+//
+// int cnt = 0;
+// Ain1D[cnt++] = A[0];
+// Ain1D[cnt++] = A[1];
+// Ain1D[cnt++] = A[2];
+// Ain1D[cnt++] = A[1];
+// Ain1D[cnt++] = A[3];
+// Ain1D[cnt++] = A[4];
+// Ain1D[cnt++] = A[2];
+// Ain1D[cnt++] = A[4];
+// Ain1D[cnt++] = A[5];
+//
+// int i;
+// std::cout << std::endl;
+// for (i = 0; i < 9; i++) {
+// std::cout << Ain1D[i] << " ";
+// if ((i == 1) || (i == 4))
+// std::cout << std::endl;
+// }
+// std::cout << std::endl << std::endl;
+//
+// mat = cvMat(3, 3, CV_32FC1, Ain1D);
+//
+// cvEigenVV(&mat, eigenVec, eigenVal, 1e-300);
+//
+// for(int i=0;i < 3;i++){
+// eigenvalues[i] = cvmGet(eigenVal,i,0); //Fetching Eigen Value
+//
+// for(int j=0;j < 3;j++){
+// eigenvectors[i][j] = cvmGet(eigenVec,i,j); //Fetching each component of Eigenvector i
+// }
+// }
+//
+// return;
+
+/*
+ const T& A11 = A[0];
+ const T& A12 = A[1];
+ const T& A13 = A[2];
+ const T& A22 = A[3];
+ const T& A23 = A[4];
+ const T& A33 = A[5];
+
+ T& eig1 = eigenvalues[0];
+ T& eig2 = eigenvalues[1];
+ T& eig3 = eigenvalues[2];
+
+ T p1 = A12 * A12 + A13 * A13 + A23 * A23;
+ if( p1 < SmallEpsilon ) { // if A is diagonal
+ eig1 = A11;
+ eig2 = A22;
+ eig3 = A33;
+ memset( eigenvectors, 0, sizeof(T) * 9 );
+ eigenvectors[0][0] = eigenvectors[1][1] = eigenvectors[2][2] = 1;
+ }
+ else { // if A is not diagonal
+ // Compute 1/3 of the trace of matrix A: trace(A)/3
+ T q = ( A11 + A22 + A33 ) / 3;
+ T p2 = (A11-q)*(A11-q) + (A22-q)*(A22-q) + (A33-q)*(A33-q) + 2 * p1;
+ T p = sqrt( p2 / 6);
+
+ // Construct matrix B
+ // B = (1 / p) * (A - q * I), where I is the identity matrix
+ T B11 = (1 / p) * (A11-q);
+ T B12 = (1 / p) * A12;
+ T& B21 = B12;
+ T B13 = (1 / p) * A13;
+ T& B31 = B13;
+ T B22 = (1 / p) * (A22-q);
+ T B23 = (1 / p) * A23;
+ T& B32 = B23;
+ T B33 = (1 / p) * (A33-q);
+
+ // Determinant of a 3 by 3 matrix B
+ // Reference: http://www.mathworks.com/help/aeroblks/determinantof3x3matrix.html
+ T detB = B11*(B22*B33-B23*B32) - B12*(B21*B33-B23*B31) + B13*(B21*B32-B22*B31);
+
+ // In exact arithmetic for a symmetric matrix -1 <= r <= 1
+ // but computation error can leave it slightly outside this range.
+ T r = detB / 2;
+ T phi;
+ const T M_PI3 = T(3.14159265 / 3);
+ if( r <= -1.0f ) {
+ phi = M_PI3;
+ }
+ else if (r >= 1.0f)
+ phi = 0;
+ else {
+ phi = acos(r) / 3;
+ }
+
+ // The eigenvalues satisfy eig3 <= eig2 <= eig1
+ // Notice that: trace(A) = eig1 + eig2 + eig3
+ eig1 = q + 2 * p * cos( phi );
+ eig3 = q + 2 * p * cos( phi + 2 * M_PI3 );
+ eig2 = 3 * q - eig1 - eig3;
+
+ // make sure that |eig1| >= |eig2| >= |eig3|
+ if( abs(eig1) < abs(eig2) ) std::swap(eig1, eig2);
+ if( abs(eig2) < abs(eig3) ) std::swap(eig2, eig3);
+ if( abs(eig1) < abs(eig2) ) std::swap(eig1, eig2);
+
+ // If eig1 is too small, it is not worth to compute the eigenvectors
+ if( abs(eig1) < BigEpsilon ) {
+ memset( eigenvectors, 0, sizeof(T) * 9 );
+ eigenvectors[0][0] = eigenvectors[1][1] = eigenvectors[2][2] = 1;
+ }
+
+ // Compute the corresponding eigenvectors
+ // Reference: [Cayley-Hamilton_theorem](http://en.wikipedia.org/wiki/Cayley-Hamilton_theorem)
+ // If eig1, eig2, eig3 are the distinct eigenvalues of the matrix A;
+ // that is: eig1 != eig2 != eig3. Then
+ // (A - eig1 * I)(A - eig2 * I)(A - eig3 * I) = 0
+ // Thus the columns of the product of any two of these matrices
+ // will contain an eigenvector for the third eigenvalue.
+ else if( abs(eig1-eig2)<BigEpsilon && abs(eig2-eig3)<BigEpsilon ) {
+ memset( eigenvectors, 0, sizeof(T) * 9 );
+ eigenvectors[0][0] = eigenvectors[1][1] = eigenvectors[2][2] = 1;
+ }
+ else if ( abs(eig1-eig2)<BigEpsilon ) {
+ // tested
+ eigenvector_from_value( A, eig1, eig2, eigenvectors[2] );
+ normal_vectors( eigenvectors[2], eigenvectors[1], eigenvectors[0] );
+ }
+ else if ( abs(eig2-eig3)<BigEpsilon ) {
+ // tested
+ eigenvector_from_value( A, eig2, eig3, eigenvectors[0] );
+ normal_vectors( eigenvectors[0], eigenvectors[1], eigenvectors[2] );
+ }
+ else {
+ // eig1!=eig2 && eig2!=eig3 && eig1!=eig3
+ // tested
+ eigenvector_from_value( A, eig2, eig3, eigenvectors[0] );
+ eigenvector_from_value( A, eig1, eig3, eigenvectors[1] );
+ eigenvector_from_value( A, eig1, eig2, eigenvectors[2] );
+ }
+ }
+*/
+}
+
+/*
+template<typename T=float>
+inline void eigenvector_from_value( const T A[6],
+ const T& eig1, const T& eig2, T eigenvector[3] )
+{
+ const T& A11 = A[0];
+ const T& A12 = A[1];
+ const T& A13 = A[2];
+ const T& A22 = A[3];
+ const T& A23 = A[4];
+ const T& A33 = A[5];
+
+ T tempLength = 0, length = 0;
+ T tempVector[3];
+ T colum[3];
+ for( int i=0; i<3; i++ ) {
+ switch( i ) {
+ case 0:
+ colum[0] = A11 - eig2;
+ colum[1] = A12;
+ colum[2] = A13;
+ break;
+ case 1:
+ colum[0] = A12;
+ colum[1] = A22 - eig2;
+ colum[2] = A23;
+ break;
+ default:
+ colum[0] = A13;
+ colum[1] = A23;
+ colum[2] = A33 - eig2;
+ break;
+ }
+ tempVector[0] = colum[0] * (A11 - eig1) + colum[1] * A12 + colum[2] * A13;
+ tempVector[1] = colum[0] * A12 + colum[1] * (A22 - eig1) + colum[2] * A23;
+ tempVector[2] = colum[0] * A13 + colum[1] * A23 + colum[2] * (A33 - eig1);
+
+ tempLength = 0;
+ for( int k=0; k<3; k++ ) {
+ tempLength += tempVector[k] * tempVector[k];
+ }
+
+ // There are three columns in each resulting matrix of ( A-lambda_i * I ) * ( A-lambda_j * I )
+ // We choose the one with the maximum length for the sake of accuracy.
+ if( length < tempLength ) {
+ length = tempLength;
+ memcpy( eigenvector, tempVector, sizeof(T)*3);
+ }
+
+ if( length > BigEpsilon ) break;
+ }
+
+ // The vector is almost zero,
+// assert( length >= SmallEpsilon && "All eigenvector are zero. " );
+
+ // TODO: for debuging
+// if( length < SmallEpsilon ) {
+// std::cerr << "Lenght of vector is too small. Maybe problematic. " << std::endl;
+// }
+
+ // nromailized the vector
+ length = sqrt( length );
+ for( int k=0; k<3; k++ ) eigenvector[k] /= length;
+}
+*/
diff --git a/brick/processing/frangi_c/src/eigen.h b/brick/processing/frangi_c/src/eigen.h
new file mode 100644
index 0000000000000000000000000000000000000000..1585f407dc309aa6c8cca77a482207f6911d6cd3
--- /dev/null
+++ b/brick/processing/frangi_c/src/eigen.h
@@ -0,0 +1,37 @@
+/* Compute the eigen value decomposition of a 3 by 3 symetric matrix
+Input:
+A = [ A11 A12 A13 = [ A[0] A[1] A[2]
+A21 A22 A23 0 A[3] A[4]
+A31 A32 A33 ]; 0 0 A[5] ];
+Output:
+eigenvalues, eigenvectors
+
+Reference: http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices
+*/
+// #include <cmath>
+// #include <assert.h>
+// #include <algorithm>
+// #include <iostream>
+
+
+
+
+// eigenvalue decomposition of a three by three symmetric matrix
+template<typename T=float>
+void eigen_decomp( const T A[6], T eigenvalues[3], T eigenvectors[3][3] );
+
+// give a 3 by 3 symmetric matrix and its eigenvalues, compute the corresponding eigenvectors
+template<typename T=float>
+inline void eigenvector_from_value( const T A[6], const T& eig1, const T& eig2, T eigenvector[3] );
+
+// normalize a vector
+template<typename T=float>
+inline void normalize( T v[3] );
+
+// compute the cross product of two vectors
+template<typename T=float>
+inline void cross_product( const T v1[3], const T v2[3], T v3[3] );
+
+// given a vector v1, compute two arbitrary normal vectors v2 and v3
+template<typename T=float>
+void normal_vectors( const T v1[3], T v2[3], T v3[3] ) ;
diff --git a/brick/processing/frangi_c/src/kernel3D.cpp b/brick/processing/frangi_c/src/kernel3D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..61c440fcabafeabf08b82051b2cc19be6d9ee73f
--- /dev/null
+++ b/brick/processing/frangi_c/src/kernel3D.cpp
@@ -0,0 +1,96 @@
+#include "kernel3D.h"
+#include <opencv2/imgproc/imgproc.hpp>
+
+
+template<typename T>
+Kernel3D<T>::Kernel3D( const cv::Vec3i& n_size )
+{
+ reset( n_size );
+}
+
+
+template<typename T>
+void Kernel3D<T>::reset( const cv::Vec3i& n_size )
+{
+ Data3D<T>::reset( n_size, 0 );
+
+ for( int i=0; i<3; i++ )
+ {
+ center[i] = n_size[i]/2;
+ min_p[i] = -center[i];
+ max_p[i] = ( n_size[i]%2!=0 ) ? ( center[i]+1 ) : ( center[i] );
+ }
+}
+
+
+template <typename U>
+std::ostream& operator<<(std::ostream& out, const Kernel3D<U>& data)
+{
+ // diable output if the size is too big
+ for( int i=0; i<3; i++ )
+ {
+ if ( data.get_size(i)>9 )
+ {
+ std::cout << "I am so sorry. data size is too big to display." << std::endl;
+ return out;
+ }
+ }
+ // output data
+ int x, y, z;
+ for ( z=0; z<data.get_size(2); z++ )
+ {
+ out << "Level " << z << std::endl;
+ for ( y=0; y<data.get_size(1); y++ )
+ {
+ std::cout << "\t";
+ for( x=0; x<data.get_size(0); x++ )
+ {
+ out.precision(3);
+ out << std::scientific << data.at( x, y, z ) << " ";
+ }
+ out << std::endl;
+ }
+ }
+ return out;
+}
+
+
+
+template<typename T>
+Kernel3D<T> Kernel3D<T>::GaussianFilter3D( cv::Vec3i size )
+{
+ smart_assert( typeid(T)==typeid(float) || typeid(T)==typeid(double),
+ "Can only use float or double" );
+
+ for( int i=0; i<3; i++ )
+ smart_assert( size[i]>0 && size[i]%2!=0, "size should be possitive odd number" );
+
+ // Calculate sigma from size:
+ // sigma = 0.3 * ( (ksize-1)/2 - 1 ) + 0.8 = 0.15*ksize + 0.35
+ // reference: OpenCv Documentation 2.6.4.0 getGaussianKernel
+ // http://docs.opencv.org/modules/imgproc/doc/filtering.html#creategaussianfilter
+ // If we calculate sigma based on 99.7% Rule:
+ // sigma = ( size - 1 ) /6 = 0.17 * size - 0.17
+ // But we are flowing the definition of OpenCv here
+ cv::Mat gaussian[3];
+ for( int i=0; i<3; i++ ) gaussian[i] = cv::getGaussianKernel( size[i], 0, CV_64F );
+
+ Kernel3D<T> kernel(size);
+ int x, y, z;
+ for( z=0; z<kernel.get_size_z(); z++ )
+ {
+ for( y=0; y<kernel.get_size_y(); y++ )
+ {
+ for( x=0; x<kernel.get_size_x(); x++ )
+ {
+ kernel.at(x,y,z)
+ = gaussian[0].at<double>(x)
+ * gaussian[1].at<double>(y)
+ * gaussian[2].at<double>(z);
+ }
+ }
+ }
+
+ return kernel;
+}
+
diff --git a/brick/processing/frangi_c/src/kernel3D.h b/brick/processing/frangi_c/src/kernel3D.h
new file mode 100644
index 0000000000000000000000000000000000000000..be82fb603c43bcb807c430ed7e8c6adbc7b95f5b
--- /dev/null
+++ b/brick/processing/frangi_c/src/kernel3D.h
@@ -0,0 +1,78 @@
+#include "data3D.h"
+#include <opencv2/imgproc/imgproc.hpp>
+
+// template<typename T> class Data3D;
+
+template<typename T = float>
+class Kernel3D : public Data3D<T>
+{
+ // operator overload
+ template <typename U>
+ friend std::ostream& operator<<(std::ostream& out, const Kernel3D<U>& data);
+
+public:
+ // constructor and destructors
+ Kernel3D() {};
+ Kernel3D( const cv::Vec3i& n_size );
+ virtual ~Kernel3D() {}
+ // getters
+ const int& min_pos(int i) const
+ {
+ return min_p[i];
+ }
+ const int& max_pos(int i) const
+ {
+ return max_p[i];
+ }
+ virtual const inline T& offset_at(const int& x, const int& y, const int& z) const
+ {
+ return Data3D<T>::at(x+center[0], y+center[1], z+center[2]);
+ }
+ virtual inline T& offset_at(const int& x, const int& y, const int& z)
+ {
+ return Data3D<T>::at(x+center[0], y+center[1], z+center[2]);
+ }
+
+ // setters
+ void reset( const cv::Vec3i& n_size );
+private:
+ // minimun and maximum possition
+ cv::Vec3i min_p, max_p, center;
+public:
+ // Some Global Functions for Getting 3D Kernel's that are commonly used.
+ static Kernel3D<T> GaussianFilter3D( cv::Vec3i size );
+ inline static Kernel3D<T> GaussianFilter3D( int ksize )
+ {
+ smart_assert( typeid(T)==typeid(float) || typeid(T)==typeid(double),
+ "Can only use float or double" );
+ return GaussianFilter3D( cv::Vec3i(ksize, ksize, ksize) );
+ }
+ inline static Kernel3D<T> dx()
+ {
+ smart_assert( typeid(T)==typeid(float) || typeid(T)==typeid(double),
+ "Can only use float or double" );
+ Kernel3D<T> dx( cv::Vec3i(3,1,1) );
+ dx.at( 0, 0, 0 ) = -0.5;
+ dx.at( 2, 0, 0 ) = 0.5;
+ return dx;
+ }
+ inline static Kernel3D<T> dy()
+ {
+ smart_assert( typeid(T)==typeid(float) || typeid(T)==typeid(double),
+ "Can only use float or double" );
+ Kernel3D<T> dy( cv::Vec3i(1,3,1) );
+ dy.at( 0, 0, 0 ) = -0.5;
+ dy.at( 0, 2, 0 ) = 0.5;
+ return dy;
+ }
+ inline static Kernel3D<T> dz()
+ {
+ smart_assert( typeid(T)==typeid(float) || typeid(T)==typeid(double),
+ "Can only use float or double" );
+ Kernel3D<T> dz( cv::Vec3i(1,1,3) );
+ dz.at( 0, 0, 0 ) = -0.5;
+ dz.at( 0, 0, 2 ) = 0.5;
+ return dz;
+ }
+
+};
diff --git a/brick/processing/frangi_c/src/main.cpp b/brick/processing/frangi_c/src/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..21584090561e5bd5894df9713389b182bea23283
--- /dev/null
+++ b/brick/processing/frangi_c/src/main.cpp
@@ -0,0 +1,64 @@
+#include <iostream>
+#include "detector.h"
+#include "data3D.h"
+#include "processing.h"
+
+using namespace std;
+using namespace cv;
+
+extern "C" {
+void vesselness(
+ float* array, const int width, const int height, const int depth,
+ float* response,
+ float sigma_from = 1.0f,
+ float sigma_to = 8.10f,
+ float sigma_step = 0.3f,
+ float alpha = 1.0e-1f,
+ float beta = 5.0e0f,
+ float gamma = 3.5e5f )
+{
+ // Sigma: Parameters for Vesselness
+ // [sigma_from, sigma_to]: the potential size rang of the vessels
+ // sigma_step: precision of computation
+ // Parameters for vesselness, please refer to Frangi's paper
+ // or this [blog](http://yzhong.co/?p=351)
+
+ Data3D<float> im_short(width, height, depth);
+ im_short.SetFromArray(array);
+
+ Data3D<Vesselness_Sig> vn_sig;
+ VesselDetector::compute_vesselness( im_short, vn_sig, sigma_from, sigma_to, sigma_step,
+ alpha, beta, gamma );
+
+// Data3D<Vesselness_Sig> vn_sig_nms;
+// IP::non_max_suppress( vn_sig, vn_sig_nms );
+//
+// Data3D<Vesselness_Sig> vn_sig_et;
+// IP::edge_tracing( vn_sig_nms, vn_sig_et, 0.45f, 0.05f );
+
+ int x,y,z;
+ float* array_cursor = response;
+ for( z=0; z<vn_sig.get_size_z(); z++ ) {
+ for( y=0; y<vn_sig.get_size_y(); y++ ) {
+ for( x=0; x<vn_sig.get_size_x(); x++ ) {
+ *array_cursor = vn_sig.at(x,y,z).rsp;
+ array_cursor++;
+ }
+ }
+ }
+// vn_sig_et.CopyToArray(array);
+// for (int idx=0; idx < 20; idx++) {
+// cout << array[idx] << endl;
+// }
+}
+}
+
+int main(void)
+{
+ float* array = new float[100*100*100];
+ float* response = new float[100*100*100];
+
+ vesselness( array, 100,100,100, response );
+
+ return 0;
+}
diff --git a/brick/processing/frangi_c/src/processing.cpp b/brick/processing/frangi_c/src/processing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..984d645050895c2a97a4e70119a6e88cd7981f51
--- /dev/null
+++ b/brick/processing/frangi_c/src/processing.cpp
@@ -0,0 +1,156 @@
+#include <opencv2/core/core.hpp>
+#include <queue>
+#include <iostream>
+#include "types.h"
+
+#include "processing.h"
+#include "processing_code.h"
+
+using namespace cv;
+using namespace std;
+
+void ImageProcessing::non_max_suppress( const Data3D<Vesselness_Sig>& src, Data3D<Vesselness_Sig>& dst )
+{
+ dst.reset( src.get_size() );
+
+ static const float s2 = sqrt(1/2.0f);
+ static const float s3 = sqrt(1/3.0f);
+
+ // 13 major orientations in 3d
+ static const int MAJOR_DIR_NUM = 13;
+ static const Vec3f major_dirs[MAJOR_DIR_NUM] = {
+ // directions along the axis, there are 3 of them
+ Vec3f(1, 0, 0), Vec3f(0, 1, 0), Vec3f(0, 0, 1),
+ // directions that lie within the 3D plane of two axis, there are 6 of them
+ Vec3f(s2, s2,0), Vec3f(0,s2, s2), Vec3f(s2,0, s2),
+ Vec3f(s2,-s2,0), Vec3f(0,s2,-s2), Vec3f(s2,0,-s2),
+ // directions that are equally in between three axis, there are 4 of them
+ Vec3f(s3,s3,s3), Vec3f(s3,s3,-s3), Vec3f(s3,-s3,s3), Vec3f(s3,-s3,-s3)
+ };
+
+ // there are 26 directions in 3D
+ static const int NUM_DIR_3D = 26;
+ Vec3i neighbor3d[NUM_DIR_3D] = {
+ // directions along the axis, there are 6 of them
+ Vec3i(0,0, 1), Vec3i(0, 1,0), Vec3i( 1,0,0),
+ Vec3i(0,0,-1), Vec3i(0,-1,0), Vec3i(-1,0,0),
+ // directions that lie within the 3D plane of two axis, there are 12 of them
+ Vec3i(0,1,1), Vec3i(0, 1,-1), Vec3i( 0,-1,1), Vec3i( 0,-1,-1),
+ Vec3i(1,0,1), Vec3i(1, 0,-1), Vec3i(-1, 0,1), Vec3i(-1, 0,-1),
+ Vec3i(1,1,0), Vec3i(1,-1, 0), Vec3i(-1, 1,0), Vec3i(-1,-1, 0),
+ // directions that are equally in between three axis, there are 8 of them
+ Vec3i( 1,1,1), Vec3i( 1,1,-1), Vec3i( 1,-1,1), Vec3i( 1,-1,-1),
+ Vec3i(-1,1,1), Vec3i(-1,1,-1), Vec3i(-1,-1,1), Vec3i(-1,-1,-1),
+ };
+
+ // cross section for the major orientation
+ vector<Vec3i> cross_section[MAJOR_DIR_NUM];
+
+ // Setting offsets that are perpendicular to dirs
+ for( int i=0; i<MAJOR_DIR_NUM; i++ ) {
+ for( int j=0; j<NUM_DIR_3D; j++ ) {
+ // multiply the two directions
+ float temp = major_dirs[i].dot( neighbor3d[j] );
+ // the temp is 0, store this direciton
+ if( abs(temp)<1.0e-5 ) {
+ cross_section[i].push_back( neighbor3d[j] );
+ }
+ }
+ }
+
+ int x,y,z;
+ for( z=0; z<src.get_size_z(); z++ ) {
+ for( y=0; y<src.get_size_y(); y++ ) {
+ for( x=0; x<src.get_size_x(); x++ ) {
+ // non-maximum surpression
+ bool isMaximum = true;
+
+ // Method I:
+ // find the major orientation
+ // assigning the orientation to one of the 13 categories
+ const Vec3f& cur_dir = src.at(x,y,z).dir;
+ int mdi = 0; // major direction id
+ float max_dot_product = 0;
+ for( int di=0; di<MAJOR_DIR_NUM; di++ ) {
+ float current_dot_product = abs( cur_dir.dot(major_dirs[di]) );
+ if( max_dot_product < current_dot_product ) {
+ max_dot_product = current_dot_product;
+ mdi = di;// update the major direction id
+ }
+ }
+ for( unsigned int i=0; i<cross_section[mdi].size(); i++ ) {
+ int ox = x + cross_section[mdi][i][0];
+ int oy = y + cross_section[mdi][i][1];
+ int oz = z + cross_section[mdi][i][2];
+ if( src.isValid(ox,oy,oz) && src.at(x,y,z).rsp < src.at(ox,oy,oz).rsp ) {
+ isMaximum = false;
+ break;
+ }
+ }
+
+ // Method II: Based on Sigma
+ //int sigma = (int) ceil( src.at(x,y,z).sigma );
+ //for( int i = -sigma; i <= sigma; i++ ) for( int j = -sigma; j <= sigma; j++ ) {
+ // Vec3i offset = src.at(x,y,z).normals[0] * i + src.at(x,y,z).normals[1] * j;
+ // int ox = x + offset[0];
+ // int oy = y + offset[1];
+ // int oz = z + offset[2];
+ // if( src.isValid(ox,oy,oz) && src.at(x,y,z).rsp < src.at(ox,oy,oz).rsp ) {
+ // isMaximum = false; break;
+ // }
+ //}
+
+ if( isMaximum ) {
+ dst.at(x,y,z) = src.at(x,y,z);
+ }
+ }
+ }
+ }
+}
+
+
+void ImageProcessing::edge_tracing( Data3D<Vesselness_Sig>& src, Data3D<Vesselness_Sig>& dst, const float& thres1, const float& thres2 )
+{
+ Data3D<float> src1d;
+ src.copyDimTo( src1d, 0 );
+ IP::normalize( src1d, 1.0f );
+
+ int x, y, z;
+ std::queue<Vec3i> q;
+ Data3D<unsigned char> mask( src.get_size() );
+ for(z=0; z<src.SZ(); z++) for (y=0; y<src.SY(); y++) for(x=0; x<src.SX(); x++) {
+ if( src1d.at(x,y,z) > thres1 ) {
+ q.push( Vec3i(x,y,z) );
+ mask.at(x,y,z) = 255;
+ }
+ }
+
+ Vec3i dir[26];
+ for( int i=0; i<26; i++ ) {
+ int index = (i + 14) % 27;
+ dir[i][0] = index/9%3 - 1;
+ dir[i][1] = index/3%3 - 1;
+ dir[i][2] = index/1%3 - 1;
+ }
+
+ while( !q.empty() ) {
+ Vec3i pos = q.front();
+ q.pop();
+ Vec3i off_pos;
+ for( int i=0; i<26; i++ ) {
+ off_pos = pos + dir[i];
+ if( src.isValid(off_pos) && !mask.at( off_pos ) && src1d.at(off_pos) > thres2 ) {
+ mask.at( off_pos ) = 255;
+ q.push( off_pos );
+ }
+ }
+ }
+
+ dst.reset( src.get_size() );
+ for(z=0; z<src.SZ(); z++) for (y=0; y<src.SY(); y++) for(x=0; x<src.SX(); x++) {
+ if( mask.at( x,y,z ) ) {
+ dst.at( x,y,z ) = src.at( x,y,z );
+ // dst.at( x,y,z ).rsp = sqrt( src1d.at( x,y,z ) );
+ }
+ }
+}
diff --git a/brick/processing/frangi_c/src/processing.h b/brick/processing/frangi_c/src/processing.h
new file mode 100644
index 0000000000000000000000000000000000000000..ae736119cc88d2dd7943759c5958e880195c0b78
--- /dev/null
+++ b/brick/processing/frangi_c/src/processing.h
@@ -0,0 +1,14 @@
+#include "data3D.h"
+
+
+namespace ImageProcessing {
+
+void non_max_suppress( const Data3D<Vesselness_Sig>& src, Data3D<Vesselness_Sig>& dst );
+
+void edge_tracing( Data3D<Vesselness_Sig>& src, Data3D<Vesselness_Sig>& dst, const float& thres1, const float& thres2 );
+
+void dir_tracing( Data3D<Vesselness_All>& src_vn, Data3D<Vesselness_Sig>& dst );
+
+void edge_tracing_mst( Data3D<Vesselness_All>& src_vn, Data3D<Vesselness_Sig>& dst, const float& thres1, const float& thres2 );
+
+}
diff --git a/brick/processing/frangi_c/src/processing_code.h b/brick/processing/frangi_c/src/processing_code.h
new file mode 100644
index 0000000000000000000000000000000000000000..6221073d2ed8a22a021de99afbc58688af4ca803
--- /dev/null
+++ b/brick/processing/frangi_c/src/processing_code.h
@@ -0,0 +1,121 @@
+#pragma once
+
+#include <vector>
+#include "data3D.h"
+#include "kernel3D.h"
+#include "smart_assert.h"
+
+
+class Vesselness;
+class Vesselness_Sig;
+class Vesselness_Nor;
+class Vesselness_All;
+
+namespace ImageProcessing {
+template<typename T1, typename T2>
+bool GaussianBlur3D( const Data3D<T1>& src, Data3D<T2>& dst, int ksize, double sigma = 0.0 );
+
+template<typename T>
+void normalize( Data3D<T>& data, T norm_max = 1);
+
+
+///////////////////////////////////////////////////////////////////////////
+// NON-MAXIMUM SUPPRESSION
+void non_max_suppress( const Data3D<Vesselness_Sig>& src, Data3D<Vesselness_Sig>& dst );
+void edge_tracing( Data3D<Vesselness_Sig>& src, Data3D<Vesselness_Sig>& dst, const float& thres1 = 0.85f, const float& thres2 = 0.10f );
+void dir_tracing( Data3D<Vesselness_All>& src, Data3D<Vesselness_Sig>& res_dt );
+void edge_tracing_mst( Data3D<Vesselness_All>& src, Data3D<Vesselness_Sig>& dst, const float& thres1 = 0.85f, const float& thres2 = 0.10f );
+}
+
+namespace IP = ImageProcessing;
+
+
+
+template<typename T1, typename T2>
+bool ImageProcessing::GaussianBlur3D( const Data3D<T1>& src, Data3D<T2>& dst, int ksize, double sigma )
+{
+ smart_return( ksize%2!=0, "kernel size should be odd number", false );
+
+ //////////////////////////////////////////////////////////////////////////////////////
+ // Relationship between Sigma and Kernal Size (ksize)
+ /////////////////////////////
+ // Option I - Based on OpenCV
+ // Calculate sigma from size:
+ // sigma = 0.3 * ( (ksize-1)/2 - 1 ) + 0.8 = 0.15*ksize + 0.35
+ // Calculate size from sigma:
+ // ksize = ( sigma - 0.35 ) / 0.15 = 6.67 * sigma - 2.33
+ // Reference: OpenCv Documentation 2.6.4.0 getGaussianKernel
+ // http://docs.opencv.org/modules/imgproc/doc/filtering.html#creategaussianfilter
+ // Option II - Based on the traditional 99.7%
+ // Calculate size from sigma:
+ // ksize = 6 * sigma + 1
+ // Calculate sigma from ksize:
+ // sigma = (size - 1)/6 = 0.17 * size - 0.17
+
+ cv::Mat gk = cv::getGaussianKernel( ksize, sigma, CV_64F );
+
+ int hsize = ksize/2;
+ smart_assert( 2*hsize+1==ksize, "Alert: Bug!" );
+
+ cv::Vec3i spos, epos;
+ spos = cv::Vec3i(0, 0, 0);
+ epos = src.get_size();
+
+ // gaussian on x-direction
+ Data3D<T2> tmp1( src.get_size(), T2(0) ); // the data will be set to zero
+ #pragma omp parallel for
+ for( int z=spos[2]; z<epos[2]; z++ ) for( int y=spos[1]; y<epos[1]; y++ ) for( int x=spos[0]; x<epos[0]; x++ ) {
+ double sum = 0.0;
+ for( int i=0; i<ksize; i++) {
+ int x2 = x+i-hsize;
+ if( x2>=0 && x2<epos[0] ) {
+ tmp1.at(x, y, z) = T2( gk.at<double>(i) * src.at(x2, y, z) + tmp1.at(x, y, z));
+ sum += gk.at<double>(i);
+ }
+ }
+ tmp1.at(x, y, z) = T2( tmp1.at(x, y, z)/sum );
+ }
+
+
+ // gaussian on y-direction
+ Data3D<T2> tmp2( src.get_size(), T2(0) ); // the data will be set to zero
+ #pragma omp parallel for
+ for( int z=spos[2]; z<epos[2]; z++ ) for( int y=spos[1]; y<epos[1]; y++ ) for( int x=spos[0]; x<epos[0]; x++ ) {
+ double sum = 0.0;
+ for( int i=0; i<ksize; i++) {
+ int y2 = y+i-hsize;
+ if( y2>=0 && y2<epos[1] ) {
+ tmp2.at(x, y, z) = T2( gk.at<double>(i) * tmp1.at(x, y2, z) + tmp2.at(x, y, z));
+ sum += gk.at<double>(i);
+ }
+ }
+ tmp2.at(x, y, z) = T2( tmp2.at(x, y, z)/sum );
+ }
+
+ tmp1.reset(); // tmp1 is no long in use. release memory
+
+ // gaussian on z-direction
+ dst.reset( src.get_size() );
+ #pragma omp parallel for
+ for( int z=spos[2]; z<epos[2]; z++ ) for( int y=spos[1]; y<epos[1]; y++ ) for( int x=spos[0]; x<epos[0]; x++ ) {
+ double sum = 0.0;
+ for( int i=0; i<ksize; i++) {
+ int z2 = z+i-hsize;
+ if( z2>=0 && z2<epos[2] ) {
+ dst.at(x, y, z) = T2( dst.at(x, y, z) + gk.at<double>(i) * tmp2.at(x, y, z2) );
+ sum += gk.at<double>(i);
+ }
+ }
+ dst.at(x, y, z) = T2( dst.at(x, y, z)/sum );
+ }
+
+ return true;
+}
+
+
+template<typename T>
+void ImageProcessing::normalize( Data3D<T>& data, T norm_max )
+{
+ cv::Vec<T, 2> min_max = data.get_min_max_value();
+ data.getMat() = 1.0 * norm_max / (min_max[1]-min_max[0]) * (data.getMat() - min_max[0]);
+}
diff --git a/brick/processing/frangi_c/src/smart_assert.h b/brick/processing/frangi_c/src/smart_assert.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6d38fddb212e1cf8bf02fcdde4f048a8044c683
--- /dev/null
+++ b/brick/processing/frangi_c/src/smart_assert.h
@@ -0,0 +1,41 @@
+/* Smart Assertion
+ * - Show assertion under debug mode
+ * - Print assertion message udner release mode
+ *
+ * Created by Yuchen on Apr 30th, 2013
+ * Last Modified by Yuchen on Jul 2nd, 2013
+ *
+ * Example:
+ Code:
+ int i = -1;
+ smart_assert(i>=0 && i<3, "Index out of bound" );
+ Output:
+ Assertion failed: i>0 && i<3
+ Messages: Index out of bound
+ Location: file main.cpp, line 17
+ */
+#ifndef SMART_ASSERT_H
+#define SMART_ASSERT_H
+#include <assert.h>
+#include <iostream>
+
+// The following color is defined using [ANSI colour codes](http://en.wikipedia.org/wiki/ANSI_escape_code)
+#define SMA_RED "\033[0;31m"
+#define SMA_BLACK "\x1b[0;49m"
+
+#define smart_assert( condition, message ) \
+ if( !(condition) ) { \
+ std::cerr << std::endl; \
+ std::cerr << SMA_RED << "Assertion failed: " << (#condition) << std::endl; \
+ std::cerr << " Messages: " << message << std::endl; \
+ std::cerr << " Location: file "<< __FILE__ << ", line " << __LINE__ << std::endl; \
+ std::cerr << SMA_BLACK << std::endl;\
+ }
+
+#define smart_return( condition, message, return_value ) \
+ smart_assert( condition, message )\
+ if( !(condition) ) { \
+ return return_value;\
+ }
+
+#endif
diff --git a/brick/processing/frangi_c/src/type_info.h b/brick/processing/frangi_c/src/type_info.h
new file mode 100644
index 0000000000000000000000000000000000000000..b339e18de228817e6fd7d2ddf27e1705b5537486
--- /dev/null
+++ b/brick/processing/frangi_c/src/type_info.h
@@ -0,0 +1,81 @@
+#pragma once
+
+#include <string>
+#include <opencv2/core/core.hpp>
+
+template <class T> struct TypeInfo {
+ static std::string str()
+ {
+ return "undefined";
+ }
+ static int CV_TYPE()
+ {
+ return -1;
+ }
+};
+
+template <> struct TypeInfo <int> {
+ static std::string str()
+ {
+ return "int";
+ }
+ static int CV_TYPE()
+ {
+ return CV_32S;
+ }
+};
+
+template <> struct TypeInfo <short> {
+ static std::string str()
+ {
+ return "short";
+ }
+ static int CV_TYPE()
+ {
+ return CV_16S;
+ }
+};
+
+template <> struct TypeInfo <float> {
+ static std::string str()
+ {
+ return "float";
+ }
+ static int CV_TYPE()
+ {
+ return CV_32F;
+ }
+};
+
+template <> struct TypeInfo <double> {
+ static std::string str()
+ {
+ return "double";
+ }
+ static int CV_TYPE()
+ {
+ return CV_64F;
+ }
+};
+
+template <> struct TypeInfo <unsigned short> {
+ static std::string str()
+ {
+ return "ushort";
+ }
+ static int CV_TYPE()
+ {
+ return CV_16U;
+ }
+};
+
+template <> struct TypeInfo <unsigned char> {
+ static std::string str()
+ {
+ return "uchar";
+ }
+ static int CV_TYPE()
+ {
+ return CV_8U;
+ }
+};
diff --git a/brick/processing/frangi_c/src/types.h b/brick/processing/frangi_c/src/types.h
new file mode 100644
index 0000000000000000000000000000000000000000..3b647840ee5491e422b164b835965e8328196ac6
--- /dev/null
+++ b/brick/processing/frangi_c/src/types.h
@@ -0,0 +1,297 @@
+#pragma once
+
+#include <opencv2/core/core.hpp>
+#include "type_info.h"
+#include "smart_assert.h"
+
+// VesselnessTypes.h
+// Define 4 different vesselness data types
+// 1) class Vesselness
+// vesselness response (rsp) + vessel orientation (dir)
+// 2) class Vesselness_Sig
+// vesselness response (rsp) + vessel orientation (dir) + vessel size (sigma)
+// 3) class Vesselness_Nor
+// vesselness response (rsp) + vessel orientation (dir) + vessel orientations normals (normals)
+// 4) class Vesselness_All
+// vesselness response (rsp) + vessel orientation (dir) + vessel orientations normals (normals)
+
+struct Vesselness_Data;
+struct Vesselness_Sig_Data;
+struct Vesselness_Nor_Data;
+struct Vesselness_All_Data;
+
+class Vesselness;
+class Vesselness_Sig;
+class Vesselness_Nor;
+class Vesselness_All;
+
+struct Vesselness_Data {
+ float rsp; // vesselness response
+ cv::Vec3f dir; // vessel dirction
+
+ Vesselness_Data( float s = 0.0f ) : rsp( s ), dir( cv::Vec3f(1,0,0) ) { }
+ bool operator>( const Vesselness_Data& right ) const
+ {
+ return ( (this->rsp) > right.rsp );
+ }
+ bool operator<( const Vesselness_Data& right ) const
+ {
+ return ( (this->rsp) < right.rsp );
+ }
+};
+
+struct Vesselness_Sig_Data : public Vesselness_Data {
+ float sigma; // relative size of the vessel
+
+ Vesselness_Sig_Data( float s = 0.0f ) : Vesselness_Data ( s ) { }
+};
+
+struct Vesselness_Nor_Data : public Vesselness_Data {
+ cv::Vec3f normals[2]; // normals of vessel orientation
+};
+
+struct Vesselness_All_Data : public Vesselness_Data {
+ float sigma; // relative size of the vessel
+ cv::Vec3f normals[2]; // normals of vessel orientation
+};
+
+
+
+
+
+class Vesselness : public Vesselness_Data {
+public:
+ Vesselness( const float& val = 0 ) : Vesselness_Data( val ) { }
+
+ // Size tells us how much number of float we are using in the structure.
+ // The returned value should be:
+ // 4 for Vesselness;
+ // 5 for Vesselness_Sig;
+ // 10 for Vesselness_Nor;
+ // 11 for Vesselness_All;
+ // This virtual function here will take 2*4 = 8 bytes of data
+ static const int _size = sizeof(Vesselness_Data)/sizeof(float);
+ int size(void) const
+ {
+ return _size;
+ }
+
+ const float& operator[]( const int& i ) const
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+
+ float& operator[]( const int& i )
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+};
+
+
+class Vesselness_Nor : public Vesselness_Nor_Data {
+public:
+ static const int _size = sizeof(Vesselness_Nor_Data)/sizeof(float);
+
+ const float& operator[]( const int& i ) const
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+
+ float& operator[]( const int& i )
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+};
+
+class Vesselness_All : public Vesselness_All_Data {
+public:
+ static const int _size = sizeof(Vesselness_All_Data)/sizeof(float);
+
+ Vesselness_All() {}
+
+ Vesselness_All( const Vesselness_Nor& v_Nor, const float& s )
+ {
+ this->dir = v_Nor.dir;
+ this->normals[0] = v_Nor.normals[0];
+ this->normals[1] = v_Nor.normals[1];
+ this->rsp = v_Nor.rsp;
+ this->sigma = s;
+ }
+
+ const float& operator[]( const int& i ) const
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+
+ float& operator[]( const int& i )
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+};
+
+
+class Vesselness_Sig : public Vesselness_Sig_Data {
+ Vesselness_Sig& operator=( const float& s );
+public:
+ static const int _size = sizeof(Vesselness_Sig_Data)/sizeof(float);
+
+ Vesselness_Sig() { }
+ Vesselness_Sig( float s ) : Vesselness_Sig_Data( s ) { }
+ operator Vesselness( )
+ {
+ Vesselness vn;
+ vn.rsp = this->rsp;
+ vn.dir = this->dir;
+ return vn;
+ }
+ operator float()
+ {
+ return this->rsp;
+ }
+
+ Vesselness_Sig( const Vesselness_All& src )
+ {
+ this->rsp = src.rsp;
+ this->dir = src.dir;
+ this->sigma = src.sigma;
+ }
+
+ Vesselness_Sig& operator=( const Vesselness_All& src )
+ {
+ this->rsp = src.rsp;
+ this->dir = src.dir;
+ this->sigma = src.sigma;
+ return (*this);
+ }
+
+ const float& operator[]( const int& i ) const
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+
+ float& operator[]( const int& i )
+ {
+ smart_return( i>=0&&i<_size, "index invalid", *(float*)(this) );
+ return *((float*)(this)+i);
+ }
+};
+
+// The namespace cv here is necessary.
+// Reference: http://stackoverflow.com/questions/2282349/specialization-of-templateclass-tp-struct-stdless-in-different-namespace
+namespace cv {
+template< > class DataType< Vesselness > : public DataType< Vec<float, Vesselness::_size> > { };
+// template< > class DataDepth< Vesselness > : public DataDepth< Vec<float, Vesselness::_size> > { };
+template< > class traits::Type< Vesselness > : public traits::Type< Vec<float, Vesselness::_size> > { };
+
+template< > class DataType< Vesselness_Sig > : public DataType< Vec<float, Vesselness_Sig::_size> > { };
+// template< > class DataDepth< Vesselness_Sig > : public DataDepth< Vec<float, Vesselness_Sig::_size> > { };
+template< > class traits::Type< Vesselness_Sig > : public traits::Type< Vec<float, Vesselness_Sig::_size> > { };
+
+template< > class DataType< Vesselness_Nor > : public DataType< Vec<float, Vesselness_Nor::_size> > { };
+// template< > class DataDepth< Vesselness_Nor > : public DataDepth< Vec<float, Vesselness_Nor::_size> > { };
+template< > class traits::Type< Vesselness_Nor > : public traits::Type< Vec<float, Vesselness_Nor::_size> > { };
+
+template< > class DataType< Vesselness_All > : public DataType< Vec<float, Vesselness_All::_size> > { };
+// template< > class DataDepth< Vesselness_All > : public DataDepth< Vec<float, Vesselness_All::_size> > { };
+template< > class traits::Type< Vesselness_All > : public traits::Type< Vec<float, Vesselness_All::_size> > { };
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////////
+// Reference: These following are defined in <core.hpp>
+//
+//template<> class DataType<int>
+//{
+//public:
+// typedef int value_type;
+// typedef value_type work_type;
+// typedef value_type channel_type;
+// typedef value_type vec_type;
+// enum { generic_type = 0, depth = DataDepth<channel_type>::value, channels = 1,
+// fmt=DataDepth<channel_type>::fmt,
+// type = CV_MAKETYPE(depth, channels) };
+//};
+//
+//template<> class DataType<float>
+//{
+//public:
+// typedef float value_type;
+// typedef value_type work_type;
+// typedef value_type channel_type;
+// typedef value_type vec_type;
+// enum { generic_type = 0, depth = DataDepth<channel_type>::value, channels = 1,
+// fmt=DataDepth<channel_type>::fmt,
+// type = CV_MAKETYPE(depth, channels) };
+//};
+//
+//template<typename _Tp, int cn> class DataType<Vec<_Tp, cn> >
+//{
+//public:
+// typedef Vec<_Tp, cn> value_type;
+// typedef Vec<typename DataType<_Tp>::work_type, cn> work_type;
+// typedef _Tp channel_type;
+// typedef value_type vec_type;
+// enum { generic_type = 0, depth = DataDepth<channel_type>::value, channels = cn,
+// fmt = ((channels-1)<<8) + DataDepth<channel_type>::fmt,
+// type = CV_MAKETYPE(depth, channels) };
+//};
+
+
+template <> struct TypeInfo <Vesselness> {
+ static std::string str()
+ {
+ std::stringstream ss;
+ ss << "float," << Vesselness::_size;
+ return ss.str();
+ }
+ static int CV_TYPE()
+ {
+ return CV_32FC( Vesselness::_size );
+ }
+};
+
+template <> struct TypeInfo <Vesselness_Sig> {
+ static std::string str()
+ {
+ std::stringstream ss;
+ ss << "float," << Vesselness_Sig::_size;
+ return ss.str();
+ }
+ static int CV_TYPE()
+ {
+ return CV_32FC( Vesselness_Sig::_size );
+ }
+};
+
+template <> struct TypeInfo <Vesselness_Nor> {
+ static std::string str()
+ {
+ std::stringstream ss;
+ ss << "float," << Vesselness_Nor::_size;
+ return ss.str();
+ }
+ static int CV_TYPE()
+ {
+ return CV_32FC( Vesselness_Nor::_size );
+ }
+};
+
+template <> struct TypeInfo <Vesselness_All> {
+ static std::string str()
+ {
+ std::stringstream ss;
+ ss << "float," << Vesselness_All::_size;
+ return ss.str();
+ }
+ static int CV_TYPE()
+ {
+ return CV_32FC( Vesselness_All::_size );
+ }
+};
diff --git a/frangi_3d.pxd b/brick/processing/frangi_py/frangi_3d.pxd
similarity index 100%
rename from frangi_3d.pxd
rename to brick/processing/frangi_py/frangi_3d.pxd
diff --git a/frangi_3d.py b/brick/processing/frangi_py/frangi_3d.py
similarity index 96%
rename from frangi_3d.py
rename to brick/processing/frangi_py/frangi_3d.py
index 1f40db5e6ad2e99131feecfbe9160e1930570619..9fbff6bf0aa14fb4f3039566bd2ccfda43856493 100644
--- a/frangi_3d.py
+++ b/brick/processing/frangi_py/frangi_3d.py
@@ -33,7 +33,7 @@
Follows: https://github.com/ellisdg/frangi3d
'''
-from type import array_t
+from brick.general.type import array_t
import multiprocessing as pl_
from ctypes import c_float as c_float_t
diff --git a/map_labeling.py b/brick/processing/map_labeling.py
similarity index 98%
rename from map_labeling.py
rename to brick/processing/map_labeling.py
index ee2dfa8609ceb2a652304e664971d05b9d20896c..ca5ff6f13ddf9407a47cbcea5882d731d316cf36 100644
--- a/map_labeling.py
+++ b/brick/processing/map_labeling.py
@@ -29,7 +29,7 @@
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.
-from type import array_t
+from brick.general.type import array_t
import numpy as np_
diff --git a/plot.py b/brick/processing/plot.py
similarity index 97%
rename from plot.py
rename to brick/processing/plot.py
index 6b4924b99a7f69504c22a587504d06eef5e0c7b3..d9f36493b7f983a0c6d0298d6d1925a26f0c363b 100644
--- a/plot.py
+++ b/brick/processing/plot.py
@@ -29,10 +29,10 @@
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.
-import dijkstra_1_to_n as dk_
-from extension import extension_t
-from soma import soma_t
-from type import array_t, py_array_picker_h
+import brick.processing.dijkstra_1_to_n as dk_
+from brick.component.extension import extension_t
+from brick.component.soma import soma_t
+from brick.general.type import array_t, py_array_picker_h
from typing import Any, Optional, Sequence, Tuple, Union
diff --git a/dijkstra_1_to_n.py b/dijkstra_1_to_n.py
deleted file mode 120000
index c28f560f82fc17f05e899094160a79e5174f378a..0000000000000000000000000000000000000000
--- a/dijkstra_1_to_n.py
+++ /dev/null
@@ -1 +0,0 @@
-../../../../brick/def/dijkstra-img/dijkstra_1_to_n.py
\ No newline at end of file
diff --git a/main.py b/nutrimorph.py
similarity index 97%
rename from main.py
rename to nutrimorph.py
index 674d423ffeaed54ef5b3ce4adb384725ad65eae6..41fdf358980bfce8bee4542bc2c14601054b6710 100644
--- a/main.py
+++ b/nutrimorph.py
@@ -37,11 +37,11 @@
# mprof run main.py
# mprof plot
-import connection as cn_
-import extension as xt_
-import feedback as fb_
+import brick.component.connection as cn_
+import brick.component.extension as xt_
+import brick.component.soma as sm_
+import brick.general.feedback as fb_
# from main_prm import *
-import soma as sm_
import os as os_
import sys as sy_
@@ -73,7 +73,7 @@ max_weighted_length_c = None
min_area_c = None
in_parallel = None
-execfile(sy_.argv[1])
+exec(open(sy_.argv[1]).read())
soma_t = sm_.soma_t
extension_t = xt_.extension_t
diff --git a/main_prm.py b/parameters.py
similarity index 100%
rename from main_prm.py
rename to parameters.py
diff --git a/auto-py-to-exe.json b/standalone/auto-py-to-exe.json
similarity index 100%
rename from auto-py-to-exe.json
rename to standalone/auto-py-to-exe.json
diff --git a/setup.py b/standalone/setup.py
similarity index 100%
rename from setup.py
rename to standalone/setup.py
diff --git a/standalone-2019-12-06.zip b/standalone/standalone-2019-12-06.zip
similarity index 100%
rename from standalone-2019-12-06.zip
rename to standalone/standalone-2019-12-06.zip