Commit 5200f8e5 authored by Jeannie Backer's avatar Jeannie Backer
Browse files

Added cnetcombinept program to the ISIS package. Fixes #3870

git-svn-id: http://subversion.wr.usgs.gov/repos/prog/isis3/trunk@6631 41f8697f-d340-4b68-9986-7bafba869bb8
parent aac3d376
Loading
Loading
Loading
Loading
+310 −0
Original line number Diff line number Diff line
#ifndef ControlPointCloudPt_h
#define ControlPointCloudPt_h
/**
 * @file
 * $Revision: 1.0 $ 
 * $Date: 2014/02/27 18:49:25 $ 
 *
 *   Unless noted otherwise, the portions of Isis written by the USGS are public
 *   domain. See individual third-party library and package descriptions for
 *   intellectual property information,user agreements, and related information.
 *
 *   Although Isis has been used by the USGS, no warranty, expressed or implied,
 *   is made by the USGS as to the accuracy and functioning of such software
 *   and related material nor shall the fact of distribution constitute any such
 *   warranty, and no responsibility is assumed by the USGS in connection
 *   therewith.
 *
 *   For additional information, launch
 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html in a browser or see
 *   the Privacy & Disclaimers page on the Isis website,
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
 *   http://www.usgs.gov/privacy.html.
 */

#include <cmath>

#include <QExplicitlySharedDataPointer>
#include <QSharedData>
#include <QVector>
#include <QtGlobal>

#include "ControlPoint.h"
#include "ControlMeasure.h"


namespace Isis {

/**
 * @brief 3-D ControlPoint class for use in PointCloud datasets
 *  
 *  The ControlPoint container is required to not change its content for the
 *  duration of use of the nanoflann kd-tree built from the points.
 *  
 * @author  2015-10-11 Kris Becker
 *  
 * @internal 
 *   @history 2015-10-11 Kris Becker - Original Version 
 */

class ControlPointCloudPt {
  public:
    enum CoordinateType { Image, Ground };
    enum Ownership { Shared,      // Indicates creator of ControlPoint retains
                                  // ownership. take() provides a unique clone
                                  // of ControlPoint.
                     Exclusive    // Creator gives ownership to ControlPointCloudPt
                                  // take() will give up pointer interally.
                   };

    ControlPointCloudPt() : m_xyz(), m_type(Image), m_serialno(),
                            m_data( new ControlPointData() ),
                            m_merged() {
      m_data->getImageCoordinates(m_xyz); 
    }

    ControlPointCloudPt(ControlPoint *point, const CoordinateType &ptype,
                        const Ownership owner,
                        const QString &serialno = "", 
                        const double &weight = 1.0) : 
                        m_xyz( ), m_type(ptype), m_serialno(serialno),
                        m_data(new ControlPointData(point, owner, weight)),
                        m_merged() {
      if ( Image == ptype ) {
        m_data->setReference(serialno);
        if ( !selectImageCoordinates() ) { 
          m_data->disable(); 
        }
      }
      else {  //  ( Ground == ptype )
        if ( !selectGroundCoordinates() ) { 
          m_data->disable(); 
        }
      }
    }

    virtual ~ControlPointCloudPt() { }

    inline bool isValid() const {
      return ( !m_data->isDisabled() );
    }

    inline void disable() {
      return ( m_data->disable() );
    }

    inline int size() const {
      return ( m_data->size() );
    }

    Ownership owner() const {
      return ( m_data->m_owner );
    }

    ControlPoint *take() const {
       return ( m_data->take() );
    }

    inline bool selectGroundCoordinates() {
      m_type = Ground;
      return ( m_data->getGroundCoordinates(m_xyz) );
    }

    inline bool selectImageCoordinates() {
      m_type = Image;
      return ( m_data->getImageCoordinates(m_xyz) );
    }

    inline CoordinateType getCoordinateType() const {
      return ( m_type );
    }

    inline QString id() const {
      return ( m_data->m_point->GetId() );
    }

    inline const ControlPoint *getPoint() const {
      return ( m_data->m_point );
    }

