Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 128b8951 authored by DEBREUVE Eric's avatar DEBREUVE Eric
Browse files

moved frangi in independent folder

parent 3b8c5fc0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 28 additions and 1859 deletions
brick/component/extension.py
+ 12
− 65
View file @ 128b8951
...@@ -31,7 +31,7 @@ ...@@ -31,7 +31,7 @@
from __future__ import annotations 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_ import brick.processing.map_labeling as ml_
from brick.component.glial_cmp import glial_cmp_t from brick.component.glial_cmp import glial_cmp_t
from brick.general.type import array_t, site_h from brick.general.type import array_t, site_h
...@@ -48,35 +48,6 @@ from scipy import ndimage as im_ ...@@ -48,35 +48,6 @@ from scipy import ndimage as im_
min_area_c = 100 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): class extension_t(glial_cmp_t):
# #
# soma_uid: connected to a soma somewhere upstream # soma_uid: connected to a soma somewhere upstream
...@@ -173,7 +144,7 @@ class extension_t(glial_cmp_t): ...@@ -173,7 +144,7 @@ class extension_t(glial_cmp_t):
# #
# import os.path as ph_ # import os.path as ph_
# if ph_.exists("./frangi.npz"): # if ph_.exists("./__runtime__/frangi.npz"):
# print("/!\\ Reading from precomputed data file") # print("/!\\ Reading from precomputed data file")
# loaded = np_.load("./frangi.npz") # loaded = np_.load("./frangi.npz")
# enhanced_img = loaded["enhanced_img"] # enhanced_img = loaded["enhanced_img"]
...@@ -185,33 +156,18 @@ class extension_t(glial_cmp_t): ...@@ -185,33 +156,18 @@ class extension_t(glial_cmp_t):
image, size=2, mode="constant", cval=0.0, origin=0 image, size=2, mode="constant", cval=0.0, origin=0
) )
preprocessed_img = preprocessed_img.astype(np_.float32) enhanced_img, scale_map = fg_.FrangiEnhancement(
enhanced_img = np_.empty(preprocessed_img.shape, dtype=np_.float32) preprocessed_img,
scale_map = np_.empty(preprocessed_img.shape, dtype=np_.float32) scale_range=(0.1, 3.1),
scale_step=1.0,
vesselness_C( alpha=0.8,
preprocessed_img.ctypes.data, beta=0.5,
*preprocessed_img.shape, frangi_c=500.0,
enhanced_img.ctypes.data, bright_on_dark=True,
0.1, in_parallel=in_parallel,
3.2, method="c",
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( # np_.savez_compressed(
# "./runtime/frangi.npz", enhanced_img=enhanced_img, scale_map=scale_map # "./runtime/frangi.npz", enhanced_img=enhanced_img, scale_map=scale_map
# ) # )
...@@ -256,15 +212,6 @@ class extension_t(glial_cmp_t): ...@@ -256,15 +212,6 @@ class extension_t(glial_cmp_t):
return result.astype(np_.int8) 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: def __HysterisisImage__(image: array_t, low: float, high: float) -> array_t:
# #
# low = 0.02 # low = 0.02
......
brick/component/glial_cmp.py
+ 14
− 0
View file @ 128b8951
...@@ -91,3 +91,17 @@ class glial_cmp_t: ...@@ -91,3 +91,17 @@ class glial_cmp_t:
def BackReferenceSoma(self, glial_cmp: glial_cmp_t) -> None: def BackReferenceSoma(self, glial_cmp: glial_cmp_t) -> None:
# #
raise NotImplementedError 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
brick/component/soma.py
+ 0
− 9
View file @ 128b8951
...@@ -238,12 +238,3 @@ def __PointsCloseTo__( ...@@ -238,12 +238,3 @@ def __PointsCloseTo__(
points_as_picker = None points_as_picker = None
return points, points_as_picker 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
brick/processing/dijkstra_1_to_n.py
+ 1
− 1
View file @ 128b8951
../../../../../../brick/def/dijkstra-img/dijkstra_1_to_n.py /home/eric/Code/brick/def/dijkstra-img/dijkstra_1_to_n.py
\ No newline at end of file \ No newline at end of file
brick/processing/frangi3.py 0 → 120000
+ 1
− 0
View file @ 128b8951
/home/eric/Code/brick/def/frangi3/frangi_py/frangi3.py
\ No newline at end of file
brick/processing/frangi_c/makefile deleted 100644 → 0
+ 0
− 43
View file @ 3b8c5fc0
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)
brick/processing/frangi_c/src/__coverage__.txt deleted 100644 → 0
+ 0
− 7
View file @ 3b8c5fc0
gcc -Wall -fprofile-arcs -ftest-coverage cov.c
Then run executable
gcov cov.c
Then check: cov.c.gcov
brick/processing/frangi_c/src/__formatting__.txt deleted 100644 → 0
+ 0
− 1
View file @ 3b8c5fc0
astyle --style=stroustrup --recursive src/*.cpp,*.h
brick/processing/frangi_c/src/data3D.h deleted 100644 → 0
+ 0
− 338
View file @ 3b8c5fc0
#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;
}
brick/processing/frangi_c/src/detector.cpp deleted 100644 → 0
+ 0
− 181
View file @ 3b8c5fc0
#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;
}
brick/processing/frangi_c/src/detector.h deleted 100644 → 0
+ 0
− 31
View file @ 3b8c5fc0
#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
brick/processing/frangi_c/src/eigen.cpp deleted 100644 → 0
+ 0
− 576
View file @ 3b8c5fc0
#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;
}
*/
brick/processing/frangi_c/src/eigen.h deleted 100644 → 0
+ 0
− 37
View file @ 3b8c5fc0
/* 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] ) ;
brick/processing/frangi_c/src/kernel3D.cpp deleted 100644 → 0
+ 0
− 96
View file @ 3b8c5fc0
#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;
}
brick/processing/frangi_c/src/kernel3D.h deleted 100644 → 0
+ 0
− 78
View file @ 3b8c5fc0
#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;
}
};
brick/processing/frangi_c/src/main.cpp deleted 100644 → 0
+ 0
− 64
View file @ 3b8c5fc0
#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;
}
brick/processing/frangi_c/src/processing.cpp deleted 100644 → 0
+ 0
− 156
View file @ 3b8c5fc0
#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 ) );
}
}
}
brick/processing/frangi_c/src/processing.h deleted 100644 → 0
+ 0
− 14
View file @ 3b8c5fc0
#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 );
}
brick/processing/frangi_c/src/processing_code.h deleted 100644 → 0
+ 0
− 121
View file @ 3b8c5fc0
#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]);
}
brick/processing/frangi_c/src/smart_assert.h deleted 100644 → 0
+ 0
− 41
View file @ 3b8c5fc0
/* 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment