PluginSpheroidManager.hpp

Go to the documentation of this file.

#ifndef PLUGINSPHEROIDMANAGER_HPP
#define PLUGINSPHEROIDMANAGER_HPP
 
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/algorithm.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/intersections.h>
#include <vector>
#include <mecacell/mecacell.h>
#include <math.h>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Polyhedron_3<K> Polyhedron_3;
// define point creator
typedef K::Point_3 Point_3;
typedef K::Line_3 Line_3;
 
struct PlaneFromFacet
{
    Polyhedron_3::Plane_3 operator()(Polyhedron_3::Facet &f)
    {
        Polyhedron_3::Halfedge_handle h = f.halfedge();
        return Polyhedron_3::Plane_3(h->vertex()->point(),
                                     h->next()->vertex()->point(),
                                     h->opposite()->vertex()->point());
    }
};
 
namespace SpheroidManager {
 
    template<typename cell_t>
    class PluginSpheroidManager {
        typedef std::vector<double> Data;
 
    private:
        MecaCell::Vec centroid; 
        double optimRad = 0.0; 
        double meanRadius; 
        inline static double Kernel(const double x) { return exp(- pow(x, 2) / 2.) / sqrt(2. * M_PI); }
 
        void  updateDistanceFromCentroid(std::vector<cell_t*> cells) {
            double r = 0;
            Data data;
            std::vector<Point_3> points;
            double max = - INFINITY;
            for (auto &c : cells) {
                r = (c->getBody().getPosition() - centroid).length();
                c->getBody().setDistanceFromCentroid(r);
                data.push_back(r);
                points.push_back(Point_3(c->getBody().getPosition().x(), c->getBody().getPosition().y(), c->getBody().getPosition().z()));
                if (r > max) {
                    max = r;
                }
            }
            optimRad = getOptimalRadius(data, max);
 
            // define polyhedron to hold convex hull
            Polyhedron_3 poly;
 
            // compute convex hull of non-collinear points
            CGAL::convex_hull_3(points.begin(), points.end(), poly);
 
            if (poly.size_of_vertices() > 5) {
                std::transform(poly.facets_begin(), poly.facets_end(), poly.planes_begin(), PlaneFromFacet());
                Point_3 pCentroid = Point_3(centroid.x(), centroid.y(), centroid.z());
                for (auto &c : cells) {
                    Point_3 pCell = Point_3(c->getBody().getPosition().x(), c->getBody().getPosition().y(),
                                            c->getBody().getPosition().z());
                    Line_3 lCellCentroid = Line_3(pCell, pCentroid);
                    Polyhedron_3::Plane_iterator pit;
                    double min = INFINITY;
                    for (pit = poly.planes_begin(); pit != poly.planes_end(); ++pit) {
                        auto result = CGAL::intersection((*pit), lCellCentroid);
                        if (result) {
                            const Point_3 *p = boost::get<Point_3>(&*result);
                            double dist = CGAL::squared_distance(pCell,(*p));
                            if (dist < min) {
                                min = dist;
                            }
                        }
                    }
                    c->getBody().setDistanceFromCentroid( c->getBody().getDistanceFromCentroid() * optimRad /
                                                        (c->getBody().getDistanceFromCentroid() + sqrt(min) ));
                }
            }
        }
 
        double getOptimalRadius(const Data &data, double max) {
            double maxDensity = -1.;
            double xMaxDensity = -1.;
            double val;
            double lambda = 10.;
            double sum;
            for (double i = 0.; i <= max; i += 0.2) {
                sum = 0;
                for (size_t j = 0; j < data.size(); ++j) {
                    sum += Kernel((i - data[j]) / lambda);
                }
                val = sum / lambda;
                if (val > maxDensity) {
                    maxDensity = val;
                    xMaxDensity = i;
                }
            }
            return (max - (max - xMaxDensity) / 2 + meanRadius);
        }
 
        MecaCell::Vec computeCentroid(std::vector<cell_t*> cells) {
            centroid = MecaCell::Vec(0, 0, 0);
            for (auto &c : cells) {
                centroid += c->getBody().getPosition();
            }
            return centroid /= cells.size();
        }
 
    public:
        inline PluginSpheroidManager() = default;
 
        inline PluginSpheroidManager(double r) : meanRadius(r){}
 
        inline MecaCell::Vec getCentroid() const{ return centroid; }
 
        inline double getSpheroidRadius() const{ return (optimRad); }
 
        inline void setMeanRadius(double r) { meanRadius = r; }
 
 
        void updateSpheroid(std::vector<cell_t*> cells) {
            if (cells.size() > 5) {
                computeCentroid(cells);
                updateDistanceFromCentroid(cells);
            }
        }
 
        template<typename world_t>
        void initSpheroid(world_t *w) {
            updateSpheroid(w->newCells);
        }
 
        template<typename world_t>
        void preBehaviorUpdate(world_t *w) {
            updateSpheroid(w->cells);
        }
    };
}
#endif