    inline ControlPoint &getPoint() {
      return ( *m_data->m_point );
    }

    inline QString getSerialNumber() const {
      return ( m_serialno );
    }

    inline ControlMeasure *getMeasure(const QString &serialno) const {
      if ( !isValid() ) { return (0); }
      if (m_data->m_point->HasSerialNumber(serialno) ) {
        return ( m_data->m_point->GetMeasure(serialno) );
      }

      return (0);
    }

    // Convenient retrieval coordinates
    inline double x() const { return ( m_xyz[0] ); }
    inline double y() const { return ( m_xyz[1] ); }
    inline double z() const { return ( m_xyz[2] ); }
    inline double w() const { return ( m_xyz[3] ); }
    inline const double *array() const { return ( m_xyz ); }

    inline bool operator==(const ControlPointCloudPt &other) const {
      return ( m_data->m_point == other.m_data->m_point );
    }

    inline bool operator!=(const ControlPointCloudPt &other) const {
      return ( m_data->m_point != other.m_data->m_point );
    }

  private:
    /** 
     *  Shared ControlPoint data pointer
     * @author  2015-10-11 Kris Becker
     *  
     * @internal 
     *   @history 2015-10-11 Kris Becker - Original Version 
     */
    class ControlPointData : public QSharedData {
      typedef ControlPointCloudPt::Ownership Ownership;
      public:
        ControlPointData() : QSharedData(), m_point(0), m_reference(0), 
                             m_weight(1.0), m_initial(0),
                             m_owner(Exclusive){ }
        ControlPointData(ControlPoint *point, const Ownership owner, 
                         const double weight = 1.0) :
                         QSharedData(), m_point(point), 
                         m_reference(point->GetRefMeasure()), 
                         m_weight(weight),
                         m_initial(point->GetNumValidMeasures()), 
                         m_owner(owner) { }
        ControlPointData(const ControlPointData &other) : QSharedData(other),
                         m_point(other.m_point),
                         m_reference(other.m_reference),
                         m_weight(other.m_weight),
                         m_initial(other.m_initial),
                         m_owner(Shared) { }
        ~ControlPointData() { 
          if ( Exclusive == m_owner ) {
             delete (m_point);
             m_point = 0;
          }
        }

        inline int size() const {
          if ( 0 == m_point ) { return (0); }
          return ( m_point->GetNumValidMeasures() );
        }

        inline bool setReference(const QString &serialno) {
          if ( 0 == m_point ) { return ( false); }
          m_reference = 0;

          if ( serialno.isEmpty() ) {
            m_reference = m_point->GetRefMeasure();
          }
          else {
            if ( m_point->HasSerialNumber(serialno) ) {
              m_reference = m_point->GetMeasure(serialno);
            }
          }

          return ( 0 != m_reference );
        }

        /** This approach allows use of 2 and 3 dimensional Euclidean
         *  distances */
        inline bool getImageCoordinates(double coords[4]) const {
          (void) getNoPointData(coords);
          if ( 0 == m_reference ) { return (false); }
          if ( 0 == m_point ) { return ( false ); }

          // ControlMeasure *refm = m_point->GetRefMeasure();
          coords[0] = m_reference->GetSample();
          coords[1] = m_reference->GetLine();
          coords[2] = 0.0;
          coords[3] = m_weight;
          return ( true );
        }


        /** This approach assumes 3 dimensional Euclidean distances */
        inline bool getGroundCoordinates(double coords[4]) const {
          (void) getNoPointData(coords);
          if ( 0 == m_point ) { return ( false ); }

          // Always get the best surface point (2015-10-27)
          SurfacePoint surfpt = m_point->GetBestSurfacePoint();
          if ( !surfpt.Valid() ) { return (false); }

          // Get the location and convert to meters!
          surfpt.ToNaifArray(&coords[0]);
          coords[0] *= 1000.0;
          coords[1] *= 1000.0;
          coords[2] *= 1000.0;
          return (true);
        }

