Merge f9c5cf6 into 03d208f

.gitignore

@@ -8,8 +8,6 @@ CMakeFiles
 ITKIOFactoryRegistration
 cmake_install.cmake
 
-/make/itk/itk-config.mk.gen
-
 /ComputeFrangi/Makefile
 /ComputeFrangi/ComputeEigenValue_Sorted_Frangi
 

.travis.yml

@@ -1,16 +1,19 @@
 language: c
-sudo: false
+os: linux
+dist: trusty
+sudo: required
+addons:
+  apt:
+    packages:
+      - libvtk6-dev
+      - libinsighttoolkit4-dev
 env:
   - PROD=1 TEST_VERBOSE=1
   - PROD=1 TEST_VERBOSE=1 COVERAGE=1
 compiler:
   - gcc
   - clang
 before_install:
-  - mkdir build
-  - git clone https://github.com/CBL-ORION/ITK-builds.git external/ITK
-  # path for CMake to search for built ITK
-  - export ITK_DIR="`( cd external/ITK/usr/local/lib/cmake/ITK-4.9 && pwd )`"
   - if [ "$COVERAGE" = 1 ]; then pip install --user cpp-coveralls;  fi
 script:
   # version of gcc on Travis-CI does not support ASAN flag

Makefile

@@ -7,6 +7,14 @@ BIN_SRC_CONFIG.c   := $(SRCDIR)/config/configurator.c
 BIN_SRC.cc  := $(SRCDIR)/compute-filter/ComputeFilter.cxx $(SRCDIR)/subsample-volume/SubsampleVolume.cxx
 #LIB_SRC.c  := $(LIBDIR)/ndarray/ndarray3.c  # $(LIBDIR)/filter/hdaf/Makefilter.c
 LIB_SRC.c  := $(shell find $(LIBDIR) -path lib/t -prune -o  -type f -name '*.c' \! -name '*_impl.c' -print)
+LIB_SRC.cc  := $(shell find $(LIBDIR) \
+	-path lib/t -prune -o \
+	-path lib/vaa3d-plugin -prune -o \
+	-path lib/config -prune -o \
+	-path lib/filter/vesselness -prune -o \
+	-path lib/integration -prune -o \
+	-type f -name '*.cxx' -print)
+$(echo $(LIB_SRC.cc))
 TEST.c := $(shell find $(TESTDIR) -path lib/t/liborion3mat -prune -o -type f -name "*.c" -print)
 ifdef FEAT_LIBORION3MAT
 TEST.c  += $(TESTDIR)/liborion3mat/test.c
@@ -16,11 +24,13 @@ include make/autodep.mk
 
 ## Target files
 OBJ_PATHSUBST.c  = $(patsubst $(LIBDIR)/%.c,$(BUILDDIR)/%.o,$(1))
+OBJ_PATHSUBST.cc  = $(patsubst $(LIBDIR)/%.cxx,$(BUILDDIR)/%.o,$(1))
 TEST_PATHSUBST.c = $(patsubst $(TESTDIR)/%.c,$(BUILDTESTDIR)/%$(EXEEXT),$(1))
 BIN_PATHSUBST.c  = $(patsubst $(SRCDIR)/%.c,$(BINDIR)/%$(EXEEXT),$(1))
 BIN_PATHSUBST.cc  = $(patsubst $(SRCDIR)/%.cxx,$(BINDIR)/%$(EXEEXT),$(1))
 
-MKDIR_BUILD = mkdir -p `dirname $(call OBJ_PATHSUBST.c,$<)`
+MKDIR_BUILD.c = mkdir -p `dirname $(call OBJ_PATHSUBST.c,$<)`
+MKDIR_BUILD.cc = mkdir -p `dirname $(call OBJ_PATHSUBST.cc,$<)`
 MKDIR_BIN   = mkdir -p `dirname $(call BIN_PATHSUBST.c,$<)`
 MKDIR_BUILDTESTDIR = mkdir -p `dirname $(call TEST_PATHSUBST.c,$<)`
 
@@ -29,9 +39,11 @@ include make/vaa3d-plugin-config.mk
 include make/test-config.mk
 include make/liborion3mat-config.mk
 include make/misc-config.mk
+include make/liborion-config.mk
 
 ## Generate the targets
 LIB_OBJ:= $(call OBJ_PATHSUBST.c,$(LIB_SRC.c))
