diff options
Diffstat (limited to 'silx/math/marchingcubes/mc.hpp')
| -rw-r--r-- | silx/math/marchingcubes/mc.hpp | 724 |
1 files changed, 0 insertions, 724 deletions
diff --git a/silx/math/marchingcubes/mc.hpp b/silx/math/marchingcubes/mc.hpp deleted file mode 100644 index 82eced9..0000000 --- a/silx/math/marchingcubes/mc.hpp +++ /dev/null @@ -1,724 +0,0 @@ -/*########################################################################## -# -# 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__*/ |