        inline bool isDisabled() const {
          if ( 0 == m_point ) return ( true );
          return (
                   m_point->IsInvalid()  || 
                   m_point->IsIgnored() || 
                   m_point->IsRejected() ||
                   m_point->IsEditLocked()  ||
                   ( 0 == m_reference )
                   ); 
        }

        inline void disable() {
          if ( !isDisabled() ) {
            m_point->SetIgnored( true );
          }
          return;
        }

        ControlPoint *take() {
          ControlPoint *p(m_point); 
          if ( 0 == m_point ) {  return (p); }

          // If someone else owns the point, clone it
          if ( Shared == m_owner ) {
            p = new ControlPoint(*m_point);
          }
          else {
            // Exclusive == Ownership
            // Relinquish ownership
            m_point = 0;
          }

          return ( p );
        }

        // Data....
        ControlPoint          *m_point; 
        ControlMeasure        *m_reference;
        double                m_weight;
        int                   m_initial;
        Ownership             m_owner;


      private:
        inline bool getNoPointData(double coords[4]) const {
          coords[0] = coords[1] = coords[2] = 0.0; 
          coords[3] = m_weight;
          return ( false );
        }
    };

    // Variables...
    double                                  m_xyz[4];
    CoordinateType                          m_type;
    QString                                 m_serialno;
    QExplicitlySharedDataPointer<ControlPointData> m_data;
    QList<ControlPointCloudPt>              m_merged;
};

};  // namespace Isis
#endif
+182 −0
Original line number Diff line number Diff line
#ifndef ControlPointMerger_h
#define ControlPointMerger_h
/**
 * @file
 * $Revision: 1.0 $ 
 * $Date: 2014/02/27 18:49:25 $ 
 *
 *   Unless noted otherwise, the portions of Isis written by the USGS are public
 *   domain. See individual third-party library and package descriptions for
 *   intellectual property information,user agreements, and related information.
 *
 *   Although Isis has been used by the USGS, no warranty, expressed or implied,
 *   is made by the USGS as to the accuracy and functioning of such software
 *   and related material nor shall the fact of distribution constitute any such
 *   warranty, and no responsibility is assumed by the USGS in connection
 *   therewith.
 *
 *   For additional information, launch
 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html in a browser or see
 *   the Privacy &amp; Disclaimers page on the Isis website,
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
 *   http://www.usgs.gov/privacy.html.
 */

#include <QList>
#include <QString>
#include <QSharedPointer>

#include <opencv2/opencv.hpp>

#include "ControlPointCloudPt.h"
#include "Statistics.h"

namespace Isis {

/**
 * @brief Combine control points based upon distance criteria 
 *  
 * This class will collect and compute ControlPoint candidates that are within 
 * a pixel tolerance for merging into a single control point. 
 *  
 * The criteria applied computes statistics on all common measures within the 
 * control point for image coordinate searches.
 *  
 * @author  2015-10-11 Kris Becker
 *  
 * @internal 
 *   @history 2015-10-11 Kris Becker - Original Version 
 */

class ControlPointMerger {
  public:
    ControlPointMerger() : m_image_tolerance(DBL_MAX),
                           m_ground_tolerance(DBL_MAX),
                           m_ground_distance(DBL_MAX),
                           m_source(), m_candidates() { }
    ControlPointMerger(const double image_tolerance,
                       const double ground_tolerance = -1.0) :
                       m_image_tolerance(image_tolerance),
                       m_ground_tolerance(ground_tolerance),
                       m_ground_distance(DBL_MAX), 
                       m_source(), m_candidates() { }

    virtual ~ControlPointMerger() { }

    int size() const {
      return ( m_candidates.size() );
    }

    Statistics getImageStatistics(int index) const {
      Q_ASSERT( index >= 0 );
      Q_ASSERT ( index < size() );
      return ( m_candidates[index].second );
    }

    double getGroundDistance() const {
      return ( m_ground_distance );
    }

    void clear() {
      m_candidates.clear();
      return;
    }