+LIB_OBJ+= $(call OBJ_PATHSUBST.cc,$(LIB_SRC.cc))
 
 TEST_OBJ:= $(call TEST_PATHSUBST.c,$(TEST.c))
 

debian-packages

@@ -0,0 +1,2 @@
+cmake
+libinsighttoolkit4-dev

lib/config/platform.h

@@ -0,0 +1,41 @@
+#ifndef CONFIG_PLATFORM_H
+#define CONFIG_PLATFORM_H 1
+
+/* define OS_* macros */
+#if defined(_WIN32)
+  #define OS_WINDOWS 1
+#else
+  #define OS_UNIX 1
+#endif
+
+
+/* Platform specific settings */
+
+/* On Windows */
+#if OS_WINDOWS
+  #if _MSC_VER && !__INTEL_COMPILER
+    /* Microsoft compiler */
+
+    /* needed for the secure rand_s function */
+    #define _CRT_RAND_S 1
+    #include <stdlib.h>
+
+  #endif
+
+  /* the extension used for executables */
+  #ifndef PLATFORM_EXE_EXT
+    #define PLATFORM_EXE_EXT ".exe"
+  #endif /* PLATFORM_EXE_EXT */
+
+#endif /* OS_WINDOWS */
+
+#if OS_UNIX
+
+  /* the extension used for executables */
+  #ifndef PLATFORM_EXE_EXT
+    #define PLATFORM_EXE_EXT ""
+  #endif /* PLATFORM_EXE_EXT */
+
+#endif /* OS_UNIX */
+
+#endif /* CONFIG_PLATFORM_H */

lib/container/array_impl.c

@@ -67,3 +67,8 @@ void TYPED_NAME(array_resize) ( TYPED_NAME(array)* array, size_t new_capacity )
 	RESIZE_COUNT(array->data, TYPE, new_capacity);
 	array->capacity = new_capacity;
 }
+
+void TYPED_NAME(array_reset) ( TYPED_NAME(array)* array ) {
+	assert( array != NULL );
+	array->length = 0;
+}

lib/container/array_impl.h

@@ -33,6 +33,8 @@ extern void TYPED_NAME(array_add) ( TYPED_NAME(array)* array, TYPE data );
 
 /* ==== memory management ==== */
 extern void TYPED_NAME(array_resize) ( TYPED_NAME(array)* array, size_t new_capacity );
+/* Resets the length of the array to zero. Does not free elements! */
+extern void TYPED_NAME(array_reset) ( TYPED_NAME(array)* array );
 
 #ifdef __cplusplus
 };

lib/io/path/path.c

