D:/Programmierung/C/TRTK/include/TRTK/SurfaceExtraction3D.hpp

Go to the documentation of this file.

/*
    Segmentation of simply connected spaces in 3D data.

    Copyright (C) 2010 - 2014 Christoph Haenisch

    Chair of Medical Engineering (mediTEC)
    RWTH Aachen University
    Pauwelsstr. 20
    52074 Aachen
    Germany

    Version 0.1.0 (2011-10-14)
*/

#ifndef SURFACE_EXTRACTION_3D_7431089643
#define SURFACE_EXTRACTION_3D_7431089643

#include <cstddef>

#include "Coordinate.hpp"


namespace TRTK
{


template <class DataType>
class SurfaceExtraction3D
{
public:
    typedef DataType data_type;

    SurfaceExtraction3D();
    SurfaceExtraction3D(const DataType * data, const unsigned width, const unsigned height, const unsigned depth);
    virtual ~SurfaceExtraction3D();

    void compute();

    void setData(const DataType * data, const unsigned width, const unsigned height, const unsigned depth);
    const DataType * getSurface() const;

private:
    void initSurfaceData();
    unsigned int toIndex(const DataType x, const DataType y, const DataType z);

    unsigned height;
    unsigned width;
    unsigned depth;

    const DataType * data;
    DataType * surfaceData;
};


template <class DataType>
SurfaceExtraction3D<DataType>::SurfaceExtraction3D() :
    height(0),
    width(0),
    depth(0),
    data(NULL),
    surfaceData(NULL)
{
}


template <class DataType>
SurfaceExtraction3D<DataType>::SurfaceExtraction3D(const DataType * data, const unsigned width, const unsigned height, const unsigned depth) :
    height(height),
    width(width),
    depth(depth),
    data(data)
{
    surfaceData = new DataType[width * height * depth];
    initSurfaceData();
}


template <class DataType>
SurfaceExtraction3D<DataType>::~SurfaceExtraction3D()
{
    if (surfaceData) delete[] surfaceData;
}


template <class DataType>
void SurfaceExtraction3D<DataType>::compute()
{
    assert(width  > 0);
    assert(height > 0);
    assert(depth  > 0);

    // Do nothing in the case of a 0-dimensional data set.

    if (width * height * depth == 0) return;

    // Copy the boundary elements (sides of the data volume) since these are
    // surface elements by definition.

    // Note: toIndex(x, y, z) = z * width * height + y * width + x

    unsigned index = 0;

    // front and back

    for (unsigned x = 0; x < width; ++x)
    {
        for (unsigned y = 0; y < height; ++y)
        {
            index = toIndex(x, y, 0);
            surfaceData[index] = data[index];

            index = toIndex(x, y, depth - 1);
            surfaceData[index] = data[index];
        }
    }

    // top and bottom

    for (unsigned x = 0; x < width; ++x)
    {
        for (unsigned z = 0; z < depth; ++z)
        {
            index = toIndex(x, 0, z);
            surfaceData[index] = data[index];

            index = toIndex(x, height - 1, z);
            surfaceData[index] = data[index];
        }
    }

    // left and right side

    for (unsigned y = 0; y < height; ++y)
    {
        for (unsigned z = 0; z < depth; ++z)
        {
            index = toIndex(0, y, z);
            surfaceData[index] = data[index];

            index = toIndex(width - 1, y, z);
            surfaceData[index] = data[index];
        }
    }

    // Compute the surface elements.

    for (unsigned z = 1; z < depth - 1; ++z)
    {
        for (unsigned y = 1; y < height - 1; ++y)
        {
            for (unsigned x = 1; x < width - 1; ++x)
            {
                // If there is a voxel in the neighborhood which does not own
                // the same label (might be zero or another region label), then
                // the current voxel must be a surface voxel.

                const data_type label = data[toIndex(x, y, z)];

                if (!(data[toIndex(x - 1, y, z)] == label &&
                      data[toIndex(x + 1, y, z)] == label &&
                      data[toIndex(x, y - 1, z)] == label &&
                      data[toIndex(x, y + 1, z)] == label &&
                      data[toIndex(x, y, z - 1)] == label &&
                      data[toIndex(x, y, z + 1)] == label ))
                {
                    // surface voxel

                    surfaceData[toIndex(x, y, z)] = data[toIndex(x, y, z)];
                }
                else
                {
                    // voxel completely surrounded by other (identically labeled) voxels
                    // --> must be inner volume voxel

                    surfaceData[toIndex(x, y, z)] = 0;
                }
            }
        }
    }
}

template <class DataType>
const DataType * SurfaceExtraction3D<DataType>::getSurface() const
{
    return surfaceData;
}


template <class DataType>
void SurfaceExtraction3D<DataType>::initSurfaceData()
{
    for (unsigned i = 0; i < width * height * depth; ++i)
    {
        surfaceData[i] = DataType(0);
    }
}


template <class DataType>
void SurfaceExtraction3D<DataType>::setData(const DataType * data, const unsigned width, const unsigned height, const unsigned depth)
{
    this->data = data;
    this->width = width;
    this->height = height;
    this->depth = depth;

    if (surfaceData) delete[] surfaceData;

    surfaceData = new DataType[width * height * depth];
    initSurfaceData();
}


template <class DataType>
inline unsigned int SurfaceExtraction3D<DataType>::toIndex(const DataType x, const DataType y, const DataType z)
{
    return z * width * height + y * width + x;
}


} // namespace TRTK

#endif // SURFACE_EXTRACTION_3D_7431089643