    void apply(ControlPointCloudPt &source, ControlPointCloudPt &candidate,
               const double distance)  {

      if ( !source.isValid() ) { return; }
      if ( !candidate.isValid() ) { return; }
      m_source = source;

      ControlPoint &point = source.getPoint();
      m_ground_distance = ground_distance( point, candidate.getPoint() );

      QList<ControlMeasure *> measures = point.getMeasures( true );
      Statistics stats;
      BOOST_FOREACH ( ControlMeasure *m, measures ) {
        if ( isValid(*m) ) {
          ControlMeasure *c = candidate.getMeasure(m->GetCubeSerialNumber()); 
          if ( (0 != c) ) {  
            if ( isValid(*c) ) { stats.AddData(image_distance(*m, *c)); }
          }
        }
      }

      // Test for conditions of a merger. If there are common image measures,
      // use the statistics. If no common measures, use the ground distance.
      if ( stats.ValidPixels() > 0 ) {
        if (stats.Average() <= m_image_tolerance) {
          m_candidates.append( qMakePair<ControlPointCloudPt, Statistics> (candidate, stats)); 
        }
      }
      else if ( m_ground_distance <= m_ground_tolerance ) {
         m_candidates.append( qMakePair<ControlPointCloudPt, Statistics> (candidate, stats)); 
      }
      return;
    }

    int merge() {
      int nMerged(0);
      ControlPoint &source = m_source.getPoint();
      for ( int i = 0 ; i < m_candidates.size() ; i++) {
        ControlPointCloudPt &cpt = m_candidates[i].first;
        ControlPoint &candidate = cpt.getPoint();
        QList<ControlMeasure *> ms = candidate.getMeasures(true);
        for ( int m = 0 ; m < ms.size() ; m++) {
          if ( !source.HasSerialNumber(ms[m]->GetCubeSerialNumber()) ) {
            source.Add( new ControlMeasure(*ms[m]) );
            nMerged++;
          }
        }
        // Essentially disables this point
        cpt.disable();
      }

      return ( nMerged );
    }

  private:
    
    double                                         m_image_tolerance;
    double                                         m_ground_tolerance;
    double                                         m_ground_distance;
    ControlPointCloudPt                            m_source;
    QList<QPair<ControlPointCloudPt, Statistics> > m_candidates;

    inline bool isValid(const ControlMeasure &m) const {
       return ( !( m.IsIgnored() || m.IsRejected() ) );
    }

    inline double image_distance(const ControlMeasure &source, 
                                 const ControlMeasure &candidate) const {
      double dx = source.GetSample() - candidate.GetSample();
      double dy = source.GetLine() - candidate.GetLine();
      return ( std::sqrt( dx*dx + dy*dy ) );
    }

    inline double ground_distance(const ControlPoint &source, 
                                  const ControlPoint &candidate) const {
      double spts[3], cpts[3];
      getGroundVector(source, spts);
      getGroundVector(candidate, cpts);

      double dx = spts[0] - cpts[0];
      double dy = spts[1] - cpts[1];
      double dz = spts[2] - cpts[2];
      return ( std::sqrt( dx*dx + dy*dy + dz*dz ) );
    }

    /** Compute ground point in meters */
    inline void getGroundVector(const ControlPoint &point, double v[3]) const {
      // Always use the best surface point available
       point.GetBestSurfacePoint().ToNaifArray(v);
       v[0] *= 1000.0;
       v[1] *= 1000.0;
       v[2] *= 1000.0;
       return;
    }
};

}  // namespace Isis
#endif
+3 −0
Original line number Diff line number Diff line
include $(ISISROOT)/make/isismake.apps

ISISCPPFLAGS += -frounding-math
+271 −0
Original line number Diff line number Diff line
#ifndef PointCloud_h
#define PointCloud_h
/**
 * @file
 * $Revision: 1.0 $ 
 * $Date: 2014/02/27 18:49:25 $ 
 *
 *   Unless noted otherwise, the portions of Isis written by the USGS are public
 *   domain. See individual third-party library and package descriptions for
 *   intellectual property information,user agreements, and related information.
 *
 *   Although Isis has been used by the USGS, no warranty, expressed or implied,
 *   is made by the USGS as to the accuracy and functioning of such software
 *   and related material nor shall the fact of distribution constitute any such
 *   warranty, and no responsibility is assumed by the USGS in connection
 *   therewith.
 *
 *   For additional information, launch
 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html in a browser or see
 *   the Privacy &amp; Disclaimers page on the Isis website,
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
 *   http://www.usgs.gov/privacy.html.
 */

