D:/Programmierung/C/TRTK/include/TRTK/RegionGrowing3D.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.4.0 (2012-07-05)
*/

#ifndef REGION_GROWING_3D_1234689743
#define REGION_GROWING_3D_1234689743

#include <cmath>
#include <cstddef>
#include <deque>
#include <list>
#include <queue>

#include "Coordinate.hpp"


namespace TRTK
{


template <class BinaryDataType, class LabelType = unsigned short, class BinaryFieldType = const BinaryDataType *>
class RegionGrowing3D
{
public:
    typedef BinaryDataType binary_value_type;
    typedef LabelType label_type;

    typedef Coordinate<unsigned> Point;
    typedef std::deque<Point> Region;
    typedef std::deque<Region> Regions;

    RegionGrowing3D();
    RegionGrowing3D(BinaryFieldType data, const unsigned width, const unsigned height, const unsigned depth);
    virtual ~RegionGrowing3D();

    void setData(BinaryFieldType data, const unsigned width, const unsigned height, const unsigned depth);
    void setNeighborhoodSize(const unsigned size);

    template <int StrideX, int StrideY, int StrideZ> void compute();
    void compute();
    void compute(const Point & seed_point);
    void compute(const std::deque<Point> & seed_points);

    const Regions & getRegions() const;
    const LabelType * getLabelMask() const;

private:

    unsigned getIndexFromPoint(const Point & point) const;
    std::list<Point> getNeighbors(const Point & point, const unsigned size = 2) const;
    Point getPointFromIndex(const unsigned index) const;
    void initLabelMask();

    unsigned height;
    unsigned width;
    unsigned depth;
    unsigned neighborhood_size;

