Filters.h

Go to the documentation of this file.

#ifndef CGAL_SURFACE_MESH_SEGMENTATION_FILTERS_H
#define CGAL_SURFACE_MESH_SEGMENTATION_FILTERS_H

#include <vector>
#include <map>
#include <algorithm>
#include <queue>

namespace CGAL {
namespace internal {

template<class Polyhedron>
class Neighbor_selector_by_edge;

template<class Polyhedron>
class Neighbor_selector_by_vertex;
  
template <class Polyhedron, class NeighborSelector = Neighbor_selector_by_edge<Polyhedron> >
class Bilateral_filtering
{
public:
    template<class ValuePropertyMap>
    void operator()(const Polyhedron& mesh,
                    int window_size,
                    ValuePropertyMap values) const
    {    
        typedef typename Polyhedron::Facet_const_handle Facet_const_handle; 
        typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator; 
        
        std::vector<double> smoothed_values; // holds smoothed values
        smoothed_values.reserve(mesh.size_of_facets());    
        
        for(Facet_const_iterator facet_it = mesh.facets_begin(); facet_it != mesh.facets_end(); ++facet_it)
        {               
            std::map<Facet_const_handle, int> neighbors;
            NeighborSelector()(facet_it, window_size, neighbors); // gather neighbors in the window
            
            // calculate deviation for range weighting.
            double current_sdf_value = values[facet_it];
            double deviation = 0.0;
            for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin(); it != neighbors.end(); ++it)
            {
                deviation += std::pow(values[it->first] - current_sdf_value, 2);
            }            
            deviation = std::sqrt(deviation / neighbors.size()); 
            if(deviation == 0.0) 
            { 
                //this might happen. In case there is no neighbors (i.e. NeighborSelector() returns just the parameter facet)
                //                 . Or all neighbor facets have same sdf value.
                smoothed_values.push_back(current_sdf_value); 
                continue;
            }    
            // smooth
            double total_sdf_value = 0.0, total_weight = 0.0;
            for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin(); it != neighbors.end(); ++it)
            {
                double spatial_weight = gaussian_function(it->second, window_size / 2.0); // window_size => 2*sigma                               
                double range_weight = gaussian_function(values[it->first] - current_sdf_value, 1.5 * deviation);
                // we can use just spatial_weight for Gauissian filtering          
                double weight = spatial_weight * range_weight;        
                     
                total_sdf_value += values[it->first] * weight;
                total_weight += weight;
            }   
            smoothed_values.push_back(total_sdf_value / total_weight);        
        }
        // put smoothed values back again to values pmap.
        std::vector<double>::iterator smoothed_value_it = smoothed_values.begin();
        for(Facet_const_iterator facet_it = mesh.facets_begin(); facet_it != mesh.facets_end();
            ++facet_it, ++smoothed_value_it)
        {
            values[facet_it] = *smoothed_value_it;
        }
    }  
private:
    double gaussian_function(double value, double deviation) const
    {
        return exp(-0.5 * (std::pow(value / deviation, 2)));
    }
};

template <class Polyhedron, class NeighborSelector = Neighbor_selector_by_vertex<Polyhedron> >
class Median_filtering
{
public:
    template<class ValuePropertyMap>
    void operator()(const Polyhedron& mesh,
                    int window_size,
                    ValuePropertyMap values) const
    {     
        typedef typename Polyhedron::Facet_const_handle Facet_const_handle; 
        typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator; 
        
        std::vector<double> smoothed_values;
        smoothed_values.reserve(mesh.size_of_facets());    
        for(Facet_const_iterator facet_it = mesh.facets_begin(); facet_it != mesh.facets_end(); ++facet_it)
        {            
            std::map<Facet_const_handle, int> neighbors;
            NeighborSelector()(facet_it, window_size, neighbors); // gather neighbors in the window
            
            std::vector<double> neighbor_values;
            neighbor_values.reserve(neighbors.size());
            for(typename std::map<Facet_const_handle, int>::iterator it = neighbors.begin(); it != neighbors.end(); ++it)
            {
                neighbor_values.push_back(values[it->first]);
            }
            // Find median.            
            int half_neighbor_count = neighbor_values.size() / 2;
            std::nth_element(neighbor_values.begin(), neighbor_values.begin() + half_neighbor_count, neighbor_values.end());
            double median_sdf = neighbor_values[half_neighbor_count];
            if(neighbor_values.size() % 2 == 0)
            {
                median_sdf += *std::max_element(neighbor_values.begin(), neighbor_values.begin() + half_neighbor_count);
                median_sdf /= 2;
            }
            smoothed_values.push_back(median_sdf);        
        }
        // put smoothed values back again to values pmap.
        std::vector<double>::iterator smoothed_value_it = smoothed_values.begin();
        for(Facet_const_iterator facet_it = mesh.facets_begin(); facet_it != mesh.facets_end(); 
            ++facet_it, ++index)
        {
            values[facet_it] = *smoothed_value_it;
        }        
    }  
};

