IntersectionPoint.h

//===========================================================================
// GoTools - SINTEF Geometry Tools
//
// GoTools module: Intersections, version 1.0
//
// Copyright (C) 2000-2007 SINTEF ICT, Applied Mathematics, Norway.
//
// This program 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 version 2 of the License. 
//
// This program 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 this program; if not, write to the Free Software            
// Foundation, Inc.,                                                      
// 59 Temple Place - Suite 330,                                           
// Boston, MA  02111-1307, USA.                                           
//
// Contact information: E-mail: tor.dokken@sintef.no                      
// SINTEF ICT, Department of Applied Mathematics,                         
// P.O. Box 124 Blindern,                                                 
// 0314 Oslo, Norway.                                                     
//
// Other licenses are also available for this software, notably licenses
// for:
// - Building commercial software.                                        
// - Building software whose source code you wish to keep private.        
//===========================================================================
#ifndef _INTERSECTIONPOINT_H
#define _INTERSECTIONPOINT_H


#include "IntersectionPointUtils.h"
#include "LinkType.h"
#include "GeoTol.h"
#include "SecondOrderProperties.h"
#include <boost/shared_ptr.hpp>
#include <vector>
#include <map>


namespace Go {



class ParamObjectInt;
class ParamGeomInt;
class ParamSurface;
class IntersectionLink;



class IntersectionPoint {
public:
    IntersectionPoint(const ParamObjectInt* obj1, 
                      const ParamObjectInt* obj2,
                      const boost::shared_ptr<GeoTol> epsge,                  
                      const double* obj1_params,
                      const double* obj2_params);

    IntersectionPoint(const ParamObjectInt* obj1,
                      const ParamObjectInt* obj2,
                      const boost::shared_ptr<IntersectionPoint> ip,
                      int missing_param);

    ~IntersectionPoint();

    IntersectionPoint() {} // create an undefined point which can be
                           // assigned or read() into.
    
    void write(std::ostream& os) const;

    // @@This function is for debugging.  Do not
    // try to do anything serious with it.  It must create 'new'
    // objects for its pointers, and thus causes memory leaks!
    void read(std::istream& is);

    void writeParams(std::ostream& os) const;

    void writeInfo() const;

    void replaceParameter(double *param);
    
    Point getPoint() const;

    Point getPoint1() const
    { return point1_; }

    Point getPoint2() const
    { return point2_; }

    void projectToParamPlanes(const Point& dir, 
                              Point& par_1_dir, Point& par_2_dir) const;

    bool hasUniqueTangentDirection() const; 

    bool tangentIsOriented() const;

    double getDist() const
    { return dist_; }

    const std::vector<double>& getPar() const
    { return par_; }

    double getPar(int idx) const
    { return par_[idx]; }

    const double* getPar1() const;

    const double* getPar2() const;

    Point getPar1Point() const;
    Point getPar2Point() const;

    int numParams1() const;

    int numParams2() const;

    Point getTangent(bool second_branch = false) const;

    Point getPar1Dir(bool second_branch = false) const;

    Point getPar2Dir(bool second_branch = false) const;

    const ParamObjectInt* getObj1() const
    { return obj1_; }

    const ParamObjectInt* getObj2() const
    { return obj2_; }

    double startParam(int pardir);

    double endParam(int pardir);

    bool pointIsSingular() const
    { return getSingularityType() != ORDINARY_POINT; }

    bool isNearSingular() const;

    double parameterTolerance(int pardir); // pardir == 0 || pardir == 1

    void fixTangentOrientation(bool flip); 

    void isBoundaryPoint(bool& first, bool& second) const;

    void tangentPointingInwards(bool& first, 
                                bool& second,
                                bool& first_along_boundary,
                                bool& second_along_boundary) const;

    bool isConnectedTo(boost::shared_ptr<IntersectionPoint> p) const
    { return isConnectedTo(p.get()); }

    bool isConnectedTo(const IntersectionPoint* point) const;

    boost::shared_ptr<IntersectionLink>
    connectTo(IntersectionPoint *const point, 
              LinkType type,
              boost::shared_ptr<IntersectionLink> model_link
              = boost::shared_ptr<IntersectionLink>(),
              int added_parameter_dir = -1);

    boost::shared_ptr<IntersectionLink>
    connectTo(boost::shared_ptr<IntersectionPoint> point, 
              LinkType type,
              boost::shared_ptr<IntersectionLink> model_link
              = boost::shared_ptr<IntersectionLink>(),
              int added_parameter_dir = -1)
    { return connectTo(point.get(), type,
                       model_link, added_parameter_dir); }

    void disconnectFrom(IntersectionPoint* point);
    
