From 128b89518ae84d9191c18ab7ad2ea4bb7f00b5cb Mon Sep 17 00:00:00 2001
From: DEBREUVE Eric <eric.debreuve@inria.fr>
Date: Tue, 7 Jan 2020 11:23:23 +0100
Subject: [PATCH] moved frangi in independent folder
---
brick/component/extension.py | 77 +--
brick/component/glial_cmp.py | 14 +
brick/component/soma.py | 9 -
brick/processing/dijkstra_1_to_n.py | 2 +-
brick/processing/frangi3.py | 1 +
brick/processing/frangi_c/makefile | 43 --
.../processing/frangi_c/src/__coverage__.txt | 7 -
.../frangi_c/src/__formatting__.txt | 1 -
brick/processing/frangi_c/src/data3D.h | 338 ----------
brick/processing/frangi_c/src/detector.cpp | 181 ------
brick/processing/frangi_c/src/detector.h | 31 -
brick/processing/frangi_c/src/eigen.cpp | 576 ------------------
brick/processing/frangi_c/src/eigen.h | 37 --
brick/processing/frangi_c/src/kernel3D.cpp | 96 ---
brick/processing/frangi_c/src/kernel3D.h | 78 ---
brick/processing/frangi_c/src/main.cpp | 64 --
brick/processing/frangi_c/src/processing.cpp | 156 -----
brick/processing/frangi_c/src/processing.h | 14 -
.../processing/frangi_c/src/processing_code.h | 121 ----
brick/processing/frangi_c/src/smart_assert.h | 41 --
brick/processing/frangi_c/src/type_info.h | 81 ---
brick/processing/frangi_c/src/types.h | 297 ---------
brick/processing/frangi_py/frangi_3d.pxd | 44 --
brick/processing/frangi_py/frangi_3d.py | 292 ---------
nutrimorph.py | 16 +-
25 files changed, 40 insertions(+), 2577 deletions(-)
create mode 120000 brick/processing/frangi3.py
delete mode 100644 brick/processing/frangi_c/makefile
delete mode 100644 brick/processing/frangi_c/src/__coverage__.txt
delete mode 100644 brick/processing/frangi_c/src/__formatting__.txt
delete mode 100644 brick/processing/frangi_c/src/data3D.h
delete mode 100644 brick/processing/frangi_c/src/detector.cpp
delete mode 100644 brick/processing/frangi_c/src/detector.h
delete mode 100644 brick/processing/frangi_c/src/eigen.cpp
delete mode 100644 brick/processing/frangi_c/src/eigen.h
delete mode 100644 brick/processing/frangi_c/src/kernel3D.cpp
delete mode 100644 brick/processing/frangi_c/src/kernel3D.h
delete mode 100644 brick/processing/frangi_c/src/main.cpp
delete mode 100644 brick/processing/frangi_c/src/processing.cpp
delete mode 100644 brick/processing/frangi_c/src/processing.h
delete mode 100644 brick/processing/frangi_c/src/processing_code.h
delete mode 100644 brick/processing/frangi_c/src/smart_assert.h
delete mode 100644 brick/processing/frangi_c/src/type_info.h
delete mode 100644 brick/processing/frangi_c/src/types.h
delete mode 100644 brick/processing/frangi_py/frangi_3d.pxd
delete mode 100644 brick/processing/frangi_py/frangi_3d.py
diff --git a/brick/component/extension.py b/brick/component/extension.py
index 0f8a80b..2141913 100644
--- a/brick/component/extension.py
+++ b/brick/component/extension.py
@@ -31,7 +31,7 @@
from __future__ import annotations
-# import brick.processing.frangi_py.frangi_3d as fg_
+import brick.processing.frangi3 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
@@ -48,35 +48,6 @@ 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
@@ -173,7 +144,7 @@ class extension_t(glial_cmp_t):
#
# import os.path as ph_
- # if ph_.exists("./frangi.npz"):
+ # if ph_.exists("./__runtime__/frangi.npz"):
# print("/!\\ Reading from precomputed data file")
# loaded = np_.load("./frangi.npz")
# enhanced_img = loaded["enhanced_img"]
@@ -185,33 +156,18 @@ class extension_t(glial_cmp_t):
image, size=2, mode="constant", cval=0.0, origin=0
)
- 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.1),
+ scale_step=1.0,
+ alpha=0.8,
+ beta=0.5,
+ frangi_c=500.0,
+ bright_on_dark=True,
+ in_parallel=in_parallel,
+ method="c",
)
- # 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
# )
@@ -256,15 +212,6 @@ class extension_t(glial_cmp_t):
return result.astype(np_.int8)
-def NormalizedImage(image: array_t) -> array_t:
- #
- middle_values = image[np_.logical_and(image > 0.0, image < image.max())]
- image_mean = middle_values.mean()
- result = image / image_mean
-
- return result
-
-
def __HysterisisImage__(image: array_t, low: float, high: float) -> array_t:
#
# low = 0.02
diff --git a/brick/component/glial_cmp.py b/brick/component/glial_cmp.py
index a4584e3..be5fa89 100644
--- a/brick/component/glial_cmp.py
+++ b/brick/component/glial_cmp.py
@@ -91,3 +91,17 @@ class glial_cmp_t:
def BackReferenceSoma(self, glial_cmp: glial_cmp_t) -> None:
#
raise NotImplementedError
+
+
+def NormalizedImage(image: array_t) -> array_t:
+ #
+ print("This normalization does not bring anything; left as is for now to avoid the need for changing prms")
+ nonextreme_values = image[np_.logical_and(image > 0.0, image < image.max())]
+
+ if nonextreme_values.size > 0:
+ nonextreme_avg = nonextreme_values.mean()
+ result = image.astype(np_.float32) / nonextreme_avg
+ else:
+ result = image.astype(np_.float32)
+
+ return result
diff --git a/brick/component/soma.py b/brick/component/soma.py
index 7bdedba..3388dfb 100644
--- a/brick/component/soma.py
+++ b/brick/component/soma.py
@@ -238,12 +238,3 @@ def __PointsCloseTo__(
points_as_picker = None
return points, points_as_picker
-
-
-def NormalizedImage(image: array_t) -> array_t:
- #
- nonextreme_values = image[np_.logical_and(image > 0.0, image < image.max())]
- nonextreme_avg = np_.mean(nonextreme_values)
- result = image / nonextreme_avg
-
- return result
diff --git a/brick/processing/dijkstra_1_to_n.py b/brick/processing/dijkstra_1_to_n.py
index a9cb6fe..be4f731 120000
--- a/brick/processing/dijkstra_1_to_n.py
+++ b/brick/processing/dijkstra_1_to_n.py
@@ -1 +1 @@
-../../../../../../brick/def/dijkstra-img/dijkstra_1_to_n.py
\ No newline at end of file
+/home/eric/Code/brick/def/dijkstra-img/dijkstra_1_to_n.py
\ No newline at end of file
diff --git a/brick/processing/frangi3.py b/brick/processing/frangi3.py
new file mode 120000
index 0000000..74fc34e
--- /dev/null
+++ b/brick/processing/frangi3.py
@@ -0,0 +1 @@
+/home/eric/Code/brick/def/frangi3/frangi_py/frangi3.py
\ No newline at end of file
diff --git a/brick/processing/frangi_c/makefile b/brick/processing/frangi_c/makefile
deleted file mode 100644
index 2f9fb57..0000000
--- a/brick/processing/frangi_c/makefile
+++ /dev/null
@@ -1,43 +0,0 @@
-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
deleted file mode 100644
index 968961f..0000000
--- a/brick/processing/frangi_c/src/__coverage__.txt
+++ /dev/null
@@ -1,7 +0,0 @@
-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
deleted file mode 100644
index 4a0ec64..0000000
--- a/brick/processing/frangi_c/src/__formatting__.txt
+++ /dev/null
@@ -1 +0,0 @@
-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
deleted file mode 100644
index 9eb64e5..0000000
--- a/brick/processing/frangi_c/src/data3D.h
+++ /dev/null
@@ -1,338 +0,0 @@
-#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
deleted file mode 100644
index 6d68f8d..0000000
--- a/brick/processing/frangi_c/src/detector.cpp
+++ /dev/null
@@ -1,181 +0,0 @@
-#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
deleted file mode 100644
index e08e6eb..0000000
--- a/brick/processing/frangi_c/src/detector.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#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
deleted file mode 100644
index 71ab97d..0000000
--- a/brick/processing/frangi_c/src/eigen.cpp
+++ /dev/null
@@ -1,576 +0,0 @@
-#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
deleted file mode 100644
index 1585f40..0000000
--- a/brick/processing/frangi_c/src/eigen.h
+++ /dev/null
@@ -1,37 +0,0 @@
-/* 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
deleted file mode 100644
index 61c440f..0000000
--- a/brick/processing/frangi_c/src/kernel3D.cpp
+++ /dev/null
@@ -1,96 +0,0 @@
-#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
deleted file mode 100644
index be82fb6..0000000
--- a/brick/processing/frangi_c/src/kernel3D.h
+++ /dev/null
@@ -1,78 +0,0 @@
-#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
deleted file mode 100644
index 2158409..0000000
--- a/brick/processing/frangi_c/src/main.cpp
+++ /dev/null
@@ -1,64 +0,0 @@
-#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
deleted file mode 100644
index 984d645..0000000
--- a/brick/processing/frangi_c/src/processing.cpp
+++ /dev/null
@@ -1,156 +0,0 @@
-#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
deleted file mode 100644
index ae73611..0000000
--- a/brick/processing/frangi_c/src/processing.h
+++ /dev/null
@@ -1,14 +0,0 @@
-#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
deleted file mode 100644
index 6221073..0000000
--- a/brick/processing/frangi_c/src/processing_code.h
+++ /dev/null
@@ -1,121 +0,0 @@
-#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
deleted file mode 100644
index a6d38fd..0000000
--- a/brick/processing/frangi_c/src/smart_assert.h
+++ /dev/null
@@ -1,41 +0,0 @@
-/* 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
deleted file mode 100644
index b339e18..0000000
--- a/brick/processing/frangi_c/src/type_info.h
+++ /dev/null
@@ -1,81 +0,0 @@
-#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
deleted file mode 100644
index 3b64784..0000000
--- a/brick/processing/frangi_c/src/types.h
+++ /dev/null
@@ -1,297 +0,0 @@
-#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/brick/processing/frangi_py/frangi_3d.pxd b/brick/processing/frangi_py/frangi_3d.pxd
deleted file mode 100644
index f5e1375..0000000
--- a/brick/processing/frangi_py/frangi_3d.pxd
+++ /dev/null
@@ -1,44 +0,0 @@
-cimport numpy as np_
-
-
-ctypedef np_.float64_t real_pix_v_t
-#ctypedef np_.ndarray[real_pix_v_t, ndim=3] real_img_t
-
-
-#cpdef tuple FrangiEnhancement(
-# np_.ndarray[real_pix_v_t, ndim=2] img,
-# tuple scale_range,
-# double scale_step = *,
-# double alpha = *,
-# double beta = *,
-# double frangi_c = *,
-# bint bright_on_dark = *,
-# bint in_parallel = *
-#)
-
-
-cdef np_.ndarray[real_pix_v_t, ndim=3] FrangiEnhancementOneScale(
- np_.ndarray[real_pix_v_t, ndim=3] img,
- double scale,
- double alpha = *,
- double beta = *,
- double frangi_c = *,
- bint bright_on_dark = *,
- object shared_enhanced = *
-)
-
-
-#cdef tuple[np_.ndarray[real_pix_v_t, ndim=3],np_.ndarray[real_pix_v_t, ndim=3],np_.ndarray[real_pix_v_t, ndim=3],np_.ndarray[real_pix_v_t, ndim=3],np_.ndarray[real_pix_v_t, ndim=3],np_.ndarray[real_pix_v_t, ndim=3]] HessianMatrix(
-# np_.ndarray[real_pix_v_t, ndim=3] img, double scale
-#)
-
-
-#cdef tuple[np_.ndarray[real_pix_v_t, ndim=3],np_.ndarray[real_pix_v_t, ndim=3],np_.ndarray[real_pix_v_t, ndim=3]] EigenValues(
-# np_.ndarray[real_pix_v_t, ndim=3] Dxx, np_.ndarray[real_pix_v_t, ndim=3] Dxy, np_.ndarray[real_pix_v_t, ndim=3] Dxz, np_.ndarray[real_pix_v_t, ndim=3] Dyy, np_.ndarray[real_pix_v_t, ndim=3] Dyz, np_.ndarray[real_pix_v_t, ndim=3] Dzz
-#)
-
-
-cdef np_.ndarray[real_pix_v_t, ndim=3] SortedByAbsValues(np_.ndarray[real_pix_v_t, ndim=3] array, int axis = *)
-
-
-cdef np_.ndarray[real_pix_v_t, ndim=3] CarefullDivision(np_.ndarray[real_pix_v_t, ndim=3] array_1, np_.ndarray[real_pix_v_t, ndim=3] array_2)
diff --git a/brick/processing/frangi_py/frangi_3d.py b/brick/processing/frangi_py/frangi_3d.py
deleted file mode 100644
index 9fbff6b..0000000
--- a/brick/processing/frangi_py/frangi_3d.py
+++ /dev/null
@@ -1,292 +0,0 @@
-# Copyright CNRS/Inria/UNS
-# Contributor(s): Eric Debreuve (since 2019)
-#
-# eric.debreuve@cnrs.fr
-#
-# This software is governed by the CeCILL license under French law and
-# abiding by the rules of distribution of free software. You can use,
-# modify and/ or redistribute the software under the terms of the CeCILL
-# license as circulated by CEA, CNRS and INRIA at the following URL
-# "http://www.cecill.info".
-#
-# As a counterpart to the access to the source code and rights to copy,
-# modify and redistribute granted by the license, users are provided only
-# with a limited warranty and the software's author, the holder of the
-# economic rights, and the successive licensors have only limited
-# liability.
-#
-# In this respect, the user's attention is drawn to the risks associated
-# with loading, using, modifying and/or developing or reproducing the
-# software by the user in light of its specific status of free software,
-# that may mean that it is complicated to manipulate, and that also
-# therefore means that it is reserved for developers and experienced
-# professionals having in-depth computer knowledge. Users are therefore
-# encouraged to load and test the software's suitability as regards their
-# requirements in conditions enabling the security of their systems and/or
-# data to be ensured and, more generally, to use and operate it in the
-# same conditions as regards security.
-#
-# The fact that you are presently reading this means that you have had
-# knowledge of the CeCILL license and that you accept its terms.
-
-'''
-Follows: https://github.com/ellisdg/frangi3d
-'''
-
-from brick.general.type import array_t
-
-import multiprocessing as pl_
-from ctypes import c_float as c_float_t
-from multiprocessing.sharedctypes import Array as shared_array_t
-from typing import Optional, Tuple
-
-import itk as ik_
-import numpy as np_
-from scipy import ndimage as im_
-
-
-_hessian_matrices = None
-
-
-def FrangiEnhancement(
- img: array_t,
- scale_range: Tuple[float, float],
- scale_step: float = 2.0,
- alpha: float = 0.5,
- beta: float = 0.5,
- frangi_c: float = 500.0,
- bright_on_dark: bool = True,
- in_parallel: bool = False,
-) -> Tuple[array_t, array_t]:
- #
- global _hessian_matrices
-
- img = img.astype(np_.float32, copy=False)
- scales = np_.linspace(
- scale_range[0],
- scale_range[1],
- num=np_.ceil((scale_range[1] - scale_range[0]) / scale_step) + 1,
- )
-
- if in_parallel and (pl_.get_start_method(allow_none=False) == "fork"):
- _hessian_matrices = None
- FrangiEnhancement_fct = _FrangiEnhancementParallel
- else:
- _hessian_matrices = np_.empty((*img.shape, 3, 3), dtype=np_.float32)
- FrangiEnhancement_fct = _FrangiEnhancementSequential
-
- enhanced, scale_map = FrangiEnhancement_fct(
- img,
- scales,
- alpha=alpha,
- beta=beta,
- frangi_c=frangi_c,
- bright_on_dark=bright_on_dark,
- )
-
- _hessian_matrices = None
-
- return enhanced, scale_map
-
-
-def _FrangiEnhancementSequential(
- img: array_t,
- scales: array_t,
- alpha: float = 0.5,
- beta: float = 0.5,
- frangi_c: float = 500.0,
- bright_on_dark: bool = True,
-) -> Tuple[array_t, array_t]:
- #
- enhanced = None
- scale_map = None
-
- for s_idx, scale in enumerate(scales):
- local_enhanced = _FrangiEnhancementOneScale(
- img,
- scale,
- alpha=alpha,
- beta=beta,
- frangi_c=frangi_c,
- bright_on_dark=bright_on_dark,
- )
-
- if s_idx > 0:
- larger_map = local_enhanced > enhanced
- enhanced[larger_map] = local_enhanced[larger_map]
- scale_map[larger_map] = scale
- else:
- enhanced = local_enhanced
- scale_map = np_.full(img.shape, scale, dtype=np_.float32)
-
- return enhanced, scale_map
-
-
-def _FrangiEnhancementParallel(
- img: array_t,
- scales: array_t,
- alpha: float = 0.5,
- beta: float = 0.5,
- frangi_c: float = 500.0,
- bright_on_dark: bool = True,
-) -> Tuple[array_t, array_t]:
- #
- fixed_args = {
- "alpha": alpha,
- "beta": beta,
- "frangi_c": frangi_c,
- "bright_on_dark": bright_on_dark,
- }
-
- shared_enhanced_s = tuple(shared_array_t(c_float_t, img.size) for ___ in scales)
- processes = tuple(
- pl_.Process(
- target=_FrangiEnhancementOneScale,
- args=(img, scale),
- kwargs={**fixed_args, "shared_enhanced": shared_enhanced,},
- )
- for scale, shared_enhanced in zip(scales, shared_enhanced_s)
- )
-
- for process in processes:
- process.start()
- for process in processes:
- process.join()
-
- enhanced = np_.array(shared_enhanced_s[0])
- scale_map = np_.full(img.size, scales[0], dtype=np_.float32)
- for s_idx, shared_enhanced in enumerate(shared_enhanced_s[1:], start=1):
- local_enhanced = np_.array(shared_enhanced)
- larger_map = local_enhanced > enhanced
- enhanced[larger_map] = local_enhanced[larger_map]
- scale_map[larger_map] = scales[s_idx]
-
- enhanced = enhanced.astype(np_.float32, copy=False).reshape(img.shape)
- scale_map = scale_map.astype(np_.float32, copy=False).reshape(img.shape)
-
- return enhanced, scale_map
-
-
-def _FrangiEnhancementOneScale(
- img: array_t,
- scale: float,
- alpha: float = 0.5,
- beta: float = 0.5,
- frangi_c: float = 500.0,
- bright_on_dark: bool = True,
- shared_enhanced: shared_array_t = None,
-) -> Optional[array_t]:
- #
- print(f" scale={scale}")
-
- sorted_eigenvalues = _EigenValuesOfHessianMatrix(img, scale)
-
- # Absolute plate, blob, and background maps
- abs_eig_val = np_.fabs(sorted_eigenvalues)
- plate_map = _CarefullDivision(abs_eig_val[...,1], abs_eig_val[...,2])
- blob_map = _CarefullDivision(abs_eig_val[...,0], np_.sqrt(np_.prod(abs_eig_val[...,1:], axis = -1)))
- bckgnd = np_.sqrt(np_.sum(abs_eig_val ** 2, axis = -1))
-
- # Normalized plate, blob, and background maps
- plate_factor = -2.0 * alpha ** 2
- blob_factor = -2.0 * beta ** 2
- bckgnd_factor = -2.0 * frangi_c ** 2
- plate_map = 1.0 - np_.exp(plate_map ** 2 / plate_factor)
- blob_map = np_.exp(blob_map ** 2 / blob_factor)
- bckgnd = 1.0 - np_.exp(bckgnd ** 2 / bckgnd_factor)
-
- local_enhanced = plate_map * blob_map * bckgnd
- if bright_on_dark:
- eig_condition = np_.any(sorted_eigenvalues[...,1:] > 0.0, axis=-1)
- else:
- eig_condition = np_.any(sorted_eigenvalues[...,1:] < 0.0, axis=-1)
- invalid_condition = np_.logical_not(np_.isfinite(local_enhanced))
- local_enhanced[np_.logical_or(eig_condition, invalid_condition)] = 0.0
-
- if shared_enhanced is None:
- return local_enhanced
- else:
- shared_enhanced[:] = local_enhanced.flatten()[:]
-
-
-def _ITKFrangiEnhancementOneScale(
- img: array_t,
- scale: float,
- alpha: float = 0.5,
- beta: float = 0.5,
- frangi_c: float = 500.0,
- bright_on_dark: bool = True,
- shared_enhanced: shared_array_t = None,
-) -> Optional[array_t]:
- '''
- Do not produce nice resultas as is. Parameter tuning?
- '''
- print(f" scale={scale}")
-
- assert bright_on_dark
-
- img = ik_.image_view_from_array(img)
- hessian_image = ik_.hessian_recursive_gaussian_image_filter(img, sigma=scale)
-
- vesselness_filter = ik_.Hessian3DToVesselnessMeasureImageFilter[ik_.ctype('float')].New()
- vesselness_filter.SetInput(hessian_image)
- # vesselness_filter.SetAlpha1(args.alpha1)
- # vesselness_filter.SetAlpha2(args.alpha2)
-
- # Do not view-version here since vesselness_filter is local
- local_enhanced = ik_.array_from_image(vesselness_filter)
-
- if shared_enhanced is None:
- return local_enhanced
- else:
- shared_enhanced[:] = local_enhanced.flatten()[:]
-
-
-def _EigenValuesOfHessianMatrix(img: array_t, scale: float) -> array_t:
- #
- global _hessian_matrices
-
- if scale > 0.0:
- img = scale ** 2 * im_.gaussian_filter(img, scale)
- else:
- img = scale ** 2 * img
-
- if _hessian_matrices is None:
- hessian_matrices = np_.empty((*img.shape, 3, 3), dtype=np_.float32)
- else:
- hessian_matrices = _hessian_matrices
- Dx, Dy, Dz = np_.gradient(img)
- (
- hessian_matrices[..., 0, 0],
- hessian_matrices[..., 0, 1],
- hessian_matrices[..., 0, 2],
- ) = np_.gradient(Dx)
- hessian_matrices[..., 1, 1], hessian_matrices[..., 1, 2] = np_.gradient(
- Dy, axis=(1, 2)
- )
- hessian_matrices[..., 2, 2] = np_.gradient(Dz, axis=2)
-
- hessian_matrices[..., 1, 0] = hessian_matrices[..., 0, 1]
- hessian_matrices[..., 2, 0] = hessian_matrices[..., 0, 2]
- hessian_matrices[..., 2, 1] = hessian_matrices[..., 1, 2]
-
- eigenvalues = np_.linalg.eigvalsh(hessian_matrices)
-
- # Sorted by abs value
- index = list(np_.ix_(*[np_.arange(size) for size in eigenvalues.shape]))
- index[-1] = np_.fabs(eigenvalues).argsort(axis=-1)
- sorted_eigenvalues = eigenvalues[tuple(index)]
-
- return sorted_eigenvalues
-
-
-def _CarefullDivision(array_1: array_t, array_2: array_t) -> array_t:
- #
- null_map = array_2 == 0.0
- if null_map.any():
- denominator = array_2.copy()
- denominator[null_map] = 1.0e-10
- else:
- denominator = array_2
-
- return array_1 / denominator
diff --git a/nutrimorph.py b/nutrimorph.py
index 41fdf35..e0d1031 100644
--- a/nutrimorph.py
+++ b/nutrimorph.py
@@ -39,6 +39,7 @@
import brick.component.connection as cn_
import brick.component.extension as xt_
+import brick.component.glial_cmp as gl_
import brick.component.soma as sm_
import brick.general.feedback as fb_
# from main_prm import *
@@ -87,13 +88,20 @@ start_time = tm_.time()
# --- Images
image = io_.imread(data_path)
-assert image.shape.__len__() == 3
-image = image[:, 512:, 512:]
+assert image.ndim == 3
+image = image[:, 512:, 512:]#562
img_shape = image.shape
+#
+# max_img = image.max()
+# image.fill(image.min())
+# image[2:3,10:11,:] = max_img
+#
+print(f"IMAGE: s.{img_shape} t.{image.dtype} m.{image.min()} M.{image.max()}")
-image_for_soma = sm_.NormalizedImage(image)
-image_for_ext = xt_.NormalizedImage(image)
+image_for_soma = gl_.NormalizedImage(image)
+image_for_ext = gl_.NormalizedImage(image)
costs = 1.0 / (image + 1.0)
+print(f"NRM-IMG: t.{image_for_soma.dtype} m.{image_for_soma.min():.2f} M.{image_for_soma.max():.2f}")
# --- Initialization
--
GitLab