    bool data_available;
    BinaryFieldType data;
    LabelType * label_mask;
    Regions regions;
};


template <class BinaryDataType, class LabelType, class BinaryFieldType>
RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::RegionGrowing3D() :
    height(0),
    width(0),
    depth(0),
    neighborhood_size(2),
    data_available(false),
    label_mask(NULL)
{
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::RegionGrowing3D(BinaryFieldType data, const unsigned width, const unsigned height, const unsigned depth) :
    height(height),
    width(width),
    depth(depth),
    neighborhood_size(2),
    data_available(true),
    data(data)
{
    label_mask = new LabelType[width * height * depth];
    initLabelMask();
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::~RegionGrowing3D()
{
    if (label_mask) delete[] label_mask;
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
void RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::compute()
{
    // Note: We use Point(x, y, z) and index = (z * width * height + y * width + x)
    //       interchangeably. Conversions are done with getPointFromIndex() and
    //       getIndexFromPoint().

    if (!data_available) return;

    regions.clear();
    initLabelMask();

    LabelType unique_label = 1;

    // For each point in the binary image, do...

    for (unsigned index = 0; index < width * height * depth; ++index)
    {
        // Skip the point if it is not set, i.e., if its value is zero.

        if (!data[index])
        {
            continue;
        }
        else
        {
            // If the point has a label it was already assigned to a region and
            // thus skip it. Otherwise, this point is part of a new region. Hence,
            // get a new unique region label and add the current point as well as
            // its neighbors to this new region. To be able to process all points
            // and to add new points at the same time, a queue is used (FIFO).

            if (label_mask[index])
            {
                continue;
            }

            LabelType label = unique_label++;
            regions.push_back(Region());

            std::queue<Point> queue;
            queue.push(getPointFromIndex(index));

            while (!queue.empty())
            {
                Point point = queue.front();
                queue.pop();

                // Assign the label of the current region to the current point.

                label_mask[getIndexFromPoint(point)] = label;

                // Save the point in a global list which only contains points of a certain region.

                regions[label - 1].push_back(point);

                // Now, add those points of the neighborhood to the queue which are
                // set within the binary mask but still don't carry any label. Then
                // set the label of these points to the current region. As a result,
                // these points won't be selected a second time.

                std::list<Point> neighbors = getNeighbors(point, neighborhood_size);

                for(std::list<Point>::iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); ++neighbor)
                {
                    const unsigned index = getIndexFromPoint(*neighbor);
                    if (data[index] && !label_mask[index])
                    {
                        label_mask[index] = label;
                        queue.push(*neighbor);
                    }
                }
            }
        }
    }
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
template <int StrideX, int StrideY, int StrideZ>
void RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::compute()
{
    // Note: We use Point(x, y, z) and index = (z * width * height + y * width + x)
    //       interchangeably. Conversions are done with getPointFromIndex() and
    //       getIndexFromPoint().

    if (!data_available) return;

    regions.clear();
    initLabelMask();

    LabelType unique_label = 1;

    // For each point in the binary image, do...

    for (unsigned z = 0; z < depth; z += StrideZ)
    {
        for (unsigned y = 0; y < height; y += StrideY)
        {
            for (unsigned x = 0; x < width; x += StrideX)
            {
                unsigned index = z * width * height + y * width + x;

                // Skip the point if it is not set, i.e., if its value is zero.

                if (!data[index])
                {
                    continue;
                }
                else
                {
                    // If the point has a label it was already assigned to a region and
                    // thus skip it. Otherwise, this point is part of a new region. Hence,
                    // get a new unique region label and add the current point as well as
                    // its neighbors to this new region. To be able to process all points
                    // and to add new points at the same time, a queue is used (FIFO).

                    if (label_mask[index])
                    {
                        continue;
                    }

                    LabelType label = unique_label++;
                    regions.push_back(Region());

                    std::queue<Point> queue;
                    queue.push(getPointFromIndex(index));

                    while (!queue.empty())
                    {
                        Point point = queue.front();
                        queue.pop();

                        // Assign the label of the current region to the current point.

                        label_mask[getIndexFromPoint(point)] = label;

                        // Save the point in a global list which only contains points of a certain region.

                        regions[label - 1].push_back(point);

                        // Now, add those points of the neighborhood to the queue which are
                        // set within the binary mask but still don't carry any label. Then
                        // set the label of these points to the current region. As a result,
                        // these points won't be selected a second time.

                        std::list<Point> neighbors = getNeighbors(point, neighborhood_size);

                        for(std::list<Point>::iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); ++neighbor)
                        {
                            const unsigned index = getIndexFromPoint(*neighbor);
                            if (data[index] && !label_mask[index])
                            {
                                label_mask[index] = label;
                                queue.push(*neighbor);
                            }
                        }
                    }
                }
            }
        }
    }
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
void RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::compute(const Point & seed_point)
{
    std::deque<Point> seed_points;
    seed_points.push_back(seed_point);
    compute(seed_points);
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
void RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::compute(const std::deque<Point> & seed_points)
{
    // Note: We use Point(x, y, z) and index = (z * width * height + y * width + x)
    //       interchangeably. Conversions are done with getPointFromIndex() and
    //       getIndexFromPoint().

    if (!data_available) return;

    regions.clear();
    initLabelMask();

    LabelType unique_label = 1;

    // For each point in the binary image, do...

    for (unsigned i = 0; i < seed_points.size(); ++i)
    {
        unsigned index = getIndexFromPoint(seed_points[i]);

        // Skip the point if it is not set, i.e., if its value is zero.

        if (!data[index])
        {
            continue;
        }
        else
        {
            // If the point has a label it was already assigned to a region and
            // thus skip it. Otherwise, this point is part of a new region. Hence,
            // get a new unique region label and add the current point as well as
            // its neighbors to this new region. To be able to process all points
            // and to add new points at the same time, a queue is used (FIFO).

            if (label_mask[index])
            {
                continue;
            }

            LabelType label = unique_label++;
            regions.push_back(Region());

            std::queue<Point> queue;
            queue.push(getPointFromIndex(index));


            while (!queue.empty())
            {
                Point point = queue.front();
                queue.pop();

                // Assign the label of the current region to the current point.

                label_mask[getIndexFromPoint(point)] = label;

                // Save the point in a global list which only contains points of a certain region.

                regions[label - 1].push_back(point);

                // Now, add those points of the neighborhood to the queue which are
                // set within the binary mask but still don't carry any label. Then
                // set the label of these points to the current region. As a result,
                // these points won't be selected a second time.

                std::list<Point> neighbors = getNeighbors(point, neighborhood_size);

                for(std::list<Point>::iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); ++neighbor)
                {
                    const unsigned index = getIndexFromPoint(*neighbor);
                    if (data[index] && !label_mask[index])
                    {
                        label_mask[index] = label;
                        queue.push(*neighbor);
                    }
                }
            }
        }
    }
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
inline unsigned RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::getIndexFromPoint(const Point & point) const
{
    return point.z() * width * height + point.y() * width + point.x();
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
const LabelType * RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::getLabelMask() const
{
    return label_mask;
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
std::list<typename RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::Point> RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::getNeighbors(const Point & point, const unsigned size) const
{
    // typical values for size are 0, 1, 2, 4, 8

    // Return a list with all adjacent neighbors. Also check,
    // whether a neighbor is still within the data.

    std::list<Point> neighbors;

    const int delta = size == 2 ? 1 : int(std::ceil(std::sqrt(float(size))));

    for (int delta_x = -delta; delta_x <= delta; ++delta_x)
    {
        for (int delta_y = -delta; delta_y <= delta; ++delta_y)
        {
            for (int delta_z = -delta; delta_z <= delta; ++delta_z)
            {
                if (unsigned(delta_x * delta_x + delta_y * delta_y + delta_z * delta_z) <= size)
                {
                    int x = point.x() + delta_x;
                    int y = point.y() + delta_y;
                    int z = point.z() + delta_z;

                    if (0 <= x && x < int(width) && 0 <= y && y < int(height) && 0 <= z && z < int(depth))
                    {
                        neighbors.push_back(Point(x, y, z));
                    }
                }
            }
        }
    }

    return neighbors;
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
inline typename RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::Point RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::getPointFromIndex(const unsigned index) const
{
    const unsigned remainder_xy = index % (width * height);
    const Point::value_type x = remainder_xy % width;
    const Point::value_type y = remainder_xy / width;
    const Point::value_type z = index / (width * height);

    return Point(x, y, z);
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
inline const typename RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::Regions & RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::getRegions() const
{
    return regions;
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
void RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::initLabelMask()
{
    for (unsigned i = 0; i < width * height * depth; ++i)
    {
        label_mask[i] = LabelType(0);
    }
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
void RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::setData(BinaryFieldType data, const unsigned width, const unsigned height, const unsigned depth)
{
    this->data = data;
    this->data_available = true;

    this->width = width;
    this->height = height;
    this->depth = depth;

    if (label_mask) delete[] label_mask;

    label_mask = new LabelType[width * height * depth];
    initLabelMask();
}


template <class BinaryDataType, class LabelType, class BinaryFieldType>
void RegionGrowing3D<BinaryDataType, LabelType, BinaryFieldType>::setNeighborhoodSize(const unsigned size)
{
    neighborhood_size = size;
}


} // namespace TRTK

#endif // REGION_GROWING_3D_1234689743