@@ -21,7 +21,6 @@ orion_filepath* orion_filepath_new_from_string(const char* fp_string) {
 	char* fp_string_norm;
 	size_t fp_string_norm_len;
 	safe_malloc_and_strcpy(&fp_string_norm, fp_string);
-	size_t fp_string_len = strlen( fp_string );
 
 	orion_filepath* fp = _orion_filepath_init();
 

lib/io/path/path.h

@@ -39,6 +39,7 @@ extern char* orion_filepath_to_string(orion_filepath* fp);
 extern orion_filepath* orion_filepath_base( orion_filepath* fp );
 extern orion_filepath* orion_filepath_dir( orion_filepath* fp );
 
+extern orion_filepath* orion_filepath_cat( orion_filepath* fp1, orion_filepath* fp2 );
 extern orion_filepath* orion_filepath_sibling( orion_filepath* base, orion_filepath* relative);
 
 #ifdef __cplusplus

lib/kitchen-sink/01_Segmentation/dendrites_main/DetectTrainingSet/multiscaleLaplacianFilter.c

@@ -6,7 +6,9 @@
 #include "kitchen-sink/01_Segmentation/dendrites_main/DetectTrainingSet/multiscaleLaplacianFilter.h"
 
 #include "util/util.h"
+#include "simple-log/simplelog.h"
 #include "kitchen-sink/01_Segmentation/dendrites_main/DetectTrainingSet/IsotropicFilter/Makefilter.h"
+#include "ndarray/ndarray3_fft.h"
 
 
 #define ORION_LAPLACIAN_HDAF_APPROX_DEGREE 60
@@ -25,6 +27,9 @@ void orion_multiscale_laplacian_output_free(orion_multiscale_laplacian_output* r
 	if( r->max_response_at_scale_idx )
 		ndarray3_free( r->max_response_at_scale_idx );
 
+	if( r->is_background )
+		ndarray3_free( r->is_background );
+
 	free(r);
 }
 
@@ -33,6 +38,7 @@ orion_multiscale_laplacian_output* orion_multiscale_laplacian_output_new() {
 	NEW(r, orion_multiscale_laplacian_output);
 	r->laplacian = NULL;
 	r->max_response_at_scale_idx = NULL;
+	r->is_background = NULL;
 	return r;
 }
 
@@ -52,8 +58,8 @@ orion_multiscale_laplacian_output* orion_multiscaleLaplacianFilter(
 	int n = ORION_LAPLACIAN_HDAF_APPROX_DEGREE;
 
 	/* Fourier transform the input image */
-	ndarray3* v_fft = ndarray3_fftn(input_volume);
-
+	LOG_DEBUG("Computing Fourier transform of input volume");
+	ndarray3_complex* v_fft = ndarray3_fftn_real(input_volume);
 
 	output->laplacian
 		= ndarray3_new_with_size_from_ndarray3(input_volume);
@@ -62,45 +68,67 @@ orion_multiscale_laplacian_output* orion_multiscaleLaplacianFilter(
 
 	int nx = v_fft->sz[0], ny = v_fft->sz[1], nz = v_fft->sz[2];
 
-	/* used to hold the Laplacian for a given scale */
-	ndarray3* cur_scale_laplacian =
-		ndarray3_new_with_size_from_ndarray3(input_volume);
+	/* used to hold the Laplacian response in the Fourier domain for a
+	 * given scale
+	 */
+	ndarray3_complex* cur_scale_laplacian =
+		ndarray3_complex_new(
+			input_volume->sz[0],
+			input_volume->sz[1],
+			input_volume->sz[2]
+		);
 
 	size_t lap_scale_len = array_length_float(laplacian_scales);
 	bool recorded_maximum_scale = false;
-	for( int lap_idx = 0; lap_idx < lap_scale_len; lap_idx++ ) {
+	for( size_t lap_idx = 0; lap_idx < lap_scale_len; lap_idx++ ) {
+		/* compute the Laplacian for each scale */
 		float laplacian_scale_factor = array_get_float( laplacian_scales, lap_idx );
+		LOG_INFO("Computing Laplacian using scale: %f", laplacian_scale_factor);
 
 		float64 norm_const = _orion_ConstantToNormalizeFilter(laplacian_scale_factor);
 		ndarray3* filt = orion_Makefilter( nx,ny,nz,
 			n, laplacian_scale_factor, orion_Makefilter_FLAG_A);
 
-
+		LOG_DEBUG("Applying Laplacian filter");
 		NDARRAY3_LOOP_OVER_START( filt, i,j,k) {
-			/* normalise the generated Laplacian filter */
+			/* normalise the generated Laplacian filter
+			 *
+			 *     filt /= norm_const;
+			 */
 			ndarray3_set(filt, i,j,k,
 				ndarray3_get(filt, i,j,k)
 				/ norm_const
 			);
 
 			/* Apply the filter to the input volume using pointwise
-			 * multiplication in the Fourier domain. */
-			ndarray3_set(cur_scale_laplacian, i,j,k,
-				  ndarray3_get(v_fft, i,j,k)
-				* ndarray3_get(filt,         i,j,k)
+			 * multiplication in the Fourier domain.
+			 *
+			 *     cur_scale_laplacian = v_fft .* filt;
+			 *
+			 * NOTE: This is a complex multiplied by a scalar.
+			 */
+			ndarray3_complex_set(cur_scale_laplacian, i,j,k,
+				  ndarray3_complex_get(v_fft,     i,j,k)
+				* ndarray3_get(filt,              i,j,k)
 			);
 		} NDARRAY3_LOOP_OVER_END;
 
 		/* Now apply the inverse Fourier transform to bring the
-		 * filtered volume back to the spatial domain */
+		 * filtered volume back to the spatial domain.
+		 *
+		 * Again, note that the response is real because the input data
+		 * is real and we treat the data in `cur_scale_laplacian` after
+		 * applying the filter as conjugate symmetric.
+		 * */
+		LOG_DEBUG("Inverse Fourier transform for scale %f", laplacian_scale_factor);
 		ndarray3* cur_scale_filt_vol =
-			ndarray3_ifftn_symmetric( cur_scale_laplacian );
+			ndarray3_ifftn_real( cur_scale_laplacian );
 
 		/* check for maximum response across all scales */
 		if( !recorded_maximum_scale ) {
-			NDARRAY3_LOOP_OVER_START( cur_scale_laplacian, i,j,k) {
+			NDARRAY3_LOOP_OVER_START( cur_scale_filt_vol, i,j,k) {
 				ndarray3_set( output->laplacian, i,j,k,
-					ndarray3_get(cur_scale_laplacian, i,j,k)
+					ndarray3_get(cur_scale_filt_vol, i,j,k)
 				    );
 				ndarray3_set( output->max_response_at_scale_idx, i,j,k,
 					lap_idx
@@ -110,13 +138,37 @@ orion_multiscale_laplacian_output* orion_multiscaleLaplacianFilter(
 			recorded_maximum_scale = true;
 		} else {
 			if( p->multiscale ) {
-				WARN_UNIMPLEMENTED;
+				/* Multiscale approach */
+				NDARRAY3_LOOP_OVER_START( cur_scale_filt_vol, i,j,k) {
+					pixel_type cur = ndarray3_get(cur_scale_filt_vol, i,j,k);
+					pixel_type lap = ndarray3_get(output->laplacian,  i,j,k);
+					if(   fabs(cur) >= fabs(lap) ) {
+						/* the current scale is recorded as the maximum at this point */
+						ndarray3_set(output->laplacian, i,j,k, lap);
+						ndarray3_set(output->max_response_at_scale_idx, i,j,k,  lap_idx);
+					}
+				} NDARRAY3_LOOP_OVER_END;
 			} else {
-				WARN_UNIMPLEMENTED;
+				/* ISBI 2014 */
+				/* Check if the values themselves are larger
+				 * than the current maximum
+				 */
+				NDARRAY3_LOOP_OVER_START( cur_scale_filt_vol, i,j,k) {
+					pixel_type cur = ndarray3_get(cur_scale_filt_vol, i,j,k);
+					pixel_type lap = ndarray3_get(output->laplacian, i,j,k);
+					if( cur > lap ) {
+						ndarray3_set(output->laplacian, i,j,k, cur);
+					}
+				} NDARRAY3_LOOP_OVER_END;
 			}
 		}
+
+		ndarray3_free(cur_scale_filt_vol);
+		ndarray3_free(filt);
 	}
-	ndarray3_free(cur_scale_laplacian);
+
+	ndarray3_complex_free(cur_scale_laplacian);
+	ndarray3_complex_free(v_fft);
 
 	return output;
 }

lib/kitchen-sink/01_Segmentation/dendrites_main/DetectTrainingSet/multiscaleLaplacianFilter.h

@@ -7,20 +7,35 @@
 
 /* local headers */
 #include "ndarray/ndarray3.h"
+#include "ndarray/ndarray3_complex.h"
 #include "container/array.h"
 #include "param/segmentation.h"
 
 /* structs, enums */
 typedef struct {
 	/* each element contains the maximum response of the Laplacian over the
 	 * scales in `laplacian_scales`.
+	 *
+	 * The elements of this are not complex because we are operating on
+	 * only real-data and this means that the Fourier domain values are
+	 * conjugate symmetric.
 	 */
 	ndarray3* laplacian;
 
 	/* each element of scale_for_max_response contains the index of the
 	 * scale in `laplacian_scales` for which the response is maximised
+	 *
+	 * NOTE: Elements of this are of type `size_t`, but are stored in a
+	 * `pixel_type` for convenience.
 	 */
 	ndarray3* max_response_at_scale_idx;
+
+	/* Each element indicates if the Laplacian response at that voxel
+	 * should be classified as in in the background.
+	 *
+	 * Element is either 0.0 or 1.0.
+	 */
+	ndarray3* is_background;
 } orion_multiscale_laplacian_output;
 
 #ifdef __cplusplus
@@ -33,6 +48,9 @@ extern orion_multiscale_laplacian_output* orion_multiscaleLaplacianFilter(
 		array_float* laplacian_scales,
 		orion_segmentation_param* p );
 
+extern orion_multiscale_laplacian_output* orion_multiscale_laplacian_output_new();
+extern void orion_multiscale_laplacian_output_free(orion_multiscale_laplacian_output* r);
+
 #ifdef __cplusplus
 };
 #endif /* __cplusplus */