path: root/silx/math/marchingcubes/mc.hpp

/*##########################################################################
#
# 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__*/