version 4.0.0

Typedefine.h → Typedefine.hpp

@@ -14,16 +14,6 @@
  * limitations under the License.
  */
 
-/*
- * Includes code derived from the voronoi algorithm by Steven Fortune
- * (http://ect.bell-labs.com/who/sjf/)
- * as modified by Derek Bradley
- * (http://zurich.disneyresearch.com/derekbradley/voronoi.html)
- *
- * Reference: Steve J. Fortune (1987) A Sweepline Algorithm for Voronoi Diagrams,
- * Algorithmica 2, 153-174.
- */
-
 #include "tiff.h"
 
 #ifndef _Typedefine_H_
@@ -37,7 +27,7 @@
 #define MaxImages 3
 
 #ifndef bool
-#define bool unsigned char
+//#define bool unsigned char
 #define true 0x1
 #define false 0x0
 #endif
@@ -155,7 +145,7 @@ typedef struct UpdateGrid{
 
 typedef struct tagVoxelinfo
 {
-    //float GNCC;
+    //float ANCC;
     float height;
     bool flag_cal;
     float INCC;
@@ -236,6 +226,7 @@ typedef struct ProjectInfo{
 	int threads_num;
     int number_of_images;
     enum SensorType sensor_type; // 1 is for RFM (default), 2 is for Collinear Equation (Frame)
+    uint8 pyramid_level;
 
 	char Imagefilename[MaxImages][500];
 	char RPCfilename[MaxImages][500];
@@ -259,6 +250,18 @@ typedef struct ProjectInfo{
 	bool IsRA, IsSP, IsRR, IsSaveStep, Overall_DEM, Affine_RA, pre_DEMtif, check_tile_array;
     bool check_Matchtag;
     bool check_selected_image[MaxImages];
+    bool check_full_cal;
+
+    //SGM test flag
+    bool check_SNCC;
+    bool check_updateheight;
+    bool check_blunderdetection;
+    bool check_NCCpeak;
+    bool check_minTH;
+    bool check_adaptive_P2;
+    bool check_orthoblunder;
+    bool check_8d;
+    int SGM_py;
 
     FrameInfo frameinfo;
 
@@ -295,7 +298,8 @@ typedef struct ArgumentInfo{
     int ortho_count;
     int RA_only;
     int number_of_images; // 2 is for stereo (default), n is for multi more than 3
-
+    uint8 pyramid_level;
+
     char Image[MaxImages][500];
 	char Outputpath[500];
 	char Outputpath_name[500];
@@ -331,6 +335,18 @@ typedef struct ArgumentInfo{
     bool check_EO;
     bool check_fl;
     bool check_ccd;
+    bool check_full_cal;
+
+    //SGM test flag
+    bool check_SNCC;
+    bool check_updateheight;
+    bool check_blunderdetection;
+    bool check_NCCpeak;
+    bool check_minTH;
+    bool check_adaptive_P2;
+    bool check_orthoblunder;
+    bool check_8d;
+    int SGM_py;
 
 } ARGINFO;
 

basic_topology_types.hpp

@@ -0,0 +1,178 @@
+#ifndef BASIC_TOPOLOGY_TYPES_H
+#define BASIC_TOPOLOGY_TYPES_H
+
+#include <cstddef>
+#include <stdint.h>
+
+/**************************
+ **************************
+ * MACROS, TYPEDEFS,      *
+ * FOREWARD DECLARATIONS, *
+ * ETC.                   *
+ **************************
+ **************************/
+#define THREADED_CUTOFF 100	// Number of points after which the algorithm should switch to a serial triangulation
+
+// TODO
+typedef int64_t INT64;
+//#ifdef __SIZEOF_INT128__ // __int128_t is a compiler extension
+//	typedef __int128_t INT128; 
+//#else // Need to implement own 128-bit int type (or live with 64-bit?)
+	typedef int64_t INT128;
+//#endif
+typedef int INDEX;
+typedef struct GridPoint GridPoint;
+typedef struct Edge Edge;
+typedef struct ExtremeEdges ExtremeEdges;
+
+/*************************************************
+ *************************************************
+ * Primary data structures used in triangulation *
+ *************************************************
+ *************************************************/
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Data structure for point in grid
+struct GridPoint
+{
+	INDEX row; // Coordinates
+	INDEX col;
+	bool operator==(const GridPoint &other) { return this->row == other.row && this->col == other.col; };
+	bool operator!=(const GridPoint &other) { return this->row != other.row || this->col != other.col; };
+};
+
+// Compare function to sort by x/col then y/row
+inline bool LessThanXY(const GridPoint &a, const GridPoint &b) { return (a.col < b.col) || ((a.col == b.col) && (a.row < b.row)); }
+// Compare function to sort by y/row then x/col
+inline bool LessThanYX(const GridPoint &a, const GridPoint &b) { return (a.row < b.row) || ((a.row == b.row) && (a.col > b.col)); } 
+// Conversion functions from GridPoint to linear index and vice versa
+inline GridPoint Convert(std::size_t index, INDEX width) { GridPoint p = { index / width, index % width }; return p; }
+inline std::size_t Convert(const GridPoint &p, INDEX width) { return p.row * width + p.col; }
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Data structure for (half-)edge in grid
+struct Edge
+{
+	Edge *dnext;	// Next edge out of destination vertex, counter-clockwise
+	Edge *oprev;	// Previous edge into origin vertex, counter-clockwise (that is, next edge into origin, clockwise)
+	Edge *twin;	// The companion edge with origin and destination switched
+	GridPoint orig;	// The origin of the edge
+};
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Data structure to hold various extreme edges of a triangulation/convex hull
+struct ExtremeEdges
+{
+	Edge *left_edge_ccw;	// The counter-clockwise edge of the hull with leftmost origin
+	Edge *right_edge_cw;	// The clockwise edge of the hull with rightmost origin
+	Edge *bottom_edge_ccw;	// The counter-clockwise edge of the hull with lowest origin
+	Edge *top_edge_cw;	// The clockwise edge of the hull with highest origin
+};
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// Forward declaration
+template <typename GridType, typename IterType>
+class Grid;
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Data structure to store the (half-)edges of a triangulation
+class EdgeList
+{
+	private:
+		Edge *edges;		// List of Edge objects
+		Edge **unused_edges;	// List of pointers to currently unused positions in 'edges'
+		std::size_t idx;	// Current length of 'unused_edges' 
+		std::size_t size;	// Maximum number of edges to store
+		bool is_part_of_split;	// Denotes whether this is part of a larger, split EdgeList
+
+		// Alternative constructor to create EdgeList from components of a base EdgeList
+		//   edges - list of edge objects from base EdgeList
+		//   unused_edges - list of unused edges from base EdgeList
+		//   num_points = maximum number of points for THIS triangulation (not for base triangulation)
+		EdgeList(Edge *edges, Edge **unused_edges, std::size_t num_points);
+	public:
+		// Standard constructor to create EdgeList from scratch
+		//   num_points - maximum number of points in corresponding triangulation
+		EdgeList(std::size_t num_points);
+		// Deconstructor
+		~EdgeList();
+		// Returns a pointer to a valid Edge object for use in the triangulation
+		Edge *GetNewEdge();
+		// Removes an Edge from triangulation
+		//   edge - the Edge object to be removed from triangulation
+		void RemoveEdge(Edge &edge);
+		// Splits the current EdgeList into two sub lists
+		//   left_list - pointer to an EdgeList; will be replaced with 'left half'
+		//     of the current EdgeList
+		//   right_list - pointer to an EdgeList; will be replaced with 'right half'
+		//     of the current EdgeList
+		//   median - the desired number of points in the triangulation corresponding
+		//     to left_list (leaving the rest for right_list)
+		// The current edge list should be empty (all edges unused) when split for the
+		// SplitEdgeList + MergeEdgeLists combination to work properly
+		void SplitEdgeList(EdgeList *&left_list, EdgeList *&right_list, std::size_t median);
+		// Merges two EdgeLists into one EdgeList (this)
+		//   left_list - the left EdgeList from SplitEdgeList call
+		//   right_list - the right EdgeList from SplitEdgeList call
+		// The input EdgeLists should have been created by calling SplitEdgeList on this
+		void MergeEdgeLists(EdgeList &left_list, EdgeList &right_list);
+
+		// Sets grid so that each grid point corresponding to a point in the
+		// triangulation stores an pointer to an Edge object out of that point
+		//   grid - the grid object for storing Edge pointers
+		// Best to call only once after EdgeList (corresponding GridTriangulation)
+		// has been fully triangulated
+		template <typename GridType, typename IterType>
+		void SetGrid(Grid<GridType, IterType> *grid);
+};
+
+template <typename GridType, typename IterType>
+void EdgeList::SetGrid(Grid<GridType, IterType> *grid)
+{
+	bool *invalid_edge = new bool[this->size];
+	#pragma omp parallel
+	{
+		// Create a companion array that holds
+		// booleans tracking whether or not the
+		// corresponding edge memory is actually
+		// used
+		#pragma omp for
+		for (size_t t = 0; t < this->size; t++)
+		{
+			invalid_edge[t] = false;
+		}
+		#pragma omp for
+		for (size_t t = 0; t < this->idx; t++)
+		{
+			invalid_edge[this->unused_edges[t] - this->edges] = true;
+		}
+
+		// Loop over edges and set the edge out of the
+		// origin of the current edge to the current edge
+		// Threading race conditions are not an issue, since
+		// and succesful write is valid (TODO ?)
+		/*
+		#pragma omp for
+		for (size_t t = 0; t < this->size; t++)
+		{
+			if (invalid_edge[t]) continue;
+			
+			grid->SetElem((this->edges[t]).orig, this->edges + t);
+		}
+		*/
+	}
+
+	for (size_t t = 0; t < this->size; t++)
+	{
+		if (invalid_edge[t]) continue;
+
+		grid->SetElem((this->edges[t]).orig, this->edges + t);
+	}
+
+	delete [] invalid_edge;
+}
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+#endif // BASIC_TOPOLOGY_TYPES_H

edge_list.cpp

@@ -0,0 +1,75 @@
+#include <cstddef>
+#include "basic_topology_types.hpp"
+
+EdgeList::EdgeList(std::size_t num_points)
+{
+	// 3F < 2E, E <= F + V <= (2/3)E + V <-> E <= 3V
+	// We store edge and twin, so need memory of 6*num_points
+	// Allocate all memory upfront and use stack structure
+	// (this->unused_edges + this->idx) to track currently
+	// unused memory
+	this->size = 6 * num_points;
+	this->idx = this->size;
+	this->edges = new Edge[this->size];
+	this->unused_edges = new Edge*[this->size];
+	for (std::size_t t = 0; t < this->size; t++)
+	{
+		this->unused_edges[t] = this->edges + t;
+	}
+	this->is_part_of_split = false;
+}
+
+EdgeList::EdgeList(Edge *edges, Edge **unused_edges, size_t num_points)
+{
+	// Constructor used to make sub-EdgeList
+	// from a larger EdgeList
+	this->size = 6 * num_points;
+	this->idx = this->size;
+	this->edges = edges;
+	this->unused_edges = unused_edges;
+	this->is_part_of_split = true;
+}
+
+EdgeList::~EdgeList()
+{
+	// In case this is not part of a
+	// split EdgeList, free edge memory
+	if (!this->is_part_of_split)
+	{
+		delete [] this->edges;
+		delete [] this->unused_edges;
+	}
+}
+
+Edge *EdgeList::GetNewEdge()
+{
+	this->idx--;
+	return this->unused_edges[this->idx];
+}
+
+void EdgeList::RemoveEdge(Edge &edge)
+{
+	this->unused_edges[this->idx] = &edge;
+	this->idx++;
+}
+
+void EdgeList::SplitEdgeList(EdgeList *&left_list, EdgeList *&right_list, size_t median)
+{
+	// Splits this into two sub EdgeLists with enough memory to store
+	// edges for a triangulation of median and (this->size / 6) - median
+	// points, respectively
+	left_list = new EdgeList(this->edges, this->unused_edges, median);
+	right_list = new EdgeList(this->edges + 6 * median, this->unused_edges + 6 * median, (this->size / 6) - median);
+}
+
+void EdgeList::MergeEdgeLists(EdgeList &left_list, EdgeList &right_list)
+{
+	// Combine lists of unused edges, updating
+	// this->idx appropriately
+	this->idx = left_list.idx;
+	for (size_t t = 0; t < right_list.idx; t++)
+	{
+		this->unused_edges[this->idx + t] = right_list.unused_edges[t];
+	}
+	this->idx += right_list.idx;
+}

grid.cpp

@@ -0,0 +1,53 @@
+#include "grid_types.hpp"
+
+template <typename GridType, typename IterType>
+Grid<GridType, IterType>::Grid(INDEX width, INDEX height, GridPoint *points[], std::size_t num_points)
+{
+	this->grid = new GridType(width, height, points, num_points);
+}
+
+template <typename GridType, typename IterType>
+Grid<GridType, IterType>::~Grid()
+{
+	delete this->grid;
+}
+
+template <typename GridType, typename IterType>
+std::size_t Grid<GridType, IterType>::GetAllTris(GridPoint tris[][3])
+{
+	std::size_t idx = 0;
+
+	// Loop over all points in triangulation, collecting adjacent triangles
+	// To prevent 'over-counting' only collect triangle from its smallest vertex
+	for (auto p_it = this->PointBegin(); p_it != this->PointEnd(); p_it++)
+	{
+		Edge *e = this->GetElem(*p_it);
+		Edge *f = e->twin->dnext->twin->dnext; // If degree of current point is 1 or 2, then f == e
+
+		GridPoint p1 = e->orig;
+		GridPoint p2 = e->twin->orig;
+		GridPoint p3 = e->twin->dnext->twin->orig;
+
+		// Check that degree > 2 or
+		//  at least edge(s) out not parallel
+		if ((f != e) || !((p1.col - p2.col)*(p1.row - p3.row) == (p1.row - p2.row)*(p1.col - p3.col)))
+		{
+			for (auto t_it = this->AdjTriBegin(*p_it); t_it != this->AdjTriEnd(*p_it); t_it++)
+			{
+				GridPoint *points = *t_it;
+				if (LessThanXY(points[0], points[1]) && LessThanXY(points[0], points[2]))
+				{
+					tris[idx][0] = points[0];
+					tris[idx][1] = points[1];
+					tris[idx][2] = points[2];
+					idx++;
+				}
+			}
+		}
+	}
+
+	return idx;
+}
+
+template class Grid<SparseGrid, SparsePointIter>;
+template class Grid<FullGrid, FullPointIter>;