    int numBranches() const;

    IntPtClassification getClassification(const ParamSurface *surf1, 
                                          const ParamSurface *surf2,
                                          int branch_num = 0) const;

    IntPtClassification getClassification(const ParamSurface *par_func,
                                          double C,
                                          int branch_num = 0) const;

    SingularityType getSingularityType() const;

    boost::shared_ptr<GeoTol> getTolerance() 
    { return epsge_; }

    boost::shared_ptr<const GeoTol> getTolerance() const
    { return epsge_; }

    boost::shared_ptr<IntersectionPoint> parentPoint()
    { return parent_point_; }

    void setParentPoint(boost::shared_ptr<IntersectionPoint> p)
    {
        //ASSERT(!parent_point_); // should only be used with orphans
        parent_point_ = p;
    }

    bool hasParentPoint()
        {
            return (parent_point_.get() != 0);
        }

    void flipObjects();

    double getInfluenceArea(int dir, bool forward, bool first_outside,
                            double aeps=0.0) const;

    int inInfluenceArea(int pardir, double par, bool first_outside) const;

    void shareInfluenceAreaWith(boost::shared_ptr<IntersectionPoint> other_pt,
                                int pardir);

    // @@@ VSK. For the time being, but needs implementing when the
    // CachedInterval::inside is OK
    bool inInfluenceAreaBracket(int pardir, double parval)
    { return false; }

    boost::shared_ptr<IntersectionLink>
    getIntersectionLink(const IntersectionPoint* point);
    
    void getNeighbours(std::vector<IntersectionPoint*>& pnts) const;

    void getNeighbourLinks(std::vector<boost::
                           shared_ptr<IntersectionLink> >& links) const;

    int numNeighbours()const { return int(intersection_links_.size()); }

    template<class bool_iterator>
    void setDifferentiateFromLeft(bool_iterator it) const  
    {
        int n = numParams1() + numParams2();
        ASSERT(n == int(diff_from_left_.size()));
        for (int i = 0; i < n; ++i, ++it) {
            diff_from_left_[i] = g2_discontinuous_params_[i] && (*it);
        }
        int key = bool2int(diff_from_left_.begin(), diff_from_left_.end());
        cur_active_2nd_order_properties_ = &sec_order_properties_[key];
    }

    bool differentiatesFromLeft(int pardir) const
    {
        ASSERT(pardir >= 0 && pardir < numParams1() + numParams2());
        return diff_from_left_[pardir];
    }

    int hasIsoLinks(int* const iso_par_dirs) const;

    int commonIsoLinks(int* const iso_par_dirs) const;

    bool containsG2Discontinuity() const 
    {
        int ki;
        for (ki=0; ki < int(g2_discontinuous_params_.size()); ki++)
            if (g2_discontinuous_params_[ki])
                return true;
        return false;
    }

    bool isDegenerate() const;

    bool isSamePoint(const IntersectionPoint* point) const;

    bool checkIntersectionPoint(IntPtInfo& int_pt_info) const;

    bool isInDomain(double* frompar, double* topar) const;

    static const double tangent_tol;

private:

    // ----- Data members -----

    // NB: the IntersectionPoint has no ownership to the objects
    // pointed to by obj1_ and obj2_
    const ParamObjectInt* obj1_; // Object number one
    const ParamObjectInt* obj2_; // Object number two

    Point point1_; // Intersection point in space in object number one
    Point point2_; // Intersection point in space in object number two
    double dist_;  // Distance between the points in the two objects

    // Parameter values of this intersection point
    std::vector<double> par_;

    // If obj1_ or obj2_ has parents, then parent_point_ is the
    // corresponding intersection point
    boost::shared_ptr<IntersectionPoint> parent_point_;

    // Contains link objects to other IntersectionPoints. The links
    // are unique.
    std::vector<boost::shared_ptr<IntersectionLink> > intersection_links_;
    
    // Tolerance-related stuff
    boost::shared_ptr<GeoTol> epsge_;

    // Caching of influence area
    mutable std::vector<CachedInterval> cached_influence_area_forwards_;
    mutable std::vector<CachedInterval> cached_influence_area_backwards_;

    // Information (that might be) influenced by the second derivative
    // of the surface.
    std::vector<bool> g2_discontinuous_params_;
    mutable std::map<int, SecondOrderProperties> sec_order_properties_;
    mutable SecondOrderProperties* cur_active_2nd_order_properties_;
    mutable std::vector<bool> diff_from_left_;

    // ------ member functions -----

    void constructor_implementation(const ParamObjectInt* obj1, 
                                    const ParamObjectInt* obj2,
                                    boost::shared_ptr<GeoTol> epsge,
                                    const double* obj_1_params,
                                    const double* obj_2_params);