template<class Polyhedron>
class Neighbor_selector_by_edge
{  
private:    
    typedef typename Polyhedron::Facet::Halfedge_around_facet_const_circulator Halfedge_around_facet_const_circulator;
public:
    typedef typename Polyhedron::Facet_const_handle Facet_const_handle; 
    void operator()(Facet_const_handle facet, 
                    int max_level, 
                    std::map<Facet_const_handle, int>& neighbors) const
    {
        typedef std::pair<Facet_const_handle, int> Facet_level_pair;       
        
        std::queue<Facet_level_pair> facet_queue;
        facet_queue.push(Facet_level_pair(facet, 0));
        neighbors.insert(facet_queue.front());
        
        while(!facet_queue.empty())
        {
            const Facet_level_pair& pair = facet_queue.front();  
                      
            Halfedge_around_facet_const_circulator facet_circulator = pair.first->facet_begin();
            do {
                if(!facet_circulator->opposite()->is_border())
                {
                    Facet_level_pair new_pair(facet_circulator->opposite()->facet(), pair.second + 1);
                    if(neighbors.insert(new_pair).second && new_pair.second < max_level) // first insert new_pair to map
                    {                                                                    // if insertion is OK, then check its level
                        facet_queue.push(new_pair);                                      // if its level is equal to max_level do not put it in
                    }                                                                    // queue since we do not want to traverse its neighbors
                }
            } while(++facet_circulator != pair.first->facet_begin());
            
            facet_queue.pop();
        }
    }
};

template<class Polyhedron>
class Neighbor_selector_by_vertex
{ 
private:    
    typedef typename Polyhedron::Facet::Halfedge_around_vertex_const_circulator Halfedge_around_vertex_const_circulator; 
    typedef typename Polyhedron::Halfedge_const_iterator Halfedge_const_iterator;
    typedef typename Polyhedron::Vertex_const_iterator   Vertex_const_iterator;
public:
    typedef typename Polyhedron::Facet_const_handle Facet_const_handle;
    void operator()(Facet_const_handle facet, 
                    int max_level, 
                    std::map<Facet_const_handle, int>& neighbors) const
    {
        typedef std::pair<Facet_const_handle, int> Facet_level_pair;
        
        std::queue<Facet_level_pair> facet_queue;
        facet_queue.push(Facet_level_pair(facet, 0));
        neighbors.insert(facet_queue.front());
        
        while(!facet_queue.empty())
        {
            const Facet_level_pair& pair = facet_queue.front();

            Facet_const_handle facet_front = pair.first;
            Halfedge_const_iterator edge = facet_front->halfedge();
            do { // loop on three vertices of the facet
                Vertex_const_iterator vertex = edge->vertex();
                Halfedge_around_vertex_const_circulator vertex_circulator = vertex->vertex_begin();
                do { // for each vertex loop on incoming edges (through those edges loop on neighbor facets which includes the vertex)
                    if(!vertex_circulator->is_border())
                    {
                        Facet_level_pair new_pair(vertex_circulator->opposite()->facet(), pair.second + 1);
                        if(neighbors.insert(new_pair).second && new_pair.second < max_level) // first insert new_pair to map
                        {                                                                    // if insertion is OK, then check its level
                            facet_queue.push(new_pair);                                      // if its level is equal to max_level do not put it in
                        }                                                                    // queue since we do not want to traverse its childs
                    }
                } while(++vertex_circulator != vertex->vertex_begin());
            } while((edge = edge->next()) != facet_front->halfedge());
            
            facet_queue.pop();
        }
    }
};
}//namespace internal
}//namespace CGAL
#endif //CGAL_SURFACE_MESH_SEGMENTATION_FILTERS_H