geodesic_cpp_mod/geodesic_algorithm_graph_base.h

#ifndef GEODESIC_ALGORITHM_GRAPH_BASE_010907
#define GEODESIC_ALGORITHM_GRAPH_BASE_010907

#include "geodesic_algorithm_base.h"
#include "geodesic_mesh_elements.h"
#include <vector>
#include <set>

namespace geodesic{

template<class Node>
class GeodesicAlgorithmGraphBase: public GeodesicAlgorithmBase
{
public:
    typedef Node* node_pointer;

    GeodesicAlgorithmGraphBase(geodesic::Mesh* mesh):
        GeodesicAlgorithmBase(mesh)
    {};

    ~GeodesicAlgorithmGraphBase(){};

    void propagate(std::vector<SurfacePoint>& sources,
                   double max_propagation_distance = GEODESIC_INF,          //propagation algorithm stops after reaching the certain distance from the source
                   std::vector<SurfacePoint>* stop_points = NULL); //or after ensuring that all the stop_points are covered

    void trace_back(SurfacePoint& destination,      //trace back piecewise-linear path
        std::vector<SurfacePoint>& path);

    unsigned best_source(SurfacePoint& point,           //quickly find what source this point belongs to and what is the distance to this source
                            double& best_source_distance); 

    void print_statistics()
    {
        GeodesicAlgorithmBase::print_statistics();

        double memory = m_nodes.size()*sizeof(Node);
        std::cout << "uses about " << memory/1e6 << "Mb of memory" <<std::endl;
    }

protected:
    unsigned node_index(vertex_pointer v)       //gives index of the node that corresponds to this vertex
    {
        return v->id();
    };

    void set_sources(std::vector<SurfacePoint>& sources)
    {
        m_sources = sources;
    }

    node_pointer best_first_node(SurfacePoint& point, double& best_total_distance)
    {
        node_pointer best_node = NULL;  
        if(point.type() == VERTEX)      
        {
            vertex_pointer v = (vertex_pointer)point.base_element();
            best_node = &m_nodes[node_index(v)];
            best_total_distance = best_node->distance_from_source(); 
        }
        else
        {
            std::vector<node_pointer> possible_nodes;
            list_nodes_visible_from_source(point.base_element(), possible_nodes);

            best_total_distance = GEODESIC_INF;
            for(unsigned i=0; i<possible_nodes.size(); ++i)
            {
                node_pointer node = possible_nodes[i];

                double distance_from_dest = node->distance(&point);
                if(node->distance_from_source() + distance_from_dest < best_total_distance)
                {
                    best_total_distance = node->distance_from_source() + distance_from_dest;
                    best_node = node;
                }
            }
        }

        //FEASSERT(best_node);
        //FEASSERT(best_total_distance<GEODESIC_INF);
        if(best_total_distance > m_propagation_distance_stopped)        //result is unreliable
        {
            best_total_distance = GEODESIC_INF;
            return NULL;
        }
        else
        {
            return best_node;
        }
    };  //quickly find what node will be the next one in geodesic path

    bool check_stop_conditions(unsigned& index);        //check when propagation should stop

    virtual void list_nodes_visible_from_source(MeshElementBase* p, 
                                                std::vector<node_pointer>& storage) = 0;        //list all nodes that belong to this mesh element

    virtual void list_nodes_visible_from_node(node_pointer node,            //list all nodes that belong to this mesh element
                                              std::vector<node_pointer>& storage,
                                              std::vector<double>& distances, 
                                              double threshold_distance) = 0;   //list only the nodes whose current distance is larger than the threshold

    std::vector<Node> m_nodes;  //nodes of the graph 

    typedef std::set<node_pointer, Node> queue_type;
    queue_type m_queue;

