/*
 * 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;
}