diff options
Diffstat (limited to 'silx/math/marchingcubes/mc.hpp')
-rw-r--r-- | silx/math/marchingcubes/mc.hpp | 724 |
1 files changed, 724 insertions, 0 deletions
diff --git a/silx/math/marchingcubes/mc.hpp b/silx/math/marchingcubes/mc.hpp new file mode 100644 index 0000000..82eced9 --- /dev/null +++ b/silx/math/marchingcubes/mc.hpp @@ -0,0 +1,724 @@ +/*########################################################################## +# +# Copyright (c) 2015-2016 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ +#ifndef __mc_HPP__ +#define __mc_HPP__ + +#include <iostream> +#include <cmath> +#include <map> +#include <stdexcept> +#include <vector> +#include <assert.h> + + +extern const int MCTriangleTable[256][16]; +extern const unsigned int MCEdgeIndexToCoordOffsets[12][4]; + +#define DEPTH_IDX 0 +#define HEIGHT_IDX 1 +#define WIDTH_IDX 2 + +/** Class Marching cubes + * + * Implements the marching cube algorithm and provides an API to process + * data image by image. + * + * Dimension convention used is (dim0, dim1, dim2) denoted as + * (depth, height, width) with dim2 (= width) being contiguous in memory. + * + * If data is provided as (depth, height, width), resulting vertices + * and normals will be stored as (z, y, x) + * + * Indices in memory for a single cube: + * + * dim 0 (depth) + * | + * | + * 4 +------+ 5 + * /| /| + * / | / | + * 6 +------+ 7| + * | | | | + * |0 +---|--+ 1 --- dim 2 (width) + * | / | / + * |/ |/ + * 2 +------+ 3 + * / + * / + * dim 1 (height) + */ +template <typename FloatIn, typename FloatOut> +class MarchingCubes { +public: + /** Create a marching cube object. + * + * @param level Level at which to build the isosurface + */ + MarchingCubes(const FloatIn level); + + ~MarchingCubes(); + + /** Process a 3D scalar field + * + * @param data Pointer to the data set + * @param depth The 1st dimension of the data set + * @param height The 2nd dimension of the data set + * @param width The 3rd dimension of the data set + * (tightly packed in memory) + */ + void process(const FloatIn * data, + const unsigned int depth, + const unsigned int height, + const unsigned int width); + + /** Init dimension of slices + * + * @param height Height in pixels of the slices + * @param width Width in pixels of the slices + */ + void set_slice_size(const unsigned int height, + const unsigned int width); + + /** Process a slice (i.e., an image) + * + * The size of the images MUST match height and width provided to + * set_slice_size. + * + * The marching cube process 2 consecutive images at a time. + * A slice provided as next parameter MUST be provided as current + * parameter for the next call. + * Example with 3 images: + * + * float * img1; + * float * img2; + * float * img3; + * ... + * mc = MarchingCubes<float>(100.); + * mc.set_slice_size(10, 10); + * mc.process_slice(img1, img2); + * mc.process_slice(img2, img3); + * mc.finish_process(); + * + * @param slice0 Pointer to the nth slice data + * @param slice1 Pointer to the (n+1)th slice of data + */ + void process_slice(const FloatIn * slice0, + const FloatIn * slice1); + + /** Clear marching cube processing internal cache. */ + void finish_process(); + + /** Reset all internal data and counters. */ + void reset(); + + /** Vertices of the isosurface (x, y, z) */ + std::vector<FloatOut> vertices; + + /** Approximation of normals at the vertices (nx, ny, nz) + * + * Current implementation provides coarse (but fast) normals computation + */ + std::vector<FloatOut> normals; + + /** Triangle indices */ + std::vector<unsigned int> indices; + + unsigned int depth; /**< Number of images currently processed */ + unsigned int height; /**< Images height in pixels */ + unsigned int width; /**< Images width in pixels */ + + /** Sampling of the data (depth, height, width) + * + * Default: 1, 1, 1 + */ + unsigned int sampling[3]; + + FloatIn isolevel; /**< Iso level to use */ + bool invert_normals; /**< True to inverse gradient as normals */ + +private: + + /** Start to build isosurface starting with first slice + * + * Bootstrap cache edge_indices + * + * @param slice The first slice of the data + * @param next The second slice + */ + void first_slice(const FloatIn * slice, + const FloatIn * next); + + /** Process an edge + * + * @param value0 Data at 'begining' of edge + * @param value Data at 'end' of edge + * @param depth Depth coordinate of the edge position + * @param row Row coordinate of the edge + * @param col Column coordinate of the edge + * @param direction Direction of the edge: 0 for x, 1 for y and 2 for z + * @param previous + * @param current + * @param next + */ + void process_edge(const FloatIn value0, + const FloatIn value, + const unsigned int depth, + const unsigned int row, + const unsigned int col, + const unsigned int direction, + const FloatIn * previous, + const FloatIn * current, + const FloatIn * next); + + /** Return the bit mask of cube corners <= the iso-value. + * + * @param slice1 1st slice of the cube to consider + * @param slice2 2nd slice of the cube to consider + * @param row Row of the cube to consider + * @param col Column of the cube to consider + * @return The bit mask of cube corners <= the iso-value + */ + unsigned char get_cell_code(const FloatIn * slice1, + const FloatIn * slice2, + const unsigned int row, + const unsigned int col); + + /** Compute an edge index from position and edge direction. + * + * @param depth Depth of the origin of the edge + * @param row Row of the origin of the edge + * @param col Column of the origin of the edge + * @param direction 0 for x, 1 for y, 2 for z + * @return The (4D) index of the edge + */ + unsigned int edge_index(const unsigned int depth, + const unsigned int row, + const unsigned int col, + const unsigned int direction); + + /** For each dimension, a map from edge index to vertex index + * + * This caches indices for previously processed slice. + * + * Edge index is the linearized position of the edge using size + 1 + * in all dimensions as coordinates plus the direction as 4th coord. + * WARNING: direction 0 for x, 1 for y and 2 for z + */ + std::map<unsigned int, unsigned int> * edge_indices; +}; + + +/* Implementation */ + +template <typename FloatIn, typename FloatOut> +MarchingCubes<FloatIn, FloatOut>::MarchingCubes(const FloatIn level) +{ + this->edge_indices = 0; + this->reset(); + this->height = 0; + this->width = 0; + this->isolevel = level; + this->invert_normals = true; + this->sampling[0] = 1; + this->sampling[1] = 1; + this->sampling[2] = 1; +} + +template <typename FloatIn, typename FloatOut> +MarchingCubes<FloatIn, FloatOut>::~MarchingCubes() +{ +} + +template <typename FloatIn, typename FloatOut> +void +MarchingCubes<FloatIn, FloatOut>::reset() +{ + this->depth = 0; + this->vertices.clear(); + this->normals.clear(); + this->indices.clear(); + if (this->edge_indices != 0) { + delete this->edge_indices; + this->edge_indices = 0; + } +} + +template <typename FloatIn, typename FloatOut> +void +MarchingCubes<FloatIn, FloatOut>::finish_process() +{ + if (this->edge_indices != 0) { + delete this->edge_indices; + this->edge_indices = 0; + } +} + + +template <typename FloatIn, typename FloatOut> +void +MarchingCubes<FloatIn, FloatOut>::process(const FloatIn * data, + const unsigned int depth, + const unsigned int height, + const unsigned int width) +{ + assert(data != NULL); + unsigned int size = height * width * this->sampling[DEPTH_IDX]; + + /* number of slices minus - 1 to process */ + const unsigned int nb_slices = (depth - 1) / this->sampling[DEPTH_IDX]; + + this->reset(); + this->set_slice_size(height, width); + + for (unsigned int index=0; index < nb_slices; index++) { + const FloatIn * slice0 = data + (index * size); + const FloatIn * slice1 = slice0 + size; + + this->process_slice(slice0, slice1); + } + this->finish_process(); + + this->depth = depth; /* Forced as it might be < depth otherwise */ +} + + +template <typename FloatIn, typename FloatOut> +void +MarchingCubes<FloatIn, FloatOut>::set_slice_size(const unsigned int height, + const unsigned int width) +{ + this->reset(); + this->height = height; + this->width = width; +} + + +template <typename FloatIn, typename FloatOut> +void +MarchingCubes<FloatIn, FloatOut>::process_slice(const FloatIn * slice0, + const FloatIn * slice1) +{ + assert(slice0 != NULL); + assert(slice1 != NULL); + unsigned int row, col; + + if (this->edge_indices == 0) { + /* No previously processed slice, bootstrap */ + this->first_slice(slice0, slice1); + } + + /* Keep reference to cache from previous slice */ + std::map<unsigned int, unsigned int> * previous_edge_indices = + this->edge_indices; + + /* Init cache for this slice */ + this->edge_indices = new std::map<unsigned int, unsigned int>(); + + /* Loop over slice to add vertices */ + for (row=0; row < this->height; row += this->sampling[HEIGHT_IDX]) { + unsigned int line_index = row * this->width; + + for (col=0; col < this->width; col += this->sampling[WIDTH_IDX]) { + unsigned int item_index = line_index + col; + + FloatIn value0 = slice1[item_index]; + + /* Test forward edges and add vertices in the current slice plane */ + if (col < (width - this->sampling[WIDTH_IDX])) { + FloatIn value = slice1[item_index + this->sampling[WIDTH_IDX]]; + + this->process_edge(value0, value, this->depth, row, col, 0, + slice0, slice1, 0); + } + + if (row < (height - this->sampling[HEIGHT_IDX])) { + /* Value from next line*/ + FloatIn value = slice1[item_index + this->width * this->sampling[HEIGHT_IDX]]; + + this->process_edge(value0, value, this->depth, row, col, 1, + slice0, slice1, 0); + } + + /* Test backward edges and add vertices in z direction */ + { + FloatIn value = slice0[item_index]; + + /* Expect forward edge, so pass: previous, current */ + this->process_edge(value, value0, + this->depth - this->sampling[DEPTH_IDX], + row, col, 2, + 0, slice0, slice1); + } + + } + } + + /* Loop over cubes to add triangle indices */ + for (row=0; row < this->height - this->sampling[HEIGHT_IDX]; row += this->sampling[HEIGHT_IDX]) { + for (col=0; col < this->width - this->sampling[WIDTH_IDX]; col += this->sampling[WIDTH_IDX]) { + unsigned char code = this->get_cell_code(slice0, slice1, + row, col); + + if (code == 0) { + continue; + } + + const int * edgeIndexPtr = &MCTriangleTable[code][0]; + for (; *edgeIndexPtr >= 0; edgeIndexPtr++) { + const unsigned int * offsets = \ + MCEdgeIndexToCoordOffsets[*edgeIndexPtr]; + + unsigned int edge_index = this->edge_index( + this->depth - this->sampling[DEPTH_IDX] + offsets[DEPTH_IDX] * this->sampling[DEPTH_IDX], + row + offsets[HEIGHT_IDX] * this->sampling[HEIGHT_IDX], + col + offsets[WIDTH_IDX] * this->sampling[WIDTH_IDX], + offsets[3]); + + /* Add vertex index to the list of indices */ + std::map<unsigned int, unsigned int>::iterator it, end; + if (offsets[DEPTH_IDX] == 0 && offsets[3] != 2) { + it = previous_edge_indices->find(edge_index); + end = previous_edge_indices->end(); + } else { + it = this->edge_indices->find(edge_index); + end = this->edge_indices->end(); + } + if (it == end) { + throw std::runtime_error( + "Internal error: cannot build triangle indices."); + } + else { + this->indices.push_back(it->second); + } + } + + } + } + + /* Clean-up previous slice cache */ + delete previous_edge_indices; + + this->depth += this->sampling[DEPTH_IDX]; +} + + +template <typename FloatIn, typename FloatOut> +void +MarchingCubes<FloatIn, FloatOut>::first_slice(const FloatIn * slice, + const FloatIn * next) +{ + assert(slice != NULL); + assert(next != NULL); + /* Init cache for this slice */ + this->edge_indices = new std::map<unsigned int, unsigned int>(); + + unsigned int row, col; + + /* Loop over slice, and add isosurface vertices in the slice plane */ + for (row=0; row < this->height; row += this->sampling[HEIGHT_IDX]) { + unsigned int line_index = row * this->width; + + for (col=0; col < this->width; col += this->sampling[WIDTH_IDX]) { + unsigned int item_index = line_index + col; + + /* For each point test forward edges */ + FloatIn value0 = slice[item_index]; + + if (col < (width - this->sampling[WIDTH_IDX])) { + FloatIn value = slice[item_index + this->sampling[WIDTH_IDX]]; + + this->process_edge(value0, value, this->depth, row, col, 0, + 0, slice, next); + } + + if (row < (height - this->sampling[HEIGHT_IDX])) { + /* Value from next line */ + FloatIn value = slice[item_index + this->width * this->sampling[HEIGHT_IDX]]; + + this->process_edge(value0, value, this->depth, row, col, 1, + 0, slice, next); + } + } + } + + this->depth += this->sampling[DEPTH_IDX]; +} + + +template <typename FloatIn, typename FloatOut> +inline unsigned int +MarchingCubes<FloatIn, FloatOut>::edge_index(const unsigned int depth, + const unsigned int row, + const unsigned int col, + const unsigned int direction) +{ + return ((depth * (this->height + 1) + row) * + (this->width + 1) + col) * 3 + direction; +} + + +template <typename FloatIn, typename FloatOut> +inline void +MarchingCubes<FloatIn, FloatOut>::process_edge(const FloatIn value0, + const FloatIn value, + const unsigned int depth, + const unsigned int row, + const unsigned int col, + const unsigned int direction, + const FloatIn * previous, + const FloatIn * current, + const FloatIn * next) +{ + assert(current != NULL); + + if ((value0 <= this->isolevel) ^ (value <= this->isolevel)) { + + /* Crossing iso-surface, store it */ + FloatIn offset = (this->isolevel - value0) / (value - value0); + + /* Store edge to vertex index correspondance */ + unsigned int edge_index = this->edge_index(depth, row, col, direction); + (*this->edge_indices)[edge_index] = this->vertices.size() / 3; + + /* Store vertex as (z, y, x) */ + if (direction == 0) { + this->vertices.push_back((FloatOut) depth); + this->vertices.push_back((FloatOut) row); + this->vertices.push_back( + (FloatOut) col + offset * this->sampling[WIDTH_IDX]); + } + else if (direction == 1) { + this->vertices.push_back((FloatOut) depth); + this->vertices.push_back( + (FloatOut) row + offset * this->sampling[HEIGHT_IDX]); + this->vertices.push_back((FloatOut) col); + } + else if (direction == 2) { + this->vertices.push_back( + (FloatOut) depth + offset * this->sampling[DEPTH_IDX]); + this->vertices.push_back((FloatOut) row); + this->vertices.push_back((FloatOut) col); + } else { + throw std::runtime_error( + "Internal error: dimension > 3, never event."); + } + + /* Store normal as (nz, ny, nx) */ + FloatOut nz, ny, nx; + const FloatIn * slice0 = (previous != 0) ? previous : current; + const FloatIn * slice1 = (previous != 0) ? current : next; + + unsigned int row_offset = this->width * this->sampling[HEIGHT_IDX]; + + if (direction == 0) { + { /* nz */ + unsigned int item, item_next_col; + + item = row * this->width + col; + if (col >= this->width - this->sampling[WIDTH_IDX]) { + /* For last column, use previous column */ + item -= this->sampling[WIDTH_IDX]; + } + item_next_col = item + this->sampling[WIDTH_IDX]; + + nz = ((1. - offset) * (slice1[item] - slice0[item]) + + offset * (slice1[item_next_col] - slice0[item_next_col])); + } + + { /* ny */ + unsigned int item, item_next_col; + + item = row * this->width + col; + if (row >= this->height - this->sampling[HEIGHT_IDX]) { + /* For last row, use previous row */ + item -= row_offset; + } + if (col >= this->width - this->sampling[WIDTH_IDX]) { + /* For last column, use previous column */ + item -= this->sampling[WIDTH_IDX]; + } + item_next_col = item + this->sampling[WIDTH_IDX]; + + ny = ((1. - offset) * (current[item + row_offset] - + current[item]) + + offset * (current[item_next_col + row_offset] - + current[item_next_col])); + } + + nx = value - value0; + + } else if (direction == 1) { + { /* nz */ + unsigned int item, item_next_row; + + item = row * this->width + col; + if (row >= this->height - this->sampling[HEIGHT_IDX]) { + /* For last row, use previous row */ + item -= row_offset; + } + item_next_row = item + row_offset; + + nz = ((1. - offset) * (slice1[item] - slice0[item]) + + offset * (slice1[item_next_row] - slice0[item_next_row])); + } + + ny = value - value0; + + { /* nx */ + unsigned int item, item_next_row; + + item = row * this->width + col; + if (row >= this->height - this->sampling[HEIGHT_IDX]) { + /* For last row, use previous row */ + item -= row_offset; + } + if (col >= this->width - this->sampling[WIDTH_IDX]) { + /* For last column, use previous column */ + item -= this->sampling[WIDTH_IDX]; + } + + item_next_row = item + row_offset; + + nx = ((1. - offset) * (current[item + this->sampling[WIDTH_IDX]] - current[item]) + + offset * (current[item_next_row + this->sampling[WIDTH_IDX]] - current[item_next_row])); + } + + } else { /* direction == 2 */ + assert(direction == 2); + /* Previous should always be 0, only here in case this changes */ + const FloatIn * other_slice = (previous != 0) ? previous : next; + + nz = value - value0; + + { /* ny */ + unsigned int item, item_next_row; + + item = row * this->width + col; + if (row >= this->height - this->sampling[HEIGHT_IDX]) { + /* For last row, use previous row */ + item -= row_offset; + } + item_next_row = item + row_offset; + + ny = ((1. - offset) * (current[item_next_row] - current[item]) + + offset * (other_slice[item_next_row] - other_slice[item])); + } + + { /* nx */ + unsigned int item; + + item = row * this->width + col; + if (col >= this->width - this->sampling[WIDTH_IDX]) { + /* For last column, use previous column */ + item -= this->sampling[WIDTH_IDX]; + } + const unsigned int item_next_col = item + this->sampling[WIDTH_IDX]; + + nx = ((1. - offset) * (current[item_next_col] - current[item]) + + offset * (other_slice[item_next_col] - other_slice[item])); + } + } + + /* apply sampling scaling */ + nz /= (FloatOut) this->sampling[0]; + ny /= (FloatOut) this->sampling[1]; + nx /= (FloatOut) this->sampling[2]; + + /* normalisation */ + FloatOut norm = sqrt(nz * nz + ny * ny + nx * nx); + if (this->invert_normals) { /* Normal inversion */ + norm *= -1.; + } + + if (norm != 0) { + nz /= norm; + ny /= norm; + nx /= norm; + } + this->normals.push_back(nz); + this->normals.push_back(ny); + this->normals.push_back(nx); + } +} + + +template <typename FloatIn, typename FloatOut> +inline unsigned char +MarchingCubes<FloatIn, FloatOut>::get_cell_code(const FloatIn * slice1, + const FloatIn * slice2, + const unsigned int row, + const unsigned int col) +{ + assert(slice1 != NULL); + assert(slice2 != NULL); + unsigned int item = row * this->width + col; + unsigned int item_next_row = item + this->width * this->sampling[HEIGHT_IDX]; + unsigned char code = 0; + + /* Cube convention for cell code: + * WARNING: This differ from layout in memory + * + * 4 +------+ 5 + * /| /| + * / | / | + * 7 +------+ 6| + * | | | | + * |0 +---|--+ 1 + * | / | / + * |/ |/ + * 3 +------+ 2 + * + */ + /* First slice */ + if (slice1[item] <= this->isolevel) { + code |= 1 << 0; + } + if (slice1[item + this->sampling[WIDTH_IDX]] <= this->isolevel) { + code |= 1 << 1; + } + if (slice1[item_next_row + this->sampling[WIDTH_IDX]] <= this->isolevel) { + code |= 1 << 2; + } + if (slice1[item_next_row] <= this->isolevel) { + code |= 1 << 3; + } + + /* Second slice */ + if (slice2[item] <= this->isolevel) { + code |= 1 << 4; + } + if (slice2[item + this->sampling[WIDTH_IDX]] <= this->isolevel) { + code |= 1 << 5; + } + if (slice2[item_next_row + this->sampling[WIDTH_IDX]] <= this->isolevel) { + code |= 1 << 6; + } + if (slice2[item_next_row] <= this->isolevel) { + code |= 1 << 7; + } + + return code; +} + +#endif /*__mc_HPP__*/ |