From 6653afb78547288c18f48bf1f4e798eaa0f97eeb Mon Sep 17 00:00:00 2001 From: Kirai Date: Thu, 19 Apr 2018 11:11:46 +0800 Subject: [PATCH 1/5] [debug] update type checker If people use `numpy.core.memmap.memmap` to save arrays then it will cresh. use 'numpy' as a check point will be ok. --- denseinference/CRFProcessor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/denseinference/CRFProcessor.py b/denseinference/CRFProcessor.py index 6cd8784..14cd579 100644 --- a/denseinference/CRFProcessor.py +++ b/denseinference/CRFProcessor.py @@ -125,7 +125,7 @@ def set_data_and_run(self, img, label): :param label: Continuous label tensor as ndarray. (W, H, D, L), [0, 1] :return: Hard labeled result as ndarray. (W, H, D), [0, L], dtype=int16 """ - if (type(img).__module__ != np.__name__) or (type(label).__module__ != np.__name__): + if (type(img).__module__[:5] != np.__name__) or (type(label).__module__[:5] != np.__name__): raise ValueError('Error. Ndarray expected. Image: (W, H, D), [0, 1]; Label: (W, H, D, L), [0, 1].') if img.ndim != 3: From 9118e605c879ad6fd2da7214eec28a94ee58cc70 Mon Sep 17 00:00:00 2001 From: kirai Date: Tue, 8 May 2018 15:50:35 +0800 Subject: [PATCH 2/5] add a 2D-denseCRF wrapper --- denseinference/CRFProcessor.py | 547 ++++++++----- .../lib/refine_3d/dense_inference.cpp | 748 ++++++++++++------ .../lib/refine_3d/dense_inference.h | 78 +- readme.md | 3 + 4 files changed, 899 insertions(+), 477 deletions(-) diff --git a/denseinference/CRFProcessor.py b/denseinference/CRFProcessor.py index 14cd579..b52bda0 100644 --- a/denseinference/CRFProcessor.py +++ b/denseinference/CRFProcessor.py @@ -1,4 +1,4 @@ -__author__ = 'Marc Bickel' +__author__ = 'Marc Bickel, Wendong Xu' import numpy as np import psutil @@ -7,221 +7,336 @@ class CRF3DProcessor(object): - # # - # - # Constructor - # - # # - def __init__(self, - max_iterations=10, - pos_x_std=3.0, - pos_y_std=3.0, - pos_z_std=3.0, - pos_w=3.0, - bilateral_x_std=60.0, - bilateral_y_std=60.0, - bilateral_z_std=60.0, - bilateral_intensity_std=20.0, - bilateral_w=10.0, - dynamic_z=False, - ignore_memory=False, - verbose=False): - """ - Initializes DenseCRF by Kraehenbuehl to refine volume labeling. - Features are voxel distance and voxel intensity. - - Can be tuned by adjusting the optional parameters of the constructor. - - Expected data input: Image tensor and label tensor. - Both image and label voxels must be normalized to [0, 1] - Image: (W, H, D); Label (W, H, D, C) - - Result: Result tensor. - Result tensor is hard labeled (W, H, D) - - :param max_iterations: maximum crf iterations - :param pos_x_std: DenseCRF Param (float) (3) - :param pos_y_std: DenseCRF Param (float) (3) - :param pos_z_std: DenseCRF Param (float) (3) - :param pos_w: DenseCRF Param (float) (3) - :param bilateral_x_std: DenseCRF Param (float) (60) - :param bilateral_y_std: DenseCRF Param (float) (60) - :param bilateral_z_std: DenseCRF Param (float) (60) - :param bilateral_intensity_std: DenseCRF Param (float) (20) - :param bilateral_w: DenseCRF Param (float) (10) - :param dynamic_z: Auto adjust z-params by image shape (bool) (False) - :param ignore_memory: If true, images requiring to much memory are skipped (bool) (False) - :param verbose: Print lot's of status information (bool) (False) - """ - - # private vars - self.__verbose = verbose - self.__ignore_mem = ignore_memory - - # adaptive z - self.__dynamic_z = dynamic_z - self.__bilateral_z_std = bilateral_z_std - - # init processor - self.__processor = di.DenseCRF3DProcessor(max_iterations, - pos_x_std, - pos_y_std, - pos_z_std, - pos_w, - bilateral_x_std, - bilateral_y_std, - bilateral_z_std, - bilateral_intensity_std, - bilateral_w, - verbose) - - # # - # - # Static memory helpers - # - # # - @staticmethod - def memory_usage(label_shape): - """ - Get anticipated memory requirements for labels_shape. - :param label_shape: 4-dimensional label shape as tuple - :return: required memory in bytes - """ - if len(label_shape) != 4: - raise ValueError('Error. 4-dimensional tensor expected. Got: ' + str(len(label_shape))) - - return di.DenseCRF3DProcessor.get_memory_requirements(label_shape[0], - label_shape[1], - label_shape[2], - label_shape[3]) - - @staticmethod - def memory_free(): - """ - Get current free system memory. - :return: free memory in bytes - """ - return psutil.virtual_memory().free - - @staticmethod - def memory_percent(label_shape): - """ - Get required memory for label_shape as fraction of free virtual memory as percent. - :param label_shape: 4-dimensional label shape as tuple - :return: required memory of free memory in percent - """ - return float(CRF3DProcessor.memory_usage(label_shape)) / float(CRF3DProcessor.memory_free()) * 100.0 - - # # - # - # Single exposed method to run CRF - # - # # - def set_data_and_run(self, img, label): - """ - Pass raw image volume and label tensor of your classifier (4d), run the CRF and receive refined hard labeled - output. All input data and the returned output are ndarrays. - :param img: Normalized input as ndarray. (W, H, D), [0, 1] - :param label: Continuous label tensor as ndarray. (W, H, D, L), [0, 1] - :return: Hard labeled result as ndarray. (W, H, D), [0, L], dtype=int16 - """ - if (type(img).__module__[:5] != np.__name__) or (type(label).__module__[:5] != np.__name__): - raise ValueError('Error. Ndarray expected. Image: (W, H, D), [0, 1]; Label: (W, H, D, L), [0, 1].') - - if img.ndim != 3: - raise ValueError('Error. 3d tensor expected. Got: ' + str(img.ndim)) - - if label.ndim != 4: - raise ValueError('Error. 4d tensor expected. Got: ' + str(label.ndim)) - - # check image to label shape consistency - if (img.shape[0] != label.shape[0]) or (img.shape[1] != label.shape[1]) or (img.shape[2] != label.shape[2]): - raise ValueError('Error. Image shape and label shape inconsistent: ' + - str(img.shape) + ', ' + str(label.shape)) - - # set image - self.__check_and_set_img(img) - - # set label - self.__check_and_set_label(label) - - # run crf - res_arr = self.__run_crf(label.shape) - - return CRF3DProcessor.__prepare_and_return_result(res_arr, img.shape) - - # # - # - # Private convenience methods - # - # # - def __check_and_set_img(self, img): - if img.ndim != 3: - raise ValueError('Error. 3d tensor expected. Got: ' + str(img.ndim)) - - if np.amin(img) < 0 or np.amax(img) > 1: - raise ValueError('Error. Image must be normalized to [0, 1].') - - # pass image to crf - self.__processor.set_image(np.array(img, dtype=np.float64), - img.shape[0], - img.shape[1], - img.shape[2], - 'float64') - - return img.shape - - def __check_and_set_label(self, label): - if label.ndim != 4: - raise ValueError('Error. 4d tensor expected. Got: ' + str(label.ndim)) - - if np.amin(label) < 0 or np.amax(label) > 1: - raise ValueError('Error. Label must be normalized to [0, 1].') - - # pass labels to crf - self.__processor.set_feature(np.array(label, dtype=np.float64), - label.shape[0], - label.shape[1], - label.shape[2], - label.shape[3], - 'float64') - - return label.shape - - @staticmethod - def __prepare_and_return_result(res_arr, shape): - if len(shape) != 3: - raise ValueError('Error. 3-d shape expected. Got: ' + str(len(shape))) - - if len(res_arr) != shape[0] * shape[1] * shape[2]: - raise ValueError('Error. Length of array inconsisten with length of shape (' + - str(len(res_arr)) + ', ' + str(shape[0] * shape[1] * shape[2]) + ').') - - # Reshape with fortran order - return np.asarray(res_arr, dtype=np.int16).reshape(shape, order='F') - - # # - # - # CRF - # - # # - def __run_crf(self, label_shape): - # check memory - if (not self.__ignore_mem) and (self.memory_percent(label_shape) > 100.0): - raise RuntimeError('Current image exceeds available memory. Free vs Needed: (' + - str(round(self.memory_free()/10**9, 2)) + 'GB, ' + - str(round(self.memory_usage(label_shape)/10**9, 2)) + 'GB).') - - if self.__dynamic_z: - xy_mean = np.mean([label_shape[0], label_shape[1]]) - if (xy_mean != label_shape[2]) and (xy_mean != 0): - bi_z_std = self.__bilateral_z_std * (float(label_shape[2]) / xy_mean) - if self.__verbose: - print('Set Bilateral Z Std to ' + str(bi_z_std) + '.') - self.__processor.set_bi_z_std(bi_z_std) - - return self.__processor.calc_res_array() + # # + # + # Constructor + # + # # + def __init__(self, + max_iterations=10, + pos_x_std=3.0, + pos_y_std=3.0, + pos_z_std=3.0, + pos_w=3.0, + bilateral_x_std=60.0, + bilateral_y_std=60.0, + bilateral_z_std=60.0, + bilateral_intensity_std=20.0, + bilateral_w=10.0, + dynamic_z=False, + ignore_memory=False, + verbose=False): + """ + Initializes DenseCRF by Kraehenbuehl to refine volume labeling. + Features are voxel distance and voxel intensity. + + Can be tuned by adjusting the optional parameters of the constructor. + + Expected data input: Image tensor and label tensor. + Both image and label voxels must be normalized to [0, 1] + Image: (W, H, D); Label (W, H, D, C) + + Result: Result tensor. + Result tensor is hard labeled (W, H, D) + + :param max_iterations: maximum crf iterations + :param pos_x_std: DenseCRF Param (float) (3) + :param pos_y_std: DenseCRF Param (float) (3) + :param pos_z_std: DenseCRF Param (float) (3) + :param pos_w: DenseCRF Param (float) (3) + :param bilateral_x_std: DenseCRF Param (float) (60) + :param bilateral_y_std: DenseCRF Param (float) (60) + :param bilateral_z_std: DenseCRF Param (float) (60) + :param bilateral_intensity_std: DenseCRF Param (float) (20) + :param bilateral_w: DenseCRF Param (float) (10) + :param dynamic_z: Auto adjust z-params by image shape (bool) (False) + :param ignore_memory: If true, images requiring to much memory are skipped (bool) (False) + :param verbose: Print lot's of status information (bool) (False) + """ + + # private vars + self.__verbose = verbose + self.__ignore_mem = ignore_memory + + # adaptive z + self.__dynamic_z = dynamic_z + self.__bilateral_z_std = bilateral_z_std + + # init processor + self.__processor = di.DenseCRF3DProcessor(max_iterations, + pos_x_std, + pos_y_std, + pos_z_std, + pos_w, + bilateral_x_std, + bilateral_y_std, + bilateral_z_std, + bilateral_intensity_std, + bilateral_w, + verbose) + + # # + # + # Static memory helpers + # + # # + @staticmethod + def memory_usage(label_shape): + """ + Get anticipated memory requirements for labels_shape. + :param label_shape: 4-dimensional label shape as tuple + :return: required memory in bytes + """ + if len(label_shape) != 4: + raise ValueError('Error. 4-dimensional tensor expected. Got: ' + str(len(label_shape))) + + return di.DenseCRF3DProcessor.get_memory_requirements(label_shape[0], + label_shape[1], + label_shape[2], + label_shape[3]) + + @staticmethod + def memory_free(): + """ + Get current free system memory. + :return: free memory in bytes + """ + return psutil.virtual_memory().free + + @staticmethod + def memory_percent(label_shape): + """ + Get required memory for label_shape as fraction of free virtual memory as percent. + :param label_shape: 4-dimensional label shape as tuple + :return: required memory of free memory in percent + """ + return float(CRF3DProcessor.memory_usage(label_shape)) / float(CRF3DProcessor.memory_free()) * 100.0 + + # # + # + # Single exposed method to run CRF + # + # # + def set_data_and_run(self, img, label): + """ + Pass raw image volume and label tensor of your classifier (4d), run the CRF and receive refined hard labeled + output. All input data and the returned output are ndarrays. + :param img: Normalized input as ndarray. (W, H, D), [0, 1] + :param label: Continuous label tensor as ndarray. (W, H, D, L), [0, 1] + :return: Hard labeled result as ndarray. (W, H, D), [0, L], dtype=int16 + """ + if (type(img).__module__ != np.__name__) or (type(label).__module__ != np.__name__): + raise ValueError('Error. Ndarray expected. Image: (W, H, D), [0, 1]; Label: (W, H, D, L), [0, 1].') + + if img.ndim != 3: + raise ValueError('Error. 3d tensor expected. Got: ' + str(img.ndim)) + + if label.ndim != 4: + raise ValueError('Error. 4d tensor expected. Got: ' + str(label.ndim)) + + # check image to label shape consistency + if (img.shape[0] != label.shape[0]) or (img.shape[1] != label.shape[1]) or (img.shape[2] != label.shape[2]): + raise ValueError('Error. Image shape and label shape inconsistent: ' + + str(img.shape) + ', ' + str(label.shape)) + + # set image + self.__check_and_set_img(img) + + # set label + self.__check_and_set_label(label) + + # run crf + res_arr = self.__run_crf(label.shape) + + return CRF3DProcessor.__prepare_and_return_result(res_arr, img.shape) + + # # + # + # Private convenience methods + # + # # + def __check_and_set_img(self, img): + if img.ndim != 3: + raise ValueError('Error. 3d tensor expected. Got: ' + str(img.ndim)) + + if np.amin(img) < 0 or np.amax(img) > 1: + raise ValueError('Error. Image must be normalized to [0, 1].') + + # pass image to crf + self.__processor.set_image(np.array(img, dtype=np.float64), + img.shape[0], + img.shape[1], + img.shape[2], + 'float64') + + return img.shape + + def __check_and_set_label(self, label): + if label.ndim != 4: + raise ValueError('Error. 4d tensor expected. Got: ' + str(label.ndim)) + + if np.amin(label) < 0 or np.amax(label) > 1: + raise ValueError('Error. Label must be normalized to [0, 1].') + + # pass labels to crf + self.__processor.set_feature(np.array(label, dtype=np.float64), + label.shape[0], + label.shape[1], + label.shape[2], + label.shape[3], + 'float64') + + return label.shape + + @staticmethod + def __prepare_and_return_result(res_arr, shape): + if len(shape) != 3: + raise ValueError('Error. 3-d shape expected. Got: ' + str(len(shape))) + + if len(res_arr) != shape[0] * shape[1] * shape[2]: + raise ValueError('Error. Length of array inconsisten with length of shape (' + + str(len(res_arr)) + ', ' + str(shape[0] * shape[1] * shape[2]) + ').') + + # Reshape with fortran order + return np.asarray(res_arr, dtype=np.int16).reshape(shape, order='F') + + # # + # + # CRF + # + # # + def __run_crf(self, label_shape): + # check memory + if (not self.__ignore_mem) and (self.memory_percent(label_shape) > 100.0): + raise RuntimeError('Current image exceeds available memory. Free vs Needed: (' + + str(round(self.memory_free()/10**9, 2)) + 'GB, ' + + str(round(self.memory_usage(label_shape)/10**9, 2)) + 'GB).') + + if self.__dynamic_z: + xy_mean = np.mean([label_shape[0], label_shape[1]]) + if (xy_mean != label_shape[2]) and (xy_mean != 0): + bi_z_std = self.__bilateral_z_std * (float(label_shape[2]) / xy_mean) + if self.__verbose: + print('Set Bilateral Z Std to ' + str(bi_z_std) + '.') + self.__processor.set_bi_z_std(bi_z_std) + + return self.__processor.calc_res_array() + + +''' +Implemented by Wendong Xu +''' +class CRF2DProcessor(object): + def __init__(self, + max_iterations=10, + pos_x_std=3.0, + pos_y_std=3.0, + pos_w=3.0, + bilateral_x_std=60.0, + bilateral_y_std=60.0, + bilateral_intensity_std=20.0, + bilateral_w=10.0, + ignore_memory=False, + verbose=False): + self.__verbose = verbose + self.__ignore_mem = ignore_memory + self.__processor = di.DenseCRF2DProcessor(max_iterations, + pos_x_std, + pos_y_std, + pos_w, + bilateral_x_std, + bilateral_y_std, + bilateral_intensity_std, + bilateral_w, + verbose) + + + @staticmethod + def memory_usage(label_shape): + if len(label_shape) != 3: + raise ValueError('Error. 3-dimensional tensor expected. Got: ' + str(len(label_shape))) + return di.DenseCRF2DProcessor.get_memory_requirements(label_shape[0], label_shape[1], label_shape[2]) + + + @staticmethod + def memory_free(): + return psutil.virtual_memory().free + + + @staticmethod + def memory_percent(label_shape): + return float(CRF2DProcessor.memory_usage(label_shape)) / float(CRF2DProcessor.memory_free()) * 100.0 + + + def set_data_and_run(self, img, label): + if (type(img).__module__ != np.__name__) or (type(label).__module__ != np.__name__): + raise ValueError('Error. Ndarray expected. Image: (W, H), [0, 1]; Label: (W, H, L), [0, 1].') + + if img.ndim != 2: + raise ValueError('Error. 2d tensor expected. Got: ' + str(img.ndim)) + + if label.ndim != 3: + raise ValueError('Error. 3d tensor expected. Got: ' + str(label.ndim)) + + # check image to label shape consistency + if (img.shape[0] != label.shape[0]) or (img.shape[1] != label.shape[1]) or (img.shape[2] != label.shape[2]): + raise ValueError('Error. Image shape and label shape inconsistent: ' + + str(img.shape) + ', ' + str(label.shape)) + + self.__check_and_set_img(img) + self.__check_and_set_label(label) + res_arr = self.__run_crf(label.shape) + return CRF2DProcessor.__prepare_and_return_result(res_arr, img.shape) + + + def __check_and_set_img(self, img): + if img.ndim != 2: + raise ValueError('Error. 2d tensor expected. Got: ' + str(img.ndim)) + if np.amin(img) < 0 or np.amax(img) > 1: + raise ValueError('Error. Image must be normalized to [0, 1].') + self.__processor.set_image(np.array(img, dtype=np.float64), + img.shape[0], + img.shape[1], + 'float64') + return img.shape + + + def __check_and_set_label(self, label): + if label.ndim != 3: + raise ValueError('Error. 3d tensor expected. Got: ' + str(label.ndim)) + + if np.amin(label) < 0 or np.amax(label) > 1: + raise ValueError('Error. Label must be normalized to [0, 1].') + + # pass labels to crf + self.__processor.set_feature(np.array(label, dtype=np.float64), + label.shape[0], + label.shape[1], + label.shape[2], + 'float64') + return label.shape + + + @staticmethod + def __prepare_and_return_result(res_arr, shape): + if len(shape) != 2: + raise ValueError('Error. 2-d shape expected. Got: ' + str(len(shape))) + + if len(res_arr) != shape[0] * shape[1]: + raise ValueError('Error. Length of array inconsisten with length of shape (' + + str(len(res_arr)) + ', ' + str(shape[0] * shape[1]) + ').') + # Reshape with fortran order + return np.asarray(res_arr, dtype=np.int16).reshape(shape, order='F') + + + def __run_crf(self, label_shape): + # check memory + if (not self.__ignore_mem) and (self.memory_percent(label_shape) > 100.0): + raise RuntimeError('Current image exceeds available memory. Free vs Needed: (' + + str(round(self.memory_free()/10**9, 2)) + 'GB, ' + + str(round(self.memory_usage(label_shape)/10**9, 2)) + 'GB).') + return self.__processor.calc_res_array() if __name__ == '__main__': - print('All modules loaded successfully.') + print('All modules loaded successfully.') diff --git a/denseinference/lib/refine_3d/dense_inference.cpp b/denseinference/lib/refine_3d/dense_inference.cpp index 6fbb0a4..5043666 100644 --- a/denseinference/lib/refine_3d/dense_inference.cpp +++ b/denseinference/lib/refine_3d/dense_inference.cpp @@ -1,11 +1,15 @@ /* - * - * Python wrapper for DenseCRF - * - * by Marc Bickel (2015) - * - */ - +* +* Python wrapper for DenseCRF +* +* by Marc Bickel (2015) +* +* Python wrapper for 2D-denseCRF +* +* by Wendong Xu (2018) +* +*/ + #include #include @@ -30,271 +34,273 @@ using namespace boost::python; ///// Initialize / Dealloc ///// ////////////////////////////////// DenseCRF3DProcessor::DenseCRF3DProcessor() { - feat_row_ = 0, feat_col_ = 0, feat_slc_ = 0, feat_channel_ = 0; - img_row_ = 0, img_col_ = 0, img_slc_ = 0; - - inp_.MaxIterations = 10; - inp_.PosXStd = 3; - inp_.PosYStd = 3; - inp_.PosZStd = 3; - inp_.PosW = 3; - inp_.BilateralXStd = 60; - inp_.BilateralYStd = 60; - inp_.BilateralZStd = 60; - inp_.BilateralGStd = 20; - inp_.BilateralW = 10; - - verb_ = false; + feat_row_ = 0, feat_col_ = 0, feat_slc_ = 0, feat_channel_ = 0; + img_row_ = 0, img_col_ = 0, img_slc_ = 0; + + inp_.MaxIterations = 10; + inp_.PosXStd = 3; + inp_.PosYStd = 3; + inp_.PosZStd = 3; + inp_.PosW = 3; + inp_.BilateralXStd = 60; + inp_.BilateralYStd = 60; + inp_.BilateralZStd = 60; + inp_.BilateralGStd = 20; + inp_.BilateralW = 10; + + verb_ = false; } DenseCRF3DProcessor::DenseCRF3DProcessor(int MaxIterations, float PosXStd, float PosYStd, float PosZStd, float PosW, float BilateralXStd, float BilateralYStd, float BilateralZStd, float BilateralGStd, float BilateralW, bool verbose) { - feat_row_ = 0, feat_col_ = 0, feat_slc_ = 0, feat_channel_ = 0; - img_row_ = 0, img_col_ = 0, img_slc_ = 0; - - inp_.MaxIterations = MaxIterations; - inp_.PosXStd = PosXStd; - inp_.PosYStd = PosYStd; - inp_.PosZStd = PosZStd; - inp_.PosW = PosW; - inp_.BilateralXStd = BilateralXStd; - inp_.BilateralYStd = BilateralYStd; - inp_.BilateralZStd = BilateralZStd; - inp_.BilateralGStd = BilateralGStd; - inp_.BilateralW = BilateralW; - - verb_ = verbose; + feat_row_ = 0, feat_col_ = 0, feat_slc_ = 0, feat_channel_ = 0; + img_row_ = 0, img_col_ = 0, img_slc_ = 0; + + inp_.MaxIterations = MaxIterations; + inp_.PosXStd = PosXStd; + inp_.PosYStd = PosYStd; + inp_.PosZStd = PosZStd; + inp_.PosW = PosW; + inp_.BilateralXStd = BilateralXStd; + inp_.BilateralYStd = BilateralYStd; + inp_.BilateralZStd = BilateralZStd; + inp_.BilateralGStd = BilateralGStd; + inp_.BilateralW = BilateralW; + + verb_ = verbose; } void DenseCRF3DProcessor::tidy_up(void) { - if(unary_3d_ != NULL) { - delete[] unary_3d_; - unary_3d_ = NULL; - } - - if(feat_ != NULL) { - delete[] feat_; - feat_ = NULL; - } - - if(img_ != NULL) { - delete[] img_; - img_ = NULL; - } - - if(map_ != NULL) { - delete[] map_; - map_ = NULL; - } - - // set shape variables to 0 - feat_row_ = 0, feat_col_ = 0, feat_slc_ = 0, feat_channel_ = 0; - img_row_ = 0, img_col_ = 0, img_slc_ = 0; + if (unary_3d_ != NULL) { + delete[] unary_3d_; + unary_3d_ = NULL; + } + + if (feat_ != NULL) { + delete[] feat_; + feat_ = NULL; + } + + if (img_ != NULL) { + delete[] img_; + img_ = NULL; + } + + if (map_ != NULL) { + delete[] map_; + map_ = NULL; + } + + // set shape variables to 0 + feat_row_ = 0, feat_col_ = 0, feat_slc_ = 0, feat_channel_ = 0; + img_row_ = 0, img_col_ = 0, img_slc_ = 0; } DenseCRF3DProcessor::~DenseCRF3DProcessor() { - tidy_up(); + tidy_up(); } //////////////////////////////////// ///// Protected Helper Methods ///// //////////////////////////////////// PyObject* DenseCRF3DProcessor::stdVecToPyListFloat(const std::vector& vec) { - boost::python::list* l = new boost::python::list(); - for(size_t i = 0; i < vec.size(); i++) - (*l).append(vec[i]); + boost::python::list* l = new boost::python::list(); + for (size_t i = 0; i < vec.size(); i++) + (*l).append(vec[i]); - return l->ptr(); + return l->ptr(); } PyObject* DenseCRF3DProcessor::stdVecToPyListShort(const std::vector& vec) { - boost::python::list* l = new boost::python::list(); - for(size_t i = 0; i < vec.size(); i++) - (*l).append(vec[i]); + boost::python::list* l = new boost::python::list(); + for (size_t i = 0; i < vec.size(); i++) + (*l).append(vec[i]); - return l->ptr(); + return l->ptr(); } -void DenseCRF3DProcessor::calc_unary( float * unary, float * feat, int feat_row, int feat_col, int feat_slc, int feat_channel ) { - // calc a -log() over all entries - for (long i = 0; i < (feat_row * feat_col * feat_slc * feat_channel); i++) { - unary[i] = -log(feat[i]); - } +void DenseCRF3DProcessor::calc_unary(float * unary, float * feat, int feat_row, int feat_col, int feat_slc, int feat_channel) { + // calc a -log() over all entries + for (long i = 0; i < (feat_row * feat_col * feat_slc * feat_channel); i++) { + unary[i] = -log(feat[i]); + } } //////////////////////////////////////////////////////////////// ///// Image and Feature Volume (and BilateralZStd) Setters ///// //////////////////////////////////////////////////////////////// -void DenseCRF3DProcessor::set_feature( boost::python::numeric::array feature, int feat_row, int feat_col, int feat_slc, int feat_channel, const std::string &numpy_data_type ) { - // first tidy up - if (feat_ != NULL) - delete [] feat_; - - if(numpy_data_type == "float64") { - boost::python::numeric::array in = feature; - - // set shape - feat_row_ = feat_row; - feat_col_ = feat_col; - feat_slc_ = feat_slc; - feat_channel_ = feat_channel; - - // allocate memory - feat_ = new float[feat_row * feat_col * feat_slc * feat_channel]; - - // iterate over numpy array, extract and cast to float - for( int k=0; k(in[make_tuple(i, j, k, c)]); - } - } else { - throw std::invalid_argument("Invalid data type: " + numpy_data_type + ". Abort."); - } +void DenseCRF3DProcessor::set_feature(boost::python::numeric::array feature, int feat_row, int feat_col, int feat_slc, int feat_channel, const std::string &numpy_data_type) { + // first tidy up + if (feat_ != NULL) + delete[] feat_; + + if (numpy_data_type == "float64") { + boost::python::numeric::array in = feature; + + // set shape + feat_row_ = feat_row; + feat_col_ = feat_col; + feat_slc_ = feat_slc; + feat_channel_ = feat_channel; + + // allocate memory + feat_ = new float[feat_row * feat_col * feat_slc * feat_channel]; + + // iterate over numpy array, extract and cast to float + for (int k = 0; k(in[make_tuple(i, j, k, c)]); + } + } + else { + throw std::invalid_argument("Invalid data type: " + numpy_data_type + ". Abort."); + } } PyObject* DenseCRF3DProcessor::get_feature() { - // transform to vector - std::vector v; - float * start = feat_; - float * end = feat_ + feat_row_ * feat_col_ * feat_slc_ * feat_channel_; - v.clear(); - v.insert(v.end(), start, end); - - return stdVecToPyListFloat(v); + // transform to vector + std::vector v; + float * start = feat_; + float * end = feat_ + feat_row_ * feat_col_ * feat_slc_ * feat_channel_; + v.clear(); + v.insert(v.end(), start, end); + + return stdVecToPyListFloat(v); } // image packing order is std fortran order, like numpy.ravel(img, order='F'), expects image with dynamic range [0, 1] -void DenseCRF3DProcessor::set_image( boost::python::numeric::array image, int img_row, int img_col, int img_slc, const std::string &numpy_data_type ) { - // first tidy up - if (img_ != NULL) - delete [] img_; - - if (numpy_data_type == "float64") { - boost::python::numeric::array in = image; - - // set shape - img_row_ = img_row; - img_col_ = img_col; - img_slc_ = img_slc; - - // allocate memory - img_ = new float[img_row * img_col * img_slc]; - - // iterate over numpy array, extract and cast to float - for( int k=0; k(in[make_tuple(i, j, k)])) * 255.0f; - } - } else { - throw std::invalid_argument("Invalid data type: " + numpy_data_type + ". Abort."); - } +void DenseCRF3DProcessor::set_image(boost::python::numeric::array image, int img_row, int img_col, int img_slc, const std::string &numpy_data_type) { + // first tidy up + if (img_ != NULL) + delete[] img_; + + if (numpy_data_type == "float64") { + boost::python::numeric::array in = image; + + // set shape + img_row_ = img_row; + img_col_ = img_col; + img_slc_ = img_slc; + + // allocate memory + img_ = new float[img_row * img_col * img_slc]; + + // iterate over numpy array, extract and cast to float + for (int k = 0; k(in[make_tuple(i, j, k)])) * 255.0f; + } + } + else { + throw std::invalid_argument("Invalid data type: " + numpy_data_type + ". Abort."); + } } PyObject* DenseCRF3DProcessor::get_image() { - // transform to vector - std::vector v; - float * start = img_; - float * end = img_ + img_row_ * img_col_ * img_slc_; - v.clear(); - v.insert(v.end(), start, end); - - return stdVecToPyListFloat(v); + // transform to vector + std::vector v; + float * start = img_; + float * end = img_ + img_row_ * img_col_ * img_slc_; + v.clear(); + v.insert(v.end(), start, end); + + return stdVecToPyListFloat(v); } void DenseCRF3DProcessor::set_pos_z_std(float PosZStd) { - inp_.PosZStd = PosZStd; + inp_.PosZStd = PosZStd; } void DenseCRF3DProcessor::set_bi_z_std(float BilateralZStd) { - inp_.BilateralZStd = BilateralZStd; + inp_.BilateralZStd = BilateralZStd; } ///////////////////////// ///// CRF Processor ///// ///////////////////////// void DenseCRF3DProcessor::calculate_map(std::vector &v_map) { - // check if image and features are set - if (img_ == NULL || feat_ == NULL) { - std::cerr << "Error. Image or features not set. Abort." << std::endl; - return; - } - - // check if shapes of image and features are equal - if ((img_row_ != feat_row_) || (img_col_ != feat_col_) || (img_slc_ != feat_slc_)) { - std::cerr << "Error. Image shape differs from features shape. Abort." << std::endl; - return; - } - - // check if shape not 0 - if ((!img_row_) || (!img_col_) || (!img_slc_)) { - std::cerr << "Error. Either rows, columns or slices are 0. Abort." << std::endl; - return; - } - // transform feature to unary - unary_3d_ = new float[feat_row_ * feat_col_ * feat_slc_ * feat_channel_]; - calc_unary(unary_3d_, feat_, feat_row_, feat_col_, feat_slc_, feat_channel_); - - // Initialize the magic.. - DenseCRF3D crf(img_slc_, img_col_, img_row_, feat_channel_); - - // Set verbose - crf.verbose(verb_); - - // Specify the unary potential as an array of size W*H*D*(#classes) - // packing order: x0y0z0l0 x0y0z0l1 x0y0z0l2 .. x1y0z0l0 x1y0z0l1 ... (row-order) (image packing order is std fortran order, like numpy.ravel(img, order='F')) - if (crf.isVerbose()) - std::cout << "pre energy" << std::endl; - crf.setUnaryEnergy(unary_3d_); - // add a gray intensity independent term (feature = pixel location 0..W-1, 0..H-1) - if (crf.isVerbose()) - std::cout << "pre gaussian" << std::endl; - crf.addPairwiseGaussian(inp_.PosXStd, inp_.PosYStd, inp_.PosZStd, inp_.PosW); - // add a gray intensity dependent term (feature = xyzg) - if (crf.isVerbose()) - std::cout << "pre bilateral" << std::endl; - crf.addPairwiseBilateral(inp_.BilateralXStd, inp_.BilateralYStd, inp_.BilateralZStd, inp_.BilateralGStd, img_, inp_.BilateralW); - - // create map (=maximum a posteriori) - map_ = new short[img_row_ * img_col_ * img_slc_]; - - if (crf.isVerbose()) - std::cout << "pre premap" << std::endl; - - crf.map(inp_.MaxIterations, map_); - - // print out some values from raw Map array configuration - if (crf.isVerbose()) { - long ctr = 0; - for (long i = 0; i < img_row_ * img_col_ * img_slc_; i++) { - if (map_[i] > 0) { - ctr++; - } - if (i % (img_row_ * img_col_ * img_slc_/20) == 0) { - std::cout << map_[i] << " "; - } - } - std::cout << std::endl << "Maximum a posteriori > 0: " << ctr << std::endl; - } - - // transform to vector - short * start = map_; - short * end = map_ + img_row_ * img_col_ * img_slc_; - v_map.clear(); - v_map.insert(v_map.end(), start, end); - - // tidy up - tidy_up(); - - if (crf.isVerbose()) - std::cout << "Done." << std::endl; + // check if image and features are set + if (img_ == NULL || feat_ == NULL) { + std::cerr << "Error. Image or features not set. Abort." << std::endl; + return; + } + + // check if shapes of image and features are equal + if ((img_row_ != feat_row_) || (img_col_ != feat_col_) || (img_slc_ != feat_slc_)) { + std::cerr << "Error. Image shape differs from features shape. Abort." << std::endl; + return; + } + + // check if shape not 0 + if ((!img_row_) || (!img_col_) || (!img_slc_)) { + std::cerr << "Error. Either rows, columns or slices are 0. Abort." << std::endl; + return; + } + // transform feature to unary + unary_3d_ = new float[feat_row_ * feat_col_ * feat_slc_ * feat_channel_]; + calc_unary(unary_3d_, feat_, feat_row_, feat_col_, feat_slc_, feat_channel_); + + // Initialize the magic.. + DenseCRF3D crf(img_slc_, img_col_, img_row_, feat_channel_); + + // Set verbose + crf.verbose(verb_); + + // Specify the unary potential as an array of size W*H*D*(#classes) + // packing order: x0y0z0l0 x0y0z0l1 x0y0z0l2 .. x1y0z0l0 x1y0z0l1 ... (row-order) (image packing order is std fortran order, like numpy.ravel(img, order='F')) + if (crf.isVerbose()) + std::cout << "pre energy" << std::endl; + crf.setUnaryEnergy(unary_3d_); + // add a gray intensity independent term (feature = pixel location 0..W-1, 0..H-1) + if (crf.isVerbose()) + std::cout << "pre gaussian" << std::endl; + crf.addPairwiseGaussian(inp_.PosXStd, inp_.PosYStd, inp_.PosZStd, inp_.PosW); + // add a gray intensity dependent term (feature = xyzg) + if (crf.isVerbose()) + std::cout << "pre bilateral" << std::endl; + crf.addPairwiseBilateral(inp_.BilateralXStd, inp_.BilateralYStd, inp_.BilateralZStd, inp_.BilateralGStd, img_, inp_.BilateralW); + + // create map (=maximum a posteriori) + map_ = new short[img_row_ * img_col_ * img_slc_]; + + if (crf.isVerbose()) + std::cout << "pre premap" << std::endl; + + crf.map(inp_.MaxIterations, map_); + + // print out some values from raw Map array configuration + if (crf.isVerbose()) { + long ctr = 0; + for (long i = 0; i < img_row_ * img_col_ * img_slc_; i++) { + if (map_[i] > 0) { + ctr++; + } + if (i % (img_row_ * img_col_ * img_slc_ / 20) == 0) { + std::cout << map_[i] << " "; + } + } + std::cout << std::endl << "Maximum a posteriori > 0: " << ctr << std::endl; + } + + // transform to vector + short * start = map_; + short * end = map_ + img_row_ * img_col_ * img_slc_; + v_map.clear(); + v_map.insert(v_map.end(), start, end); + + // tidy up + tidy_up(); + + if (crf.isVerbose()) + std::cout << "Done." << std::endl; } ///////////////////////// ///// Public Helper ///// ///////////////////////// long DenseCRF3DProcessor::get_memory_requirements(const int W, const int H, const int D, const int C) { - long M = (long)C; - long N = (long)(W*H*D); + long M = (long)C; + long N = (long)(W*H*D); - long FL = (long)sizeof(float); - long SH = (long)sizeof(short); + long FL = (long)sizeof(float); + long SH = (long)sizeof(short); - long dCRF = 6 * N * M * FL; - long d3CRF = (2 * N * M + N) * FL + (N * M * SH); - long bil = 6 * N; + long dCRF = 6 * N * M * FL; + long d3CRF = (2 * N * M + N) * FL + (N * M * SH); + long bil = 6 * N; - return dCRF + d3CRF + bil; + return dCRF + d3CRF + bil; } @@ -302,31 +308,275 @@ long DenseCRF3DProcessor::get_memory_requirements(const int W, const int H, cons ///// Python Wrapper ///// ////////////////////////// PyObject* DenseCRF3DProcessor::calc_res_array() { - std::vector map; - calculate_map(map); - - // transform to numpy array and return it - if (map.size() > 0) { - return stdVecToPyListShort(map); - } else { - throw std::runtime_error("Invalid map. Calculate_map() did not return successfully."); - } + std::vector map; + calculate_map(map); + + // transform to numpy array and return it + if (map.size() > 0) { + return stdVecToPyListShort(map); + } + else { + throw std::runtime_error("Invalid map. Calculate_map() did not return successfully."); + } +} + + + +//////////////////////////////////////// +///// 2D-denseCRF python wrapper ///// +///// Implemented by Wendong Xu ///// +//////////////////////////////////////// +DenseCRF2DProcessor::DenseCRF2DProcessor() { + feat_row_ = 0, feat_col_ = 0, feat_channel_ = 0; + img_row_ = 0, img_col_ = 0; + + inp_.MaxIterations = 10; + inp_.PosXStd = 3; + inp_.PosYStd = 3; + inp_.PosW = 3; + inp_.BilateralXStd = 60; + inp_.BilateralYStd = 60; + inp_.BilateralGStd = 20; + inp_.BilateralW = 10; + + verb_ = false; +} + +DenseCRF2DProcessor::DenseCRF2DProcessor(int MaxIterations, float PosXStd, float PosYStd, float PosW, float BilateralXStd, float BilateralYStd, float BilateralGStd, float BilateralW, bool verbose) { + feat_row_ = 0, feat_col_ = 0, feat_channel_ = 0; + img_row_ = 0, img_col_ = 0; + + inp_.MaxIterations = MaxIterations; + inp_.PosXStd = PosXStd; + inp_.PosYStd = PosYStd; + inp_.PosW = PosW; + inp_.BilateralXStd = BilateralXStd; + inp_.BilateralYStd = BilateralYStd; + inp_.BilateralGStd = BilateralGStd; + inp_.BilateralW = BilateralW; + + verb_ = verbose; +} + +void DenseCRF2DProcessor::tidy_up(void) { + if (unary_2d_ != NULL) { + delete[] unary_2d_; + unary_2d_ = NULL; + } + + if (feat_ != NULL) { + delete[] feat_; + feat_ = NULL; + } + + if (img_ != NULL) { + delete[] img_; + img_ = NULL; + } + + if (map_ != NULL) { + delete[] map_; + map_ = NULL; + } + + // set shape variables to 0 + feat_row_ = 0, feat_col_ = 0, feat_channel_ = 0; + img_row_ = 0, img_col_ = 0; +} + +DenseCRF2DProcessor::~DenseCRF2DProcessor() { + tidy_up(); +} + +PyObject* DenseCRF2DProcessor::stdVecToPyListShort(const std::vector& vec) { + boost::python::list* l = new boost::python::list(); + for (size_t i = 0; i < vec.size(); i++) + (*l).append(vec[i]); + + return l->ptr(); +} + +PyObject* DenseCRF2DProcessor::stdVecToPyListFloat(const std::vector& vec) { + boost::python::list* l = new boost::python::list(); + for (size_t i = 0; i < vec.size(); i++) + (*l).append(vec[i]); + + return l->ptr(); +} + +void DenseCRF2DProcessor::calc_unary(float * unary, float * feat, int feat_row, int feat_col, int feat_channel) { + for (long i = 0; i < (feat_row * feat_col * feat_channel); i++) { + unary[i] = -log(feat[i]); + } +} + +void DenseCRF2DProcessor::set_feature(boost::python::numeric::array feature, int feat_row, int feat_col, int feat_channel, const std::string &numpy_data_type) { + if (feat_ != NULL) { + delete[]feat_; + } + if (numpy_data_type == "float64") { + boost::python::numeric::array in = feature; + feat_row_ = feat_row; + feat_col_ = feat_col; + feat_channel_ = feat_channel; + feat_ = new float[feat_row * feat_col * feat_channel]; + for (int j = 0; j < feat_col; j++) { + for (int i = 0; i < feat_row; i++) { + for (int c = 0; c < feat_channel; c++) { + feat_[(feat_col*feat_row + j * feat_row + i)*feat_channel + c] = (float)extract(in[make_tuple(i, j, c)]); + } + } + } + } + else { + throw std::invalid_argument("Invalid data type: " + numpy_data_type + ". Abort."); + } +} + +PyObject* DenseCRF2DProcessor::get_feature() { + std::vector v; + float * start = feat_; + float * end = feat_ + feat_row_ * feat_col_ * feat_channel_; + v.clear(); + v.insert(v.end(), start, end); + return stdVecToPyListFloat(v); +} + +void DenseCRF2DProcessor::set_image(boost::python::numeric::array image, int img_row, int img_col, const std::string &numpy_data_type) { + if (img_ != NULL) { + delete[] img_; + } + if (numpy_data_type == "float64") { + boost::python::numeric::array in = image; + img_row_ = img_row; + img_col_ = img_col; + img_ = new float[img_row * img_col]; + for (int j = 0; j < img_col; j++) { + for (int i = 0; i < img_row; i++) { + img_[img_row + j * img_row + i] = ((float)extract(in[make_tuple(i, j)])) * 255.0f; + } + } + } + else { + throw std::invalid_argument("Invalid data type: " + numpy_data_type + ". Abort."); + } +} + +PyObject* DenseCRF2DProcessor::get_image() { + std::vector v; + float * start = img_; + float * end = img_ + img_row_ * img_col_; + v.clear(); + v.insert(v.end(), start, end); + + return stdVecToPyListFloat(v); +} + +void DenseCRF2DProcessor::calculate_map(std::vector &v_map) { + if (img_ == NULL || feat_ == NULL) { + std::cerr << "Error. Image or features not set. Abort." << std::endl; + return; + } + if ((img_row_ != feat_row_) || (img_col_ != feat_col_)) { + std::cerr << "Error. Image shape differs from features shape. Abort." << std::endl; + return; + } + if ((!img_row_) || (!img_col_)) { + std::cerr << "Error. Either rows, columns or slices are 0. Abort." << std::endl; + return; + } + unary_2d_ = new float[feat_row_ * feat_col_ * feat_channel_]; + calc_unary(unary_2d_, feat_, feat_row_, feat_col_, feat_channel_); + DenseCRF2D crf(img_col_, img_row_, feat_channel_); + crf.verbose(verb_); + if (crf.isVerbose()) { + std::cout << "pre energy" << std::endl; + } + crf.setUnaryEnergy(unary_2d_); + if (crf.isVerbose()) { + std::cout << "pre gaussian" << std::endl; + } + crf.addPairwiseGaussian(inp_.PosXStd, inp_.PosYStd, inp_.PosW); + if (crf.isVerbose()) { + std::cout << "pre bilateral" << std::endl; + } + map_ = new short[img_row_ * img_col_]; + if (crf.isVerbose()) { + std::cout << "pre premap" << std::endl; + } + crf.map(inp_.MaxIterations, map_); + if (crf.isVerbose()) { + long ctr = 0; + for (long i = 0; i < img_row_ * img_col_; i++) { + if (map_[i] > 0) { + ctr++; + } + if (i % (img_row_ * img_col_ / 20) == 0) { + std::cout << map_[i] << " "; + } + } + std::cout << std::endl << "Maximum a posteriori > 0: " << ctr << std::endl; + } + short * start = map_; + short * end = map_ + img_row_ * img_col_; + v_map.clear(); + v_map.insert(v_map.end(), start, end); + tidy_up(); + if (crf.isVerbose()) { + std::cout << "Done." << std::endl; + } +} + +long DenseCRF2DProcessor::get_memory_requirements(const int W, const int H, const int C) { + long M = (long)C; + long N = (long)(W*H); + + long FL = (long)sizeof(float); + long SH = (long)sizeof(short); + + long dCRF = 6 * N * M * FL; + long d3CRF = (2 * N * M + N) * FL + (N * M * SH); + long bil = 6 * N; + + return dCRF + d3CRF + bil; +} + +PyObject* DenseCRF2DProcessor::calc_res_array() { + std::vector map; + calculate_map(map); + if (map.size() > 0) { + return stdVecToPyListShort(map); + } + else { + throw std::runtime_error("Invalid map. Calculate_map() did not return successfully."); + } } // Expose classes and methods to Python BOOST_PYTHON_MODULE(dense_inference) { - boost::python::numeric::array::set_module_and_type("numpy", "ndarray"); - - class_("DenseCRF3DProcessor", init<>()) - .def(init()) - .def("set_pos_z_std", &DenseCRF3DProcessor::set_pos_z_std) - .def("set_bi_z_std", &DenseCRF3DProcessor::set_bi_z_std) - .def("set_feature", &DenseCRF3DProcessor::set_feature) - .def("get_feature", &DenseCRF3DProcessor::get_feature) - .def("set_image", &DenseCRF3DProcessor::set_image) - .def("get_image", &DenseCRF3DProcessor::get_image) - .def("get_memory_requirements", &DenseCRF3DProcessor::get_memory_requirements) - .staticmethod("get_memory_requirements") - .def("calc_res_array", &DenseCRF3DProcessor::calc_res_array) - ; + boost::python::numeric::array::set_module_and_type("numpy", "ndarray"); + + class_("DenseCRF3DProcessor", init<>()) + .def(init()) + .def("set_pos_z_std", &DenseCRF3DProcessor::set_pos_z_std) + .def("set_bi_z_std", &DenseCRF3DProcessor::set_bi_z_std) + .def("set_feature", &DenseCRF3DProcessor::set_feature) + .def("get_feature", &DenseCRF3DProcessor::get_feature) + .def("set_image", &DenseCRF3DProcessor::set_image) + .def("get_image", &DenseCRF3DProcessor::get_image) + .def("get_memory_requirements", &DenseCRF3DProcessor::get_memory_requirements) + .staticmethod("get_memory_requirements") + .def("calc_res_array", &DenseCRF3DProcessor::calc_res_array) + ; + + class_("DenseCRF2DProcessor", init<>()) + .def(init()) + .def("set_feature", &DenseCRF2DProcessor::set_feature) + .def("get_feature", &DenseCRF2DProcessor::get_feature) + .def("set_image", &DenseCRF2DProcessor::set_image) + .def("get_image", &DenseCRF2DProcessor::get_image) + .def("get_memory_requirements", &DenseCRF2DProcessor::get_memory_requirements) + .staticmethod("get_memory_requirements") + .def("calc_res_array", &DenseCRF2DProcessor::calc_res_array) + ; } diff --git a/denseinference/lib/refine_3d/dense_inference.h b/denseinference/lib/refine_3d/dense_inference.h index 0c2a84c..07831a4 100644 --- a/denseinference/lib/refine_3d/dense_inference.h +++ b/denseinference/lib/refine_3d/dense_inference.h @@ -9,18 +9,18 @@ #include struct InputData3D { - int MaxIterations; - float PosXStd; - float PosYStd; - float PosZStd; - float PosW; - float BilateralXStd; - float BilateralYStd; - float BilateralZStd; - float BilateralGStd; // gray - float BilateralW; - - void print_settings(void); + int MaxIterations; + float PosXStd; + float PosYStd; + float PosZStd; + float PosW; + float BilateralXStd; + float BilateralYStd; + float BilateralZStd; + float BilateralGStd; // gray + float BilateralW; + + void print_settings(void); }; class DenseCRF3DProcessor { @@ -78,6 +78,60 @@ class DenseCRF3DProcessor { PyObject* calc_res_array(void); }; + +struct InputData2D { + int MaxIterations; + float PosXStd; + float PosYStd; + float PosW; + float BilateralXStd; + float BilateralYStd; + float BilateralGStd; // gray + float BilateralW; + + void print_settings(void); +}; + +class DenseCRF2DProcessor { +protected: + InputData2D inp_; + + bool verb_ = false; + + float *unary_2d_ = NULL; + float *feat_ = NULL; + float *img_ = NULL; + short *map_ = NULL; + + int feat_row_, feat_col_, feat_channel_; + int img_row_, img_col_; + + PyObject* stdVecToPyListShort(const std::vector& vec); + PyObject* stdVecToPyListFloat(const std::vector& vec); + + void calc_unary(float * unary, float * feat, int feat_row, int feat_col, int feat_channel); + + void tidy_up(void); + +public: + DenseCRF2DProcessor(); + DenseCRF2DProcessor(int MaxIterations, float PosXStd, float PosYStd, float PosW, float BilateralXStd, float BilateralYStd, float BilateralGStd, float BilateralW, bool verbose); + ~DenseCRF2DProcessor(); + + void set_feature(boost::python::numeric::array feature, int feat_row, int feat_col, int feat_channel, const std::string &numpy_data_type); + PyObject* get_feature(); + + void set_image(boost::python::numeric::array image, int img_row, int img_col, const std::string &numpy_data_type); + PyObject* get_image(); + + void calculate_map(std::vector &v_map); + + static long get_memory_requirements(const int W, const int H, const int C); + + PyObject* calc_res_array(void); + +}; + #endif diff --git a/readme.md b/readme.md index c01b4db..8e2a7fa 100644 --- a/readme.md +++ b/readme.md @@ -1,3 +1,6 @@ +## 2D-denseCRF python wrapper +Here's a 2D-denseCRF wrapper. + # Dense Inference Wrapper ## Overview From 0fc796b4c84aeed08835bae41f1f04b9d22a5516 Mon Sep 17 00:00:00 2001 From: kirai Date: Tue, 8 May 2018 15:52:30 +0800 Subject: [PATCH 3/5] update readme --- readme.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 8e2a7fa..27babab 100644 --- a/readme.md +++ b/readme.md @@ -1,5 +1,6 @@ ## 2D-denseCRF python wrapper -Here's a 2D-denseCRF wrapper. +Here's a 2D-denseCRF wrapper written by Wendong Xu +Any question please connect kirai [dot] wendong [at] gmail [dot] com # Dense Inference Wrapper From 879b1c9559a6b50ee81ca2bdfa87e42051cab0c9 Mon Sep 17 00:00:00 2001 From: kirai Date: Tue, 8 May 2018 19:30:24 +0800 Subject: [PATCH 4/5] fixed: corrupted double-linked list --- denseinference/CRFProcessor.py | 2 +- denseinference/lib/refine_3d/dense_inference.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/denseinference/CRFProcessor.py b/denseinference/CRFProcessor.py index b52bda0..401cc90 100644 --- a/denseinference/CRFProcessor.py +++ b/denseinference/CRFProcessor.py @@ -279,7 +279,7 @@ def set_data_and_run(self, img, label): raise ValueError('Error. 3d tensor expected. Got: ' + str(label.ndim)) # check image to label shape consistency - if (img.shape[0] != label.shape[0]) or (img.shape[1] != label.shape[1]) or (img.shape[2] != label.shape[2]): + if (img.shape[0] != label.shape[0]) or (img.shape[1] != label.shape[1]): raise ValueError('Error. Image shape and label shape inconsistent: ' + str(img.shape) + ', ' + str(label.shape)) diff --git a/denseinference/lib/refine_3d/dense_inference.cpp b/denseinference/lib/refine_3d/dense_inference.cpp index 5043666..db59c8f 100644 --- a/denseinference/lib/refine_3d/dense_inference.cpp +++ b/denseinference/lib/refine_3d/dense_inference.cpp @@ -423,7 +423,7 @@ void DenseCRF2DProcessor::set_feature(boost::python::numeric::array feature, int for (int j = 0; j < feat_col; j++) { for (int i = 0; i < feat_row; i++) { for (int c = 0; c < feat_channel; c++) { - feat_[(feat_col*feat_row + j * feat_row + i)*feat_channel + c] = (float)extract(in[make_tuple(i, j, c)]); + feat_[(j * feat_row + i) * feat_channel + c] = (float)extract(in[make_tuple(i, j, c)]); } } } @@ -453,7 +453,7 @@ void DenseCRF2DProcessor::set_image(boost::python::numeric::array image, int img img_ = new float[img_row * img_col]; for (int j = 0; j < img_col; j++) { for (int i = 0; i < img_row; i++) { - img_[img_row + j * img_row + i] = ((float)extract(in[make_tuple(i, j)])) * 255.0f; + img_[j * img_row + i] = ((float)extract(in[make_tuple(i, j)])) * 255.0f; } } } From 45c70988e8c9e4836cd28cb6e710de89168835fb Mon Sep 17 00:00:00 2001 From: kirai Date: Wed, 9 May 2018 10:11:40 +0800 Subject: [PATCH 5/5] add a demo for comparing 2DdenceCRF and 3DdenseCRF --- demo.py | 141 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 demo.py diff --git a/demo.py b/demo.py new file mode 100644 index 0000000..c45157b --- /dev/null +++ b/demo.py @@ -0,0 +1,141 @@ +__author__ = 'Wendong Xu' + +from denseinference import CRFProcessor +import nibabel as nib +import numpy as np + + +slice_shape = (512, 512) +crf2d = CRFProcessor.CRF2DProcessor() +crf3d = CRFProcessor.CRF3DProcessor() + + +def histeq_processor(img): + """Histogram equalization""" + nbr_bins=256 + #get image histogram + imhist,bins = np.histogram(img.flatten(),nbr_bins,normed=True) + cdf = imhist.cumsum() #cumulative distribution function + cdf = 255 * cdf / cdf[-1] #normalize + #use linear interpolation of cdf to find new pixel values + original_shape = img.shape + img = np.interp(img.flatten(),bins[:-1],cdf) + img=img/255.0 + return img.reshape(original_shape) + + +def norm_hounsfield_dyn(arr, c_min=0.1, c_max=0.3): + """ Converts from hounsfield units to float64 image with range 0.0 to 1.0 """ + # calc min and max + min, max = np.amin(arr), np.amax(arr) + if min <= 0: + arr = np.clip(arr, min * c_min, max * c_max) + # right shift to zero + arr = np.abs(min * c_min) + arr + else: + arr = np.clip(arr, min, max * c_max) + # left shift to zero + arr = arr - min + # normalization + norm_fac = np.amax(arr) + if norm_fac != 0: + norm = np.divide( + np.multiply(arr, 255), + np.amax(arr)) + else: # don't divide through 0 + norm = np.multiply(arr, 255) + + norm = np.clip(np.multiply(norm, 0.00390625), 0, 1) + return norm + + +def process_img_label(imgvol, segvol): + """ + Process a given image volume and its label and return arrays as a new copy + :param imgvol: + :param label_vol: + :return: + """ + imgvol_downscaled = np.zeros((slice_shape[0], slice_shape[1], imgvol.shape[2])) + segvol_downscaled = np.zeros((slice_shape[0], slice_shape[1], imgvol.shape[2])) + imgvol[imgvol > 1200] = 0 + + for i in range(imgvol.shape[2]): + # Get the current slice, normalize and downscale + slice = np.copy(imgvol[:, :, i]) + slice = norm_hounsfield_dyn(slice) + slice = histeq_processor(slice) + imgvol_downscaled[:, :, i] = slice + # downscale the label slice for the crf + segvol_downscaled[:, :, i] = segvol[:,:,i] # to_scale(segvol[:, :, i], slice_shape) + + return [imgvol_downscaled, segvol_downscaled] + + +def crf_2d_processing(imgvol, probvol): + ''' + 2DCRF: set_data_and_run(img, prob) + img: [W, H], [0, 1] + prob: [W, H, L], [0, 1] + :param imgvol: [W, H, D] + :param probvol: [W, H, D, L] + :return: + ''' + global crf2d + assert imgvol.shape == probvol.shape, 'Invalid imgvol and probvol length' + vol_size = imgvol.shape[2] + new_probvol = np.zeros((probvol.shape[0], probvol.shape[1], probvol.shape[3], probvol.shape[2])) + new_probvol = probvol.transpose((0, 1, 3, 2)) + + hard_cls = [] + for i in range(0, vol_size): + print('slice {} finished.'.format(i)) + temp_cls = crf2d.set_data_and_run(imgvol[..., i], new_probvol[..., i]) + hard_cls.append(temp_cls) + return hard_cls + + +def crf_3d_processing(imgvol, probvol): + ''' + 3DCRF: set_data_and_run(imgvol, probvol) + imgvol: [W, H, D], [0, 1] + probvol: [W, H, D, L], [0, 1] + + :param imgvol: [W, H, D] + :param probvol: [W, H, D, L] + :return: + ''' + global crf3d + assert imgvol.shape == probvol.shape, 'Invalid imgvol and probvol length' + return crf3d.set_data_and_run(imgvol, probvol) + + +def get_validate_files(imgvol, segvol, probvol): + ''' + :param imgvol: + :param segvol: + :param probvol: + :return: + ''' + new_probvol = np.zeros((probvol.shape[0], probvol.shape[1], probvol.shape[3], probvol.shape[2])) + new_probvol = probvol.transpose((0, 1, 3, 2)) + new_probvol = new_probvol.astype(np.float64) + new_probvol[new_probvol == 0] = 1e-18 + + imgvol, segvol = process_img_label(imgvol, segvol) + return imgvol, segvol, new_probvol + + +if __name__ == '__main__': + imgvol_path = './data/volume-0.nii' + segvol_path = './data/segmentation-0.nii' + probvol_path = './data/test-logits-0.npz' + + imgvol = nib.load(imgvol_path).get_data() + segvol = nib.load(segvol_path).get_data() + probvol = np.load(probvol_path)['logits'] + + imgvol, segvol, probvol = get_validate_files(imgvol, segvol, probvol) + print(imgvol.shape, segvol.shape, probvol.shape) + seg2dcrf = crf_2d_processing(imgvol, probvol) + seg3dcrf = crf_3d_processing(imgvol, probvol)