    inline bool tangent_2d_is_cached() const
    { return cur_active_2nd_order_properties_->tangent_2d_is_cached; }

    inline bool singularity_info_is_cached() const
    { return cur_active_2nd_order_properties_->singularity_info_is_cached; }

    // Expecting object(s) to be 2-parametric.
    void compute_parametrical_tangents() const; 
    // The two following function are called (only) by
    // compute_parametrical_tangents())
    void compute_surface_parametrical_tangents_safe() const;
    void compute_function_parametrical_tangents() const; 

    // Calculates the tangent at a discovered singular point, and
    // return the point's SingularityType
    SingularityType calculate_tangent_at_singular_point() const;

    // Check each parameter for a second order discontinuity in the
    // IntersectionPoint
    //boost::shared_ptr<ParamObjectInt> obj2,
    static std::vector<bool> 
    detect_2nd_order_discontinuities(const ParamObjectInt* o1,
                                     const ParamObjectInt* o2,
                                     const boost::shared_ptr<GeoTol> gtol,
                                     const double* obj1_par,
                                     const double* obj2_par);

    void calculate_singularity_type() const;
    void calculate_singularity_type_function() const;
    void calculate_singularity_type_surface() const; // Assuming space is 3-d
    void calculate_singularity_type_curves() const ;
    void calculate_singularity_type_point() const;
    void calculate_singularity_type_curve_surface() const;

    static double determinant(const Point&, const Point&, const Point&);
    
    // We suppose that 'T' lies in the plane spanned by 'b1' and 'b2'.
    // We want to calculate the coefficients in order to express 'T'
    // as a linear combination of 'b1' and 'b2'.

    // Calculate the vector 'T', supposedly in the linear span of 'b1'
    // and 'b2', as a linear combination of the two.  The appropriate
    // coefficients are returned in 'b1_coef' and 'b2_coef'.  If 'T'
    // is _not_ in the plane, the projected vector is used.
    static void decompose(const Point& T, 
                          const Point& b1, 
                          const Point& b2, 
                          double& b1_coef,
                          double& b2_coef);

    // Project a vector to the tangent space of a geometric object,
    // and find the parametric direction in the object's parameter
    // space that correspond to this projection.  'dir' is the vector
    // to be projected.  The object is given by 'go', and the
    // parametrization for the point where the tangent space is
    // defined, is given by 'par'.  Whether the tangent space should
    // be computed used derivation "from left" or "from right" in each
    // parameter direction is dictated by 'from_left'.  The resulting
    // direction in parametric space is returned by 'par_dir'.
    static void proj_2_pplane(const Point& dir, 
                              const double* par,
                              std::vector<bool>::const_iterator from_left,
                              const ParamGeomInt* go, 
                              Point& par_dir);

    // Takes a range of bools and converts them to an int. The
    // resulting int has the binary representation that corresponds to
    // the range of bools.
    template<class iterator>
    static int bool2int(iterator begin_range, iterator end_range)
    {
        int result = 0;
        for (int pos = 0; begin_range != end_range; ++pos, ++begin_range) {
            // If the bool at the current position is true, then set the
            // bit at position 'pos' from the right in integer 'result'
            if (*begin_range) {
                result |= (1 << pos); 
            }
        }
        return result;
    }
    
    // Takes an int and converts it to a range of bools. The resulting
    // bools correspond to the binary representation of the int.
    template<class iterator>
    static void int2bool(int number, iterator cur_pos, iterator end_range)
    {
        for (int pos = 0; cur_pos != end_range; ++pos, ++cur_pos) {
            // Set the bool at the current position to true if 'number'
            // has the bit at position 'pos' from the right turned on
            *cur_pos = number & (1 << pos);
        }
    }

    std::map<int, SecondOrderProperties> 
    generate_2nd_order_property_map(const std::vector<bool>& disc_params);

    inline bool exactly_on_point(int pardir, double par) const
    { return (fabs(par_[pardir] - par) < epsge_->getRelParRes()); }

    // Should be called whenever the IntersectionPoint is changed
    void clearCache() const
    {
        std::map<int, SecondOrderProperties>::iterator it;
        for (it = sec_order_properties_.begin();
             it != sec_order_properties_.end(); ++it) {
            (it->second).clear();
        }
        cached_influence_area_forwards_
            = cached_influence_area_backwards_
            = std::vector<CachedInterval>(par_.size());
        diff_from_left_ = std::vector<bool>(par_.size(), false);
        cur_active_2nd_order_properties_ = &sec_order_properties_[0];
    }

};


}; // namespace Go


#endif  // _INTERSECTIONPOINT_H

---