    std::vector<SurfacePoint> m_sources;        //for simplicity, we keep sources as they are
};

template<class Node>
void GeodesicAlgorithmGraphBase<Node>::propagate(std::vector<SurfacePoint>& sources,
                                                 double max_propagation_distance,           //propagation algorithm stops after reaching the certain distance from the source
                                                 std::vector<SurfacePoint>* stop_points) //or after ensuring that all the stop_points are covered
{
    set_stop_conditions(stop_points, max_propagation_distance);
    set_sources(sources);
    
    m_queue.clear();
    m_propagation_distance_stopped = GEODESIC_INF;
    for(unsigned i=0; i<m_nodes.size(); ++i)
    {
        m_nodes[i].clear();
    }

    clock_t start = clock();

    std::vector<node_pointer> visible_nodes;        //initialize vertices directly visible from sources
    for(unsigned i=0; i<m_sources.size(); ++i)
    {
        SurfacePoint* source = &m_sources[i];
        list_nodes_visible_from_source(source->base_element(), 
                                       visible_nodes);          

        for(unsigned j=0; j<visible_nodes.size(); ++j)
        {
            node_pointer node = visible_nodes[j];
            double distance = node->distance(source);
            if(distance < node->distance_from_source())
            {
                node->distance_from_source() = distance;
                node->source_index() = i;
                node->previous() = NULL;
            }
        }
        visible_nodes.clear();
    }

    for(unsigned i=0; i<m_nodes.size(); ++i)        //initialize the queue
    {
        if(m_nodes[i].distance_from_source() < GEODESIC_INF)
        {
            m_queue.insert(&m_nodes[i]);
        }
    }

    unsigned counter = 0;
    unsigned satisfied_index = 0;
    
    std::vector<double> distances_between_nodes;
    while(!m_queue.empty())                 //main cycle
    {
        if(counter++ % 10 == 0)     //check if we covered all required vertices
        {
            if (check_stop_conditions(satisfied_index))
            {
                break;
            }
        }

        node_pointer min_node = *m_queue.begin();
        m_queue.erase(m_queue.begin());
        FEASSERT(min_node->distance_from_source() < GEODESIC_INF);

        visible_nodes.clear();
        distances_between_nodes.clear();
        list_nodes_visible_from_node(min_node, 
                                     visible_nodes, 
                                     distances_between_nodes, 
                                     min_node->distance_from_source());

        for(unsigned i=0; i<visible_nodes.size(); ++i)      //update all the adgecent vertices
        {
            node_pointer next_node = visible_nodes[i];

            if(next_node->distance_from_source() > min_node->distance_from_source() +
                                                   distances_between_nodes[i])
            {
                if(next_node->distance_from_source() < GEODESIC_INF)        //remove it from the queue
                {
                    typename queue_type::iterator iter = m_queue.find(next_node);
                    FEASSERT(iter != m_queue.end());
                    m_queue.erase(iter);
                }
                next_node->distance_from_source() = min_node->distance_from_source() +
                                                    distances_between_nodes[i];
                next_node->source_index() = min_node->source_index();
                next_node->previous() = min_node;
                m_queue.insert(next_node);
            }
        }
    } 

    m_propagation_distance_stopped = m_queue.empty() ? GEODESIC_INF : (*m_queue.begin())->distance_from_source();
    clock_t finish = clock();
    m_time_consumed = (static_cast<double>(finish)-static_cast<double>(start))/CLOCKS_PER_SEC;
    //std::cout << std::endl;
}

template<class Node>
inline bool GeodesicAlgorithmGraphBase<Node>::check_stop_conditions(unsigned& index)
{
    double queue_min_distance = (*m_queue.begin())->distance_from_source();

    if(queue_min_distance < m_max_propagation_distance)
    {
        return false;
    }

    while(index < m_stop_vertices.size())
    {
        vertex_pointer v = m_stop_vertices[index].first;
        Node& node = m_nodes[node_index(v)];
        if(queue_min_distance < node.distance_from_source() + m_stop_vertices[index].second)
        {
            return false;
        }
        ++index;
    }
    return true;
}

template<class Node>
inline void GeodesicAlgorithmGraphBase<Node>::trace_back(SurfacePoint& destination,     //trace back piecewise-linear path
                                                         std::vector<SurfacePoint>& path) 
{
    path.clear();

    double total_path_length;
    node_pointer node = best_first_node(destination, total_path_length);

    if(total_path_length>GEODESIC_INF/2.0)      //unable to find the path
    {
        return;
    }

    path.push_back(destination);

    if(node->distance(&destination) > 1e-50)
    {
        path.push_back(node->surface_point());
    }

    while(node->previous())         //follow the path
    {
        node = node->previous();
        path.push_back(node->surface_point());
    }

    SurfacePoint& source = m_sources[node->source_index()];     //add source to the path if it is not already there
    if(node->distance(&source) > 1e-50)
    {
        path.push_back(source);
    }
}


template<class Node>
inline unsigned GeodesicAlgorithmGraphBase<Node>::best_source(SurfacePoint& point,          //quickly find what source this point belongs to and what is the distance to this source
                                                                 double& best_source_distance)
{
    node_pointer node = best_first_node(point, best_source_distance);
    return node ? node->source_index() : 0;
};

}       //geodesic

#endif //GEODESIC_ALGORITHM_GRAPH_BASE_010907