#include <QtGlobal>
#include <QVector>

namespace Isis {

/** 
 *  Functor to compute 3-d Euclidean distances
 *  
 * @author 2014-02-17 Kris Becker
 *  
 * @internal 
 *   @history 2014-02-17 Kris Becker - Original Version 
 *  
 */
template <class T> class Dist3d {
public:
  enum { Dimension = 3 };
  Dist3d() { }
  ~Dist3d() { }

  inline int dimension() const {
    return ( Dimension );
  }

  inline double operator()(const T &datum1, const T &datum2) const {
    double dx = datum1.x() - datum2.x();
    double dy = datum1.y() - datum2.y();
    double dz = datum1.z() - datum2.z();
    return ( dx*dx + dy*dy + dz*dz );
  }

  inline double operator()(const double *datum1, const T &datum2) const {
    double dx = datum1[0] - datum2.x();
    double dy = datum1[1] - datum2.y();
    double dz = datum1[2] - datum2.z();
    return ( dx*dx + dy*dy + dz*dz );
  }
};

/** 
 *  Functor to compute 2-d Euclidean distances
 *  
 * @author 2014-02-17 Kris Becker
 *  
 * @internal 
 *   @history 2014-02-17 Kris Becker - Original Version 
 */
template <class T> class Dist2d {
public:
  enum { Dimension = 2 };
  Dist2d() { }
  ~Dist2d() { }

  inline int dimension() const {
    return ( Dimension );
  }

  inline double operator()(const T &datum1, const T &datum2) const {
    double dx = datum1.x() - datum2.x();
    double dy = datum1.y() - datum2.y();
    return ( dx*dx + dy*dy );
  }

  inline double operator()(const double *datum1, const T &datum2) const {
    double dx = datum1[0] - datum2.x();
    double dy = datum1[1] - datum2.y();
    return ( dx*dx + dy*dy );
  }
};

/** 
 *  Functor to compute 1-d Manhattan distances
 *  
 * @author 2014-02-17 Kris Becker
 *  
 * @internal 
 *   @history 2014-02-17 Kris Becker - Original Version 
 */
template <class T> class Dist1d {
public:
  enum { Dimension = 1 };
  Dist1d() { }
  ~Dist1d() { }

  inline int dimension() const {
    return ( Dimension );
  }


  inline double operator()(const T &datum1, const T &datum2) const {
    double dx = datum1.x() - datum2.x();
    return ( dx*dx );
  }

  inline double operator()(const double *datum1, const T &datum2) const {
    double dx = datum1[0] - datum2.x();
    return ( dx*dx );
  }
};


/**
 * @brief Point cloud adapter class for nanoflann kd-tree interface 
 *  
 * This class provides the point cloud class for Point3d clouds that interface 
 * with the nanoflann (http://code.google.com/p/nanoflann/) kd-tree fast search 
 * query library. The reference for this approach can be found in an example 
 * found at 
 * http://nanoflann.googlecode.com/svn/trunk/examples/pointcloud_kdd_radius.cpp.
 *  
 * This class is designed to accept a reference to a container of 3-d points. 
 * In addition, this class accesses individual points from PointCloud using a 
 * vector component operator so as to standardize and complete this interface: 
 *  
 *   @code
 *   double x() const;  // X component of point
 *   double y() const;  // Y component of point
 *   double z() const;  // Z component of point
 *   double w() const;  // Optional weight of point (default should be 1.0)
 *  @endcode
 *  
 *  The point container is required to not change its content for the duration
 *  of use of the nanoflann kd-tree built from the points. Because of this,
 *  there is no clear() method to discard existing points.
 *  
 *  This class supports 2D and 3D Euclidean distant calculations. This option is
 *  specified in the PointCloudTree class constructor when building the kd-tree
 *  index.
 *  
 *  The routines kdtree_get_point_count(), kdtree_distance(), kdtree_get_pt()
 *  and kdtree_get_bbox() (default implementation) satisfy the needs of the
 *  Nanoflann kd-tree template library.
 *  
 *  This class is not a template, but all its methods are inlined for efficiency
 *  reasons. This implemetation approach allows us to take best advantage of the
 *  optimization that the Nanoflann library offers.
 *  
 *  The Point3d class is also optimized in the same fashion and offers a
 *  flexible point base class to complete the implementation.
 *  
 *  This point cloud class is designed with the body-fixed coordinate system in
 *  mind. Therfore, the units of the point vectors is assumed to be kilometers
 *  but this is not required - as long as Euclidean distances apply to the point
 *  dataset, any 3D vector representation could utilize this class.
 *  
 * @author 2014-02-17 Kris Becker
 *  
 * @internal 
 *   @history 2014-02-17 Kris Becker - Original Version 
 */
template <class T, class Distance = Dist3d<T> > class PointCloud {
  public:
    PointCloud() : m_points()  { }
    PointCloud(const int &npoints) : m_points() { 
      m_points.reserve(npoints); 
    }

    PointCloud(const QVector<T> &points) : m_points(points) { }
    virtual ~PointCloud() { }


    /** Standard size method */
    inline int size() const {
      return ( m_points.size() );
    }

    /** Add a new point to the list */
    inline void addPoint(const T &point) {
      m_points.push_back(point);
    }

    /** Return a reference to the point at index idx */
    inline const T &point(const size_t idx) const {
      Q_ASSERT( idx >= 0 );
      Q_ASSERT( idx < (size_t) size() );
      return (m_points[idx]);
    }

    inline double distance( const T &first, const T &second ) const {
      return ( m_distance(first, second) );
    }

        /** Return number of points in cloud */
    inline size_t kdtree_get_point_count() const {
      return ( m_points.size() );
    }

  /**
   * @brief Return Euclidean distance from a source point to the indexed point
   *  
   * This method returns the Euclidean (L2) distance from a dataset point and a
   * source point (p1).  We only use X and Y components to compute this 
   * distance since our objective is to identify points in 2-d space. 
   * 
   * @author 2014-02-27 Kris Becker
   * 
   * @param p1      3-vector of the source point
   * @param idx_p2  Index into point data set contained herein
   * 
   * @return double Returns the squared distance - not the square root!
   */
    inline double kdtree_distance(const double *p1, const size_t idx_p2,
                                  size_t p_size) const {
      Q_ASSERT( idx_p2 >= 0 );
      Q_ASSERT( idx_p2 < (size_t) size() );
      return ( m_distance(p1, point(idx_p2)) );
    }

  /**
   * @brief Returns a value for a single dimemsion of a vector 
   *  
   * This method provides a simple interface to each vector element in the point 
   * cloud.
   * 
   * @history 2014-03-03 Kris Becker
   * 
   * @param idx Index of the point to get element from
   * @param dim The index (0 - 2)of the ith element of the vector 
   * 
   * @return double Value at the specfied vector point index
   */
    inline double kdtree_get_pt(const size_t idx, int dim) const {
      if ( dim == 0 )  return ( point(idx).x() );
      if ( dim == 1 )  return ( point(idx).y() );
      return ( point(idx).z() );
    }

    /** 
     *  Let nanoflann range algorithm compute the bounding box
     *  
     * @author 2014-02-17 Kris Becker
     *  
     * @internal 
     *   @history 2014-02-17 Kris Becker - Original Version 
     */
    template <class BBOX> bool kdtree_get_bbox(BBOX &bb) const { 
      return (false); 
    }

  private:
    QVector<T> m_points;   ///!< points in the point cloud
    Distance   m_distance; ///!< Instantiation of distance functor

};

};  // namespace Isis
#endif
+182 −0

File added.

Preview size limit exceeded, changes collapsed.

Loading