source: src/Filling/Generators/NodeGenerator.cpp

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * NodeGenerator.cpp
 *
 *
 *  Created on: Jan 16, 2012
 *      Author: heber
 */


// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

//#include "CodePatterns/MemDebug.hpp"

#include "NodeGenerator.hpp"

#include <cmath>
#include <vector>

#include "Shapes/Shape.hpp"
#include "Shapes/ShapeOps.hpp"

/** Constructor for NodeGenerator.
 *
 * @param _shape Shape which is to be filled with nodes
 */
NodeGenerator::NodeGenerator(const Shape &_shape) :
        shape(_shape)
{}

/** Destructor for NodeGenerator.
 *
 */
NodeGenerator::~NodeGenerator()
{}

/** Returns a set of points contained in the \a NodeGenerator::_shape
 * with a homogeneous density.
 *
 * The central idea of the algorithm is to make use of the present function
 * for obtaining homogeneously distributed points on all surfaces of
 * presently implemented Shape's.
 *
 * The given \a _shape is shrinked such that eventually the volume is filled
 * with points. Points filled in that actually reside outside the given
 * Shape are eventually extracted as all are checked for inclusion.
 *
 * \note Although a certain density of points given as it has to be converted
 * to an discrete regime this density can in general be matched only
 * approximately.
 *
 * @param density desired density of points, the unit is the inverse of the
 * required volume per point
 * @return set of nodes
 */
NodeSet  NodeGenerator::operator()(const double density)
{
        // calculate volume and surface of the given Shape
        const double volume = shape.getVolume();
        const double surface = shape.getSurfaceArea();

        // calculate the number of shrinking operations
        const double radius = shape.getRadius();
        const double fraction = surface/volume;
        const int factor = floor(radius * fraction);

        // calculate correction for surface density due to discrete number of layers
        const double surfaceDensity = volume/(double)factor;

        // fill the shrinking vector
        std::vector<double> shrinking_factors;
        for(int f=0; f < factor; ++f)
                shrinking_factors.push_back(radius*((double)f/(double)factor));

        // go through the shrinking operations
        NodeSet nodes;
        for (std::vector<double>::const_iterator iter = shrinking_factors.begin();
                        iter != shrinking_factors.end(); ++iter) {
                const Shape currentShape = resize(shape, *iter);
                std::vector<Vector> pointsOnSurface =
                                currentShape.getHomogeneousPointsOnSurface(surfaceDensity);
                nodes.insert(nodes.end(), pointsOnSurface.begin(), pointsOnSurface.end());
        }

        // check each point whether its inside the surface
        return filterOutsidePoints(nodes);
}

/** Filters out all points that are not contained inside NodeGenerator::shape.
 *
 * @param nodes nodes to check
 * @return subset of nodes that fulfill Shape::isInisde().
 */
NodeSet NodeGenerator::filterOutsidePoints(const NodeSet &nodes) const
{
        NodeSet return_nodes;
        return_nodes.reserve(nodes.size());
        for (NodeSet::const_iterator iter = nodes.begin(); iter != nodes.end(); ++iter) {
                if (shape.isInside(*iter))
                        return_nodes.push_back(*iter);
        }

        return return_nodes;
}