Revision: 77315
http://sourceforge.net/p/brlcad/code/77315
Author: starseeker
Date: 2020-10-01 21:41:50 +0000 (Thu, 01 Oct 2020)
Log Message:
-----------
gdiam is another easy one, without any particular upstream activity.
Modified Paths:
--------------
brlcad/trunk/INSTALL
brlcad/trunk/configure
brlcad/trunk/doc/legal/embedded/CMakeLists.txt
brlcad/trunk/doc/legal/other/CMakeLists.txt
brlcad/trunk/misc/win32-msvc/Dll/CMakeLists.txt
brlcad/trunk/src/external/Creo/CMakeLists.txt
brlcad/trunk/src/external/Cubit/CMakeLists.txt
brlcad/trunk/src/external/Unigraphics/CMakeLists.txt
brlcad/trunk/src/librt/CMakeLists.txt
brlcad/trunk/src/other/CMakeLists.txt
Added Paths:
-----------
brlcad/trunk/doc/legal/embedded/gdiam.txt
brlcad/trunk/src/librt/gdiam/
brlcad/trunk/src/librt/gdiam/LICENSE.MIT
brlcad/trunk/src/librt/gdiam/README
brlcad/trunk/src/librt/gdiam/gdiam.cpp
brlcad/trunk/src/librt/gdiam/gdiam.hpp
brlcad/trunk/src/librt/gdiam/gdiam_test.cpp
Removed Paths:
-------------
brlcad/trunk/doc/legal/other/libgdiam.txt
brlcad/trunk/src/other/libgdiam/
brlcad/trunk/src/other/libgdiam.dist
Modified: brlcad/trunk/INSTALL
===================================================================
--- brlcad/trunk/INSTALL 2020-10-01 21:19:02 UTC (rev 77314)
+++ brlcad/trunk/INSTALL 2020-10-01 21:41:50 UTC (rev 77315)
@@ -681,17 +681,6 @@
Aliases: ENABLE_OPENSCENEGRAPH
---- BRLCAD_GDIAM ---
-
-Option for enabling and disabling compilation of the libgdiam approximate
-tight bounding box library provided with BRL-CAD's source code. Default
-is AUTO, responsive to the toplevel BRLCAD_BUNDLED_LIBS option and
-testing first for a system version if BRLCAD_BUNDLED_LIBS is also
-AUTO.
-
-Aliases: ENABLE_GDIAM
-
-
--- BRLCAD_PROJ4 ---
Option for enabling and disabling compilation of the PROJ.4 geographic
Modified: brlcad/trunk/configure
===================================================================
--- brlcad/trunk/configure 2020-10-01 21:19:02 UTC (rev 77314)
+++ brlcad/trunk/configure 2020-10-01 21:41:50 UTC (rev 77315)
@@ -178,10 +178,6 @@
shift;;
--disable-openscenegraph) options="$options
-DBRLCAD_OSG=SYSTEM";
shift;;
- --enable-gdiam) options="$options -DBRLCAD_GDIAM=BUNDLED";
- shift;;
- --disable-gdiam) options="$options -DBRLCAD_GDIAM=SYSTEM";
- shift;;
--enable-proj4) options="$options -DBRLCAD_PROJ4=BUNDLED";
shift;;
--disable-proj4) options="$options -DBRLCAD_PROJ4=SYSTEM";
Modified: brlcad/trunk/doc/legal/embedded/CMakeLists.txt
===================================================================
--- brlcad/trunk/doc/legal/embedded/CMakeLists.txt 2020-10-01 21:19:02 UTC
(rev 77314)
+++ brlcad/trunk/doc/legal/embedded/CMakeLists.txt 2020-10-01 21:41:50 UTC
(rev 77315)
@@ -23,6 +23,7 @@
gecode.txt
gdal_gcv_plugin1.txt
gdal_gcv_plugin2.txt
+ gdiam.txt
halfedge.txt
humanize.txt
hv3.txt
Copied: brlcad/trunk/doc/legal/embedded/gdiam.txt (from rev 77314,
brlcad/trunk/doc/legal/other/libgdiam.txt)
===================================================================
--- brlcad/trunk/doc/legal/embedded/gdiam.txt (rev 0)
+++ brlcad/trunk/doc/legal/embedded/gdiam.txt 2020-10-01 21:41:50 UTC (rev
77315)
@@ -0,0 +1,23 @@
+Copyright 2001 Sariel Har-Peled ([email protected])
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+file:/src/librt/gdiam/gdiam.cpp
+file:/src/librt/gdiam/gdiam.hpp
+
Modified: brlcad/trunk/doc/legal/other/CMakeLists.txt
===================================================================
--- brlcad/trunk/doc/legal/other/CMakeLists.txt 2020-10-01 21:19:02 UTC (rev
77314)
+++ brlcad/trunk/doc/legal/other/CMakeLists.txt 2020-10-01 21:41:50 UTC (rev
77315)
@@ -4,7 +4,6 @@
incrTcl.txt
libbson.txt
libbson_yajl.txt
- libgdiam.txt
libnetpbm.txt
libpng.txt
libregex.txt
Deleted: brlcad/trunk/doc/legal/other/libgdiam.txt
===================================================================
--- brlcad/trunk/doc/legal/other/libgdiam.txt 2020-10-01 21:19:02 UTC (rev
77314)
+++ brlcad/trunk/doc/legal/other/libgdiam.txt 2020-10-01 21:41:50 UTC (rev
77315)
@@ -1,20 +0,0 @@
-Copyright 2001 Sariel Har-Peled ([email protected])
-
-Permission is hereby granted, free of charge, to any person obtaining a copy of
-this software and associated documentation files (the "Software"), to deal in
-the Software without restriction, including without limitation the rights to
-use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-of the Software, and to permit persons to whom the Software is furnished to do
-so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
-
Modified: brlcad/trunk/misc/win32-msvc/Dll/CMakeLists.txt
===================================================================
--- brlcad/trunk/misc/win32-msvc/Dll/CMakeLists.txt 2020-10-01 21:19:02 UTC
(rev 77314)
+++ brlcad/trunk/misc/win32-msvc/Dll/CMakeLists.txt 2020-10-01 21:41:50 UTC
(rev 77315)
@@ -38,7 +38,6 @@
libbn-static
libbrep-static
libbu-static
- gdiam-static
libged-static
regex-static
librt-static
@@ -118,7 +117,6 @@
libbn-static
libbrep-static
libbu-static
- gdiam-static
libged-static
regex-static
librt-static
Modified: brlcad/trunk/src/external/Creo/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/external/Creo/CMakeLists.txt 2020-10-01 21:19:02 UTC
(rev 77314)
+++ brlcad/trunk/src/external/Creo/CMakeLists.txt 2020-10-01 21:41:50 UTC
(rev 77315)
@@ -38,7 +38,7 @@
set(CREO_DEFS "-DPRO_MACHINE=36 -DPRO_OS=4 -DPRO_USE_VAR_ARGS
-D_USING_V100_SDK71_")
# BRL-CAD definitions
- set(BRLCAD_DEFS "-DHAVE_CONFIG_H -DBRLCAD_DLL -DBRLCADBUILD -DBU_DLL_IMPORTS
-DBN_DLL_IMPORTS -DNMG_DLL_IMPORTS -DBG_DLL_IMPORTS -DBREP_DLL_IMPORTS
-DRT_DLL_IMPORTS -DWDB_DLL_IMPORTS -DTIE_DLL_IMPORTS -DDB5_DLL_IMPORTS
-DGDIAM_DLL_IMPORTS -DON_DLL_IMPORTS")
+ set(BRLCAD_DEFS "-DHAVE_CONFIG_H -DBRLCAD_DLL -DBRLCADBUILD -DBU_DLL_IMPORTS
-DBN_DLL_IMPORTS -DNMG_DLL_IMPORTS -DBG_DLL_IMPORTS -DBREP_DLL_IMPORTS
-DRT_DLL_IMPORTS -DWDB_DLL_IMPORTS -DTIE_DLL_IMPORTS -DDB5_DLL_IMPORTS
-DON_DLL_IMPORTS")
# These settings are global to all the configs.
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${CREO_DEFS} ${BRLCAD_DEFS}" CACHE
STRING "" FORCE)
@@ -207,7 +207,6 @@
libnmg
librt
libwdb
- gdiam
openNURBS
poly2tri
regex_brl
@@ -369,7 +368,6 @@
set_property(TARGET creo-brl APPEND PROPERTY COMPILE_DEFINITIONS
"DB5_DLL_IMPORTS")
set_property(TARGET creo-brl APPEND PROPERTY COMPILE_DEFINITIONS
"WDB_DLL_IMPORTS")
set_property(TARGET creo-brl APPEND PROPERTY COMPILE_DEFINITIONS
"TIE_DLL_IMPORTS")
- set_property(TARGET creo-brl APPEND PROPERTY COMPILE_DEFINITIONS
"GDIAM_DLL_IMPORTS")
set_property(TARGET creo-brl APPEND PROPERTY COMPILE_DEFINITIONS
"ON_DLL_IMPORTS")
endif(HIDE_INTERNAL_SYMBOLS)
Modified: brlcad/trunk/src/external/Cubit/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/external/Cubit/CMakeLists.txt 2020-10-01 21:19:02 UTC
(rev 77314)
+++ brlcad/trunk/src/external/Cubit/CMakeLists.txt 2020-10-01 21:41:50 UTC
(rev 77315)
@@ -35,7 +35,6 @@
set_property(TARGET g-sat APPEND PROPERTY COMPILE_DEFINITIONS
"DB5_DLL_IMPORTS")
set_property(TARGET g-sat APPEND PROPERTY COMPILE_DEFINITIONS
"WDB_DLL_IMPORTS")
set_property(TARGET g-sat APPEND PROPERTY COMPILE_DEFINITIONS
"TIE_DLL_IMPORTS")
- set_property(TARGET g-sat APPEND PROPERTY COMPILE_DEFINITIONS
"GDIAM_DLL_IMPORTS")
set_property(TARGET g-sat APPEND PROPERTY COMPILE_DEFINITIONS
"ON_DLL_IMPORTS")
set_property(TARGET g-sat APPEND PROPERTY COMPILE_DEFINITIONS
"TINYCTHREAD_DLL_IMPORTS")
endif(HIDE_INTERNAL_SYMBOLS)
Modified: brlcad/trunk/src/external/Unigraphics/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/external/Unigraphics/CMakeLists.txt 2020-10-01
21:19:02 UTC (rev 77314)
+++ brlcad/trunk/src/external/Unigraphics/CMakeLists.txt 2020-10-01
21:41:50 UTC (rev 77315)
@@ -40,10 +40,8 @@
set_property(TARGET ug-g APPEND PROPERTY COMPILE_DEFINITIONS
"DB5_DLL_IMPORTS")
set_property(TARGET ug-g APPEND PROPERTY COMPILE_DEFINITIONS
"WDB_DLL_IMPORTS")
set_property(TARGET ug-g APPEND PROPERTY COMPILE_DEFINITIONS
"TIE_DLL_IMPORTS")
- set_property(TARGET ug-g APPEND PROPERTY COMPILE_DEFINITIONS
"GDIAM_DLL_IMPORTS")
set_property(TARGET ug-g APPEND PROPERTY COMPILE_DEFINITIONS
"ON_DLL_IMPORTS")
set_property(TARGET ug-g APPEND PROPERTY COMPILE_DEFINITIONS
"TINYCTHREAD_DLL_IMPORTS")
- set_property(TARGET ug-g APPEND PROPERTY COMPILE_DEFINITIONS
"LZ4_DLL_IMPORT=1")
endif(HIDE_INTERNAL_SYMBOLS)
endif (BRLCAD_ENABLE_TCL)
Modified: brlcad/trunk/src/librt/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/librt/CMakeLists.txt 2020-10-01 21:19:02 UTC (rev
77314)
+++ brlcad/trunk/src/librt/CMakeLists.txt 2020-10-01 21:41:50 UTC (rev
77315)
@@ -20,8 +20,8 @@
set(RT_LOCAL_INCLUDE_DIRS
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/vds
+ ${CMAKE_CURRENT_SOURCE_DIR}/gdiam
${REGEX_INCLUDE_DIRS}
- ${GDIAM_INCLUDE_DIRS}
)
BRLCAD_LIB_INCLUDE_DIRS(rt RT_INCLUDE_DIRS RT_LOCAL_INCLUDE_DIRS)
@@ -68,6 +68,7 @@
dir.c
dspline.c
fortray.c
+ gdiam/gdiam.cpp
globals.c
htbl.c
ls.c
@@ -363,21 +364,15 @@
set(SPR_LIB libSPR)
endif(BRLCAD_ENABLE_SPR)
-BRLCAD_ADDLIB(librt "${LIBRT_SOURCES}"
"${OPENCL_LIBS};${GDIAM_LIBRARY};libbrep;libnmg;libbg;libbn;libbu;${OPENNURBS_LIBRARIES};${REGEX_LIBRARIES};${WINSOCK_LIB};${RPCRT_LIB};${STDCXX_LIBRARIES}")
+BRLCAD_ADDLIB(librt "${LIBRT_SOURCES}"
"${OPENCL_LIBS};libbrep;libnmg;libbg;libbn;libbu;${OPENNURBS_LIBRARIES};${REGEX_LIBRARIES};${WINSOCK_LIB};${RPCRT_LIB};${STDCXX_LIBRARIES}")
set_target_properties(librt PROPERTIES VERSION 20.0.1 SOVERSION 20)
if(HIDE_INTERNAL_SYMBOLS)
set_property(TARGET librt APPEND PROPERTY COMPILE_DEFINITIONS
"TIE_DLL_EXPORTS")
set_property(TARGET librt APPEND PROPERTY COMPILE_DEFINITIONS
"DB5_DLL_EXPORTS")
- if(TARGET gdiam OR HIDE_INTERNAL_SYMBOLS_EXT)
- set_property(TARGET librt APPEND PROPERTY COMPILE_DEFINITIONS
"GDIAM_DLL_IMPORTS")
- endif(TARGET gdiam OR HIDE_INTERNAL_SYMBOLS_EXT)
if (TARGET librt-obj)
set_property(TARGET librt-obj APPEND PROPERTY COMPILE_DEFINITIONS
"TIE_DLL_EXPORTS")
set_property(TARGET librt-obj APPEND PROPERTY COMPILE_DEFINITIONS
"DB5_DLL_EXPORTS")
- if(TARGET gdiam OR HIDE_INTERNAL_SYMBOLS_EXT)
- set_property(TARGET librt-obj APPEND PROPERTY COMPILE_DEFINITIONS
"GDIAM_DLL_IMPORTS")
- endif(TARGET gdiam OR HIDE_INTERNAL_SYMBOLS_EXT)
endif (TARGET librt-obj)
endif(HIDE_INTERNAL_SYMBOLS)
@@ -403,6 +398,10 @@
vds/README
vds/COPYING
vds/vds.h
+ gdiam/gdiam.hpp
+ gdiam/gdiam_test.cpp
+ gdiam/LICENSE.MIT
+ gdiam/README
)
# Local Variables:
Copied: brlcad/trunk/src/librt/gdiam/LICENSE.MIT (from rev 77314,
brlcad/trunk/src/other/libgdiam/LICENSE.MIT)
===================================================================
--- brlcad/trunk/src/librt/gdiam/LICENSE.MIT (rev 0)
+++ brlcad/trunk/src/librt/gdiam/LICENSE.MIT 2020-10-01 21:41:50 UTC (rev
77315)
@@ -0,0 +1,20 @@
+Copyright 2001 Sariel Har-Peled ([email protected])
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
Copied: brlcad/trunk/src/librt/gdiam/README (from rev 77314,
brlcad/trunk/src/other/libgdiam/README)
===================================================================
--- brlcad/trunk/src/librt/gdiam/README (rev 0)
+++ brlcad/trunk/src/librt/gdiam/README 2020-10-01 21:41:50 UTC (rev 77315)
@@ -0,0 +1,61 @@
+Source code implementing the algorithms described in:
+
+A Practical Approach for Computing the Diameter of a Point-Set
+Sariel Har-Peled
+
+Copyright 2001 Sariel Har-Peled ([email protected])
+http://valis.cs.uiuc.edu/~sariel/research/papers/00/diameter/diam_prog.html
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of either:
+
+ * the GNU General Public License as published by the Free Software
+ Foundation; either version 2, or (at your option) any later version.
+
+or
+
+ * the GNU Lesser General Public License as published by the Free
+ Software Foundation; either version 2.1, or (at your option)
+ any later version.
+
+or
+
+ * MIT License https://opensource.org/licenses/MIT
+
+If you need the source code under a different open-source license, let
+me know.
+
+Program is provdied without any guarantee. Use at your own risk.
+
+*--------------------------------------------------------------
+ * History
+ *
+ * 3/6/18
+ * - Tobias Stohr reported & fixed a bug wis missing
+ * constructor for
+ *
+ * 8/18/16 -
+ *
+ * Apply various code tweaks from BRL-CAD (C. Yapp)
+ *
+ * 12/22/14 -
+ *
+ * Apply GSoC patch switching the convex hull algorithm to
+ * monotone chain (P. Amidon)
+ *
+ * 8/7/13 -
+ *
+ * Add DLL import/export logic for Windows (C. Yapp)
+ *
+ * 8/4/13 -
+ * Add get_vertex method for low-level data translation. (C. Yapp)
+ *
+ * 8/3/13 -
+ * Update licensing - can now use either GPLv2 or LGPLv2.1.
+ * Added CMake build.
+ *
+ * 3/28/01 -
+ * Original code updated to be more robust. It should now
+ * handle really abnoxious inputs well (i.e., points with equal
+ * coordinates, etc.
+
Copied: brlcad/trunk/src/librt/gdiam/gdiam.cpp (from rev 77314,
brlcad/trunk/src/other/libgdiam/gdiam.cpp)
===================================================================
--- brlcad/trunk/src/librt/gdiam/gdiam.cpp (rev 0)
+++ brlcad/trunk/src/librt/gdiam/gdiam.cpp 2020-10-01 21:41:50 UTC (rev
77315)
@@ -0,0 +1,2402 @@
+/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
+ * gdiam.cpp -
+ * Computes diameter, and a tight-fitting bounding box of a
+ * point-set in 3d.
+ *
+ * Copyright 2001 Sariel Har-Peled ([email protected])
+ * https://sarielhp.org/research/papers/00/diameter/diam_prog.html
+ *
+ * Note: gdiam is available under multiple licenses - we make
+ * use of it under the MIT license:
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
copy of
+ * this software and associated documentation files (the "Software"), to deal
in
+ * the Software without restriction, including without limitation the rights to
+ * use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies
+ * of the Software, and to permit persons to whom the Software is furnished to
do
+ * so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE
+ * SOFTWARE.
+ *
+ * Code is based on the paper:
+ * A Practical Approach for Computing the Diameter of a
+ * Point-Set (in ACM Sym. on Computation Geometry SoCG 2001)
+ * Sariel Har-Peled (http://www.uiuc.edu/~sariel)
+ *--------------------------------------------------------------
+ * History
+ * 3/6/18
+ * - Tobias Stohr reported & fixed a bug was missing
+ * constructor for
+ * 8/19/16 - Clifford Yapp updated the source
+ * 3/28/01 -
+ * This is a more robust version of the code. It should
+ * handlereally abnoxious inputs well (i.e., points with equal
+ * coordinates, etc.
+\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/
+
+#include "common.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <memory.h>
+#include <math.h>
+
+#include <iostream>
+#include <vector>
+#include <algorithm>
+
+#include "gdiam.hpp"
+
+/* for g++ to quell warnings */
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic push /* start new diagnostic pragma */
+# pragma GCC diagnostic ignored "-Wfloat-equal"
+#elif defined(__clang__)
+# pragma clang diagnostic push /* start new diagnostic pragma */
+# pragma clang diagnostic ignored "-Wfloat-equal"
+#endif
+
+/*--- Constants ---*/
+
+#define HEAP_LIMIT 10000
+#define HEAP_DELTA 10000
+#define PI 3.14159265958979323846
+
+/*--- Start of Code ---*/
+
+typedef long double ldouble;
+
+class GFSPPair;
+
+class GFSPTreeNode {
+private:
+ GBBox bbox;
+ gdiam_point * p_pnt_left, * p_pnt_right;
+ gdiam_real bbox_diam, bbox_diam_proj;
+ GFSPTreeNode * left, * right;
+ gdiam_point_t center;
+
+public:
+ void dump() const {
+ printf( "dump...\n" );
+ }
+ int size() const {
+ return (int)( p_pnt_right - p_pnt_left );
+ }
+ gdiam_point getCenter() const {
+ return (gdiam_point)center;
+ }
+ int nodes_number() const {
+ int sum;
+
+ if ( ( left == NULL ) && ( right == NULL ) )
+ return 1;
+ sum = 1;
+
+ if ( left != NULL )
+ sum += left->nodes_number();
+ if ( right != NULL )
+ sum += right->nodes_number();
+
+ return sum;
+ }
+
+ GFSPTreeNode( gdiam_point * _p_pnt_left,
+ gdiam_point * _p_pnt_right ) {
+ bbox.init();
+ left = right = NULL;
+ p_pnt_left = _p_pnt_left;
+ p_pnt_right = _p_pnt_right;
+ bbox_diam = bbox_diam_proj = -1;
+ }
+
+ gdiam_real maxDiam() const {
+ return bbox_diam;
+ }
+ gdiam_real maxDiamProj() const {
+ return bbox_diam_proj;
+ }
+
+ GFSPTreeNode * get_left() {
+ return left;
+ }
+ GFSPTreeNode * get_right() {
+ return right;
+ }
+
+ bool isSplitable() const {
+ return (size() > 1);
+ }
+
+ const gdiam_point * ptr_pnt_left() {
+ return p_pnt_left;
+ }
+ const gdiam_point * ptr_pnt_rand() {
+ return p_pnt_left + ( rand() % ( p_pnt_right - p_pnt_left + 1 ) );
+ }
+ const gdiam_point * ptr_pnt_right() {
+ return p_pnt_right;
+ }
+
+
+ friend class GFSPTree;
+
+ const GBBox & getBBox() const {
+ return bbox;
+ }
+};
+
+
+class GFSPTree {
+private:
+ gdiam_point * arr;
+ GFSPTreeNode * root;
+
+public:
+ void init( const gdiam_point * _arr,
+ int size ) {
+ arr = (gdiam_point *)malloc( sizeof( gdiam_point ) * size );
+ assert( arr != NULL );
+
+ for ( int ind = 0; ind < size; ind++ ) {
+ arr[ ind ] = _arr[ ind ];
+ }
+
+ root = build_node( arr, arr + size - 1 );
+ }
+
+ static void terminate( GFSPTreeNode * node ) {
+ if ( node == NULL )
+ return;
+ terminate( node->left );
+ terminate( node->right );
+
+ node->left = node->right = NULL;
+
+ delete node;
+ }
+ void term() {
+ free( arr );
+ arr = NULL;
+
+ GFSPTree::terminate( root );
+ root = NULL;
+ }
+
+ const gdiam_point * getPoints() const {
+ return arr;
+ }
+
+ int nodes_number() const {
+ return root->nodes_number();
+ }
+
+ gdiam_real brute_diameter( int a_lo, int a_hi,
+ int b_lo, int b_hi,
+ GPointPair & diam ) const {
+ for ( int ind = a_lo; ind <= a_hi; ind++ )
+ for ( int jnd = b_lo; jnd <= b_hi; jnd++ )
+ diam.update_diam( arr[ ind ], arr[ jnd ] );
+
+ return diam.distance;
+ }
+ gdiam_real brute_diameter( int a_lo, int a_hi,
+ int b_lo, int b_hi,
+ GPointPair & diam,
+ const gdiam_point dir ) const {
+ for ( int ind = a_lo; ind <= a_hi; ind++ )
+ for ( int jnd = b_lo; jnd <= b_hi; jnd++ )
+ diam.update_diam( arr[ ind ], arr[ jnd ], dir );
+
+ return diam.distance;
+ }
+ gdiam_real brute_diameter( const gdiam_point * a_lo,
+ const gdiam_point * a_hi,
+ const gdiam_point * b_lo,
+ const gdiam_point * b_hi,
+ GPointPair & diam ) const {
+ for ( const gdiam_point * ind = a_lo; ind <= a_hi; ind++ )
+ for ( const gdiam_point * jnd = b_lo; jnd <= b_hi; jnd++ )
+ diam.update_diam( *ind, *jnd );
+
+ return diam.distance;
+ }
+ gdiam_real brute_diameter( const gdiam_point * a_lo,
+ const gdiam_point * a_hi,
+ const gdiam_point * b_lo,
+ const gdiam_point * b_hi,
+ GPointPair & diam,
+ const gdiam_point dir ) const {
+ for ( const gdiam_point * ind = a_lo; ind <= a_hi; ind++ )
+ for ( const gdiam_point * jnd = b_lo; jnd <= b_hi; jnd++ )
+ diam.update_diam( *ind, *jnd, dir );
+
+ return diam.distance;
+ }
+
+ gdiam_point getPoint( int ind ) const {
+ return arr[ ind ];
+ }
+
+ GFSPTreeNode * getRoot() {
+ return root;
+ }
+
+ GFSPTreeNode * build_node( gdiam_point * left,
+ gdiam_point * right ) {
+ if ( left > right ) {
+ printf( "what!?\n" );
+ fflush( stdout );
+ assert( left <= right );
+ }
+ while ( ( right > left )
+ && ( pnt_isEqual( *right, *left ) ) )
+ right--;
+ GFSPTreeNode * node = new GFSPTreeNode( left, right );
+
+ node->bbox.init();
+ for ( const gdiam_point * ind = left; ind <= right; ind++ )
+ node->bbox.bound( *ind );
+ node->bbox_diam = node->bbox.get_diam();
+
+ node->bbox.center( node->center );
+
+ return node;
+ }
+
+ // return how many elements on left side.
+ int split_array( gdiam_point * left,
+ gdiam_point * right,
+ int dim, const gdiam_real & val ) {
+ const gdiam_point * start_left = left;
+
+ while ( left < right ) {
+ if ( (*left)[ dim ] < val )
+ left++;
+ else
+ if ( (*right)[ dim ] >= val )
+ right--;
+ else {
+ gdiam_exchange( *right, *left );
+ }
+ }
+
+ return left - start_left;
+ }
+
+
+ void split_node( GFSPTreeNode * node )
+ {
+ int dim, left_size;
+ gdiam_real split_val;
+
+ if ( node->left != NULL )
+ return;
+
+ GBBox & bb( node->bbox );
+
+ dim = bb.getLongestDim();
+ split_val = ( bb.min_coord( dim ) + bb.max_coord( dim ) ) / 2.0;
+
+ left_size = split_array( node->p_pnt_left, node->p_pnt_right,
+ dim, split_val );
+ if ( left_size <= 0.0 ) {
+ printf( "bb: %g %g\n",
+ bb.min_coord( dim ), bb.max_coord( dim ) );
+ assert( left_size > 0 );
+ }
+ if ( left_size >= (node->p_pnt_right - node->p_pnt_left + 1 ) ) {
+ printf( "left size too large?\n" );
+ fflush( stdout );
+ assert( left_size < (node->p_pnt_right - node->p_pnt_left + 1 ) );
+ }
+
+ node->left = build_node( node->p_pnt_left,
+ node->p_pnt_left + left_size - 1 );
+ node->right = build_node( node->p_pnt_left + left_size,
+ node->p_pnt_right );
+ }
+
+ void findExtremeInProjection( GFSPPair & pair,
+ GPointPair & max_pair_diam );
+};
+
+void bbox_vertex( const GBBox & bb, gdiam_point_t & ver,
+ int vert ) {
+ ver[ 0 ] = ((vert & 0x1) != 0)?
+ bb.min_coord( 0 ) : bb.max_coord( 0 );
+ ver[ 1 ] = ((vert & 0x2) != 0)?
+ bb.min_coord( 1 ) : bb.max_coord( 1 );
+ ver[ 2 ] = ((vert & 0x4) != 0)?
+ bb.min_coord( 2 ) : bb.max_coord( 2 );
+}
+
+
+gdiam_real bbox_proj_dist( const GBBox & bb1,
+ const GBBox & bb2,
+ gdiam_point_cnt proj )
+{
+ gdiam_point_t a, b;
+ gdiam_real rl;
+
+ rl = 0;
+ for ( int ind = 0; ind < 8; ind++ ) {
+ bbox_vertex( bb1, a, ind );
+
+ for ( int jnd = 0; jnd < 8; jnd++ ) {
+ bbox_vertex( bb2, b, jnd );
+ rl = max( rl, pnt_distance( a, b, proj ) );
+ }
+ }
+
+ return rl;
+}
+
+
+
+class GFSPPair
+{
+private:
+ GFSPTreeNode * left, * right;
+ double max_diam;
+
+public:
+ void init( GFSPTreeNode * _left, GFSPTreeNode * _right )
+ {
+ left = _left;
+ right = _right;
+
+ GBBox bb;
+
+ bb.init( left->getBBox(), right->getBBox() );
+ max_diam = bb.get_diam();
+ }
+
+ void init( GFSPTreeNode * _left, GFSPTreeNode * _right,
+ gdiam_point proj_dir, gdiam_real UNUSED(dist) )
+ {
+ left = _left;
+ right = _right;
+
+ GBBox bb;
+
+ bb.init( left->getBBox(), right->getBBox() );
+ max_diam = max( max( bbox_proj_dist( _left->getBBox(),
+ _right->getBBox(),
+ proj_dir ),
+ bbox_proj_dist( _left->getBBox(),
+ _left->getBBox(),
+ proj_dir ) ),
+ bbox_proj_dist( _right->getBBox(),
+ _right->getBBox(),
+ proj_dir ) );
+
+ }
+
+ // error 2r/l - approximate cos of the max possible angle in the pair
+ double getProjectionError() const
+ {
+ double l, two_r;
+
+ l = pnt_distance( left->getCenter(), right->getCenter() );
+ two_r = max( left->maxDiam(), right->maxDiam() );
+
+ if ( l == 0.0 )
+ return 10;
+
+ return two_r / l;
+ }
+
+ GFSPTreeNode * get_left() {
+ return left;
+ }
+ GFSPTreeNode * get_right() {
+ return right;
+ }
+
+ void dump() const {
+ printf( "[\n" );
+ left->dump();
+ right->dump();
+ printf( "\tmax_diam; %g\n", max_diam );
+ }
+
+ bool isBelowThreshold( int threshold ) const {
+ if ( left->size() > threshold )
+ return false;
+ if ( right->size() > threshold )
+ return false;
+
+ return true;
+ }
+
+ double directDiameter( const GFSPTree & tree,
+ GPointPair & diam ) const {
+ return tree.brute_diameter( left->ptr_pnt_left(),
+ left->ptr_pnt_right(),
+ right->ptr_pnt_left(),
+ right->ptr_pnt_right(),
+ diam );
+ }
+
+ double directDiameter( const GFSPTree & tree,
+ GPointPair & diam,
+ const gdiam_point dir ) const {
+ return tree.brute_diameter( left->ptr_pnt_left(),
+ left->ptr_pnt_right(),
+ right->ptr_pnt_left(),
+ right->ptr_pnt_right(),
+ diam, dir );
+ }
+
+ double maxDiam() const {
+ return max_diam;
+ }
+ double minAprxDiam() const {
+ double sub_diam;
+
+ sub_diam = left->maxDiam() + right->maxDiam();
+ if ( sub_diam > (max_diam / 10.0) )
+ return max_diam;
+
+ return max_diam - sub_diam / 2.0;
+ }
+};
+
+
+class g_heap_pairs_p;
+
+#define UDM_SIMPLE 1
+#define UDM_BIG_SAMPLE 2
+
+class GTreeDiamAlg
+{
+private:
+ GFSPTree tree;
+ //double diam;
+ GPointPair pair_diam;
+ int points_num;
+ int threshold_brute;
+ int update_diam_mode;
+public:
+ GTreeDiamAlg( gdiam_point * arr, int size, int mode ) {
+ tree.init( arr, size );
+ //diam = 0;
+ pair_diam.init( arr[ 0 ] );
+ points_num = size;
+ threshold_brute = 40;
+ update_diam_mode = mode;
+ }
+
+ void term() {
+ tree.term();
+ }
+
+ int size() const {
+ return points_num;
+ }
+ const GPointPair & getDiameter() const {
+ return pair_diam;
+ }
+
+ int nodes_number() const {
+ return tree.nodes_number();
+ }
+
+ void addPairHeap( g_heap_pairs_p & heap,
+ GFSPTreeNode * left,
+ GFSPTreeNode * right );
+ void addPairHeap( g_heap_pairs_p & heap,
+ GFSPTreeNode * left,
+ GFSPTreeNode * right,
+ gdiam_point proj,
+ GFSPPair & father );
+ void split_pair( GFSPPair & pair,
+ g_heap_pairs_p & heap, double eps );
+ void split_pair_proj( GFSPPair & pair,
+ g_heap_pairs_p & heap, double eps,
+ gdiam_point proj );
+ void compute_by_heap( double eps );
+ void compute_by_heap_proj( double eps, gdiam_point proj );
+
+ double diameter() {
+ return pair_diam.distance;
+ }
+};
+
+
+/*--- Heap implementation ---*/
+
+typedef int ( *ptrCompareFunc )( void * aPtr, void * bPtr );
+typedef void * voidPtr_t;
+
+struct heap_t
+{
+ voidPtr_t * pArr;
+ int curr_size, max_size;
+ ptrCompareFunc pCompFunc;
+};
+
+void heap_init( heap_t * pHeap, ptrCompareFunc _pCompFunc )
+{
+ assert( pHeap != NULL );
+ assert( _pCompFunc != NULL );
+
+ memset( pHeap, 0, sizeof( heap_t ) );
+
+ pHeap->pCompFunc = _pCompFunc;
+ pHeap->max_size = 100;
+ pHeap->pArr = (voidPtr_t *)malloc( sizeof( void * )
+ * pHeap->max_size );
+ assert( pHeap->pArr != NULL );
+ pHeap->curr_size = 0;
+}
+
+
+static void resize( heap_t * pHeap, int size )
+{
+ int max_sz;
+ voidPtr_t * pTmp;
+
+ if ( size <= pHeap->max_size )
+ return;
+
+ max_sz = size * 2;
+ pTmp = (voidPtr_t *)malloc( max_sz * sizeof( void * ) );
+ assert( pTmp != NULL );
+ memset( pTmp, 0, max_sz * sizeof( void * ) );
+ memcpy( pTmp, pHeap->pArr, pHeap->curr_size * sizeof( void * ) );
+ free( pHeap->pArr );
+ pHeap->pArr = pTmp;
+}
+
+
+void heap_term( heap_t * pHeap )
+{
+ assert( pHeap != NULL );
+
+ if ( pHeap->pArr != NULL )
+ free( pHeap->pArr );
+ memset( pHeap, 0, sizeof( heap_t ) );
+}
+
+
+void heap_insert( heap_t * pHeap, void * pData )
+{
+ int ind, father;
+ voidPtr_t pTmp;
+
+ resize( pHeap, pHeap->curr_size + 1 );
+
+ ind = pHeap->curr_size;
+ pHeap->curr_size++;
+
+ pHeap->pArr[ ind ] = pData;
+
+ while ( ind > 0 ) {
+ father = ( ind - 1 ) / 2;
+
+ if ( (*pHeap->pCompFunc)( pHeap->pArr[ ind ],
+ pHeap->pArr[ father ] ) <= 0 )
+ break;
+
+ pTmp = pHeap->pArr[ ind ];
+ pHeap->pArr[ ind ] = pHeap->pArr[ father ];
+ pHeap->pArr[ father ] = pTmp;
+
+ ind = father;
+ }
+}
+
+
+void * heap_top( heap_t * pHeap )
+{
+ return pHeap->pArr[ 0 ];
+}
+
+
+void * heap_delete_max( heap_t * pHeap )
+{
+ void * res;
+ void * pTmp;
+ int ind, left, right, max_ind;
+
+ if ( pHeap->curr_size <= 0 )
+ return NULL;
+
+ res = pHeap->pArr[ 0 ];
+
+ pHeap->curr_size--;
+ pHeap->pArr[ 0 ] = pHeap->pArr[ pHeap->curr_size ];
+ pHeap->pArr[ pHeap->curr_size ] = NULL;
+ ind = 0;
+
+ while ( ind < pHeap->curr_size ) {
+ left = 2 * ind + 1;
+ right = 2 * ind + 2;
+ if ( left >= pHeap->curr_size )
+ break;
+ if ( right >= pHeap->curr_size )
+ right = left;
+
+ if ( (*pHeap->pCompFunc)( pHeap->pArr[ left ],
+ pHeap->pArr[ right ] ) <= 0 )
+ max_ind = right;
+ else
+ max_ind = left;
+ if ( (*pHeap->pCompFunc)( pHeap->pArr[ ind ],
+ pHeap->pArr[ max_ind ] ) >= 0 )
+ break;
+
+ pTmp = pHeap->pArr[ ind ];
+ pHeap->pArr[ ind ] = pHeap->pArr[ max_ind ];
+ pHeap->pArr[ max_ind ] = pTmp;
+
+ ind = max_ind;
+ }
+
+ return res;
+}
+
+
+bool heap_is_empty( heap_t * pHeap )
+{
+ assert( pHeap != NULL );
+
+ return( pHeap->curr_size <= 0 );
+}
+
+
+int compare_pairs( void * _a, void * _b )
+{
+ const GFSPPair & a( *(GFSPPair *)_a );
+ const GFSPPair & b( *(GFSPPair *)_b );
+
+ if ( a.maxDiam() < b.maxDiam() )
+ return -1;
+ if ( a.maxDiam() > b.maxDiam() )
+ return 1;
+
+ return 0;
+}
+
+
+
+class g_heap_pairs_p
+{
+private:
+ heap_t heap;
+
+public:
+ g_heap_pairs_p() {
+ heap_init( &heap, compare_pairs );
+ }
+
+ ~g_heap_pairs_p() {
+ while ( size() > 0 ) {
+ pop();
+ }
+ heap_term( &heap );
+ }
+
+ void push( GFSPPair & pair )
+ {
+ GFSPPair * p_pair = new GFSPPair( pair );
+ heap_insert( &heap, (void *)p_pair );
+ }
+
+ int size() const {
+ return heap.curr_size;
+ }
+ GFSPPair top() {
+ return *(GFSPPair *)heap_top( &heap );
+ }
+ void pop() {
+ GFSPPair * ptr = (GFSPPair *)heap_delete_max( &heap );
+
+ delete ptr;
+ }
+};
+
+
+/* heap.implementation end */
+
+
+void computeExtremePair( const gdiam_point * arr, int size,
+ int dim, GPointPair & pp )
+{
+ pp.init( arr[ 0 ] );
+
+ for ( int ind = 1; ind < size; ind++ ) {
+ const gdiam_point pnt( arr[ ind ] );
+
+ if ( pnt[ dim ] < pp.p[ dim ] )
+ pp.p = pnt;
+ if ( pnt[ dim ] > pp.q[ dim ] )
+ pp.q = pnt;
+ }
+
+ pp.distance = pnt_distance( pp.p, pp.q );
+}
+
+
+
+void GTreeDiamAlg::compute_by_heap( double eps )
+{
+ g_heap_pairs_p heap;
+ int heap_limit, heap_delta;
+
+ heap_limit = HEAP_LIMIT;
+ heap_delta = HEAP_DELTA;
+
+ GFSPTreeNode * root = tree.getRoot();
+
+ computeExtremePair( tree.getPoints(), points_num,
+ root->getBBox().getLongestDim(),
+ pair_diam );
+
+ GFSPPair root_pair;
+
+ root_pair.init( root, root );
+
+ heap.push( root_pair );
+
+ int count =0;
+ while ( heap.size() > 0 ) {
+ GFSPPair pair = heap.top();
+ heap.pop();
+ split_pair( pair, heap, eps );
+ if ( ( count & 0x3ff ) == 0 ) {
+ if ( ((int)heap.size()) > heap_limit ) {
+ threshold_brute *= 2;
+ printf( "threshold_brute: %d\n", threshold_brute );
+ heap_limit += heap_delta;
+ }
+ }
+ count++;
+ }
+}
+
+
+void GTreeDiamAlg::compute_by_heap_proj( double eps,
+ gdiam_point proj )
+{
+ g_heap_pairs_p heap;
+ int heap_limit, heap_delta;
+ GPointPair pair_diam_x;
+
+ heap_limit = HEAP_LIMIT;
+ heap_delta = HEAP_DELTA;
+
+ GFSPTreeNode * root = tree.getRoot();
+
+ computeExtremePair( tree.getPoints(), points_num,
+ root->getBBox().getLongestDim(),
+ pair_diam_x );
+ pair_diam.init( pair_diam_x.p, pair_diam_x.q, proj );
+
+ GFSPPair root_pair;
+
+ root_pair.init( root, root, proj, 0 );
+
+ heap.push( root_pair );
+
+ int count =0;
+ while ( heap.size() > 0 ) {
+ GFSPPair pair = heap.top();
+ heap.pop();
+ split_pair_proj( pair, heap, eps, proj );
+ if ( ( count & 0x3ff ) == 0 ) {
+ if ( ((int)heap.size()) > heap_limit ) {
+ threshold_brute *= 2;
+ printf( "threshold_brute: %d\n", threshold_brute );
+ heap_limit += heap_delta;
+ }
+ }
+ count++;
+ }
+}
+
+
+void GTreeDiamAlg::addPairHeap( g_heap_pairs_p & heap,
+ GFSPTreeNode * left,
+ GFSPTreeNode * right )
+{
+ if ( update_diam_mode == UDM_SIMPLE ) {
+ const gdiam_point p( *(left->ptr_pnt_left()) );
+ const gdiam_point q( *(right->ptr_pnt_left()) );
+
+ pair_diam.update_diam( p, q );
+ } else
+ if ( update_diam_mode == UDM_BIG_SAMPLE ) {
+ if ( ( left->size() < 100 ) || ( right->size() < 100 ) ) {
+ const gdiam_point p( *(left->ptr_pnt_left()) );
+ const gdiam_point q( *(right->ptr_pnt_left()) );
+
+ pair_diam.update_diam( p, q );
+ } else
+ {
+#define UMD_SIZE 5
+ gdiam_point arr_left[ UMD_SIZE ],
+ arr_right[ UMD_SIZE ];
+ for ( int ind = 0; ind < UMD_SIZE; ind++ ) {
+ const gdiam_point p( *(left->ptr_pnt_rand()) );
+ const gdiam_point q( *(right->ptr_pnt_rand()) );
+ arr_left[ ind ] = p;
+ arr_right[ ind ] = q;
+ }
+ for ( int ind = 0; ind < UMD_SIZE; ind++ )
+ for ( int jnd = 0; jnd < UMD_SIZE; jnd++ )
+ pair_diam.update_diam( arr_left[ ind ],
+ arr_right[ jnd ] );
+ }
+ } else {
+ assert( false );
+ }
+
+ GFSPPair pair;
+
+ pair.init( left, right );
+ if ( pair.maxDiam() <= pair_diam.distance )
+ return;
+ heap.push( pair );
+}
+
+
+void GTreeDiamAlg::addPairHeap( g_heap_pairs_p & heap,
+ GFSPTreeNode * left,
+ GFSPTreeNode * right,
+ gdiam_point proj,
+ GFSPPair &UNUSED(father) )
+{
+ const gdiam_point p( *(left->ptr_pnt_left()) );
+ const gdiam_point q( *(right->ptr_pnt_left()) );
+
+ pair_diam.update_diam( p, q, proj );
+
+ GFSPPair pair;
+
+ pair.init( left, right, proj,
+ pnt_distance( p, q, proj ) );
+ if ( pair.maxDiam() <= pair_diam.distance )
+ return;
+ heap.push( pair );
+}
+
+
+void GTreeDiamAlg::split_pair_proj( GFSPPair & pair,
+ g_heap_pairs_p & heap, double eps,
+ gdiam_point proj )
+{
+ bool f_is_left_splitable, f_is_right_splitable;
+
+ if ( pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance ))
+ return;
+
+ if ( pair.isBelowThreshold( threshold_brute ) ) {
+ pair.directDiameter( tree, pair_diam, proj );
+ return;
+ }
+
+ f_is_left_splitable = pair.get_left()->isSplitable();
+ f_is_right_splitable = pair.get_right()->isSplitable();
+ assert( f_is_left_splitable || f_is_right_splitable );
+ if ( f_is_left_splitable )
+ tree.split_node( pair.get_left() );
+ if ( f_is_right_splitable )
+ tree.split_node( pair.get_right() );
+
+ if ( f_is_left_splitable
+ && f_is_right_splitable ) {
+ addPairHeap( heap,
+ pair.get_left()->get_left(),
+ pair.get_right()->get_left(),
+ proj, pair );
+ addPairHeap( heap,
+ pair.get_left()->get_left(),
+ pair.get_right()->get_right(),
+ proj, pair );
+ // to avoid exponential blowup
+ if ( pair.get_left() != pair.get_right() )
+ addPairHeap( heap,
+ pair.get_left()->get_right(),
+ pair.get_right()->get_left(),
+ proj, pair );
+ addPairHeap( heap,
+ pair.get_left()->get_right(),
+ pair.get_right()->get_right(),
+ proj, pair );
+ return;
+ }
+ if ( f_is_left_splitable ) {
+ addPairHeap( heap,
+ pair.get_left()->get_left(),
+ pair.get_right(), proj, pair );
+ addPairHeap( heap,
+ pair.get_left()->get_right(),
+ pair.get_right(), proj, pair );
+ return;
+ }
+ if ( f_is_right_splitable ) {
+ addPairHeap( heap,
+ pair.get_left(),
+ pair.get_right()->get_left(), proj, pair );
+ addPairHeap( heap,
+ pair.get_left(),
+ pair.get_right()->get_right(), proj, pair );
+ return;
+ }
+}
+
+
+void GTreeDiamAlg::split_pair( GFSPPair & pair,
+ g_heap_pairs_p & heap,
+ double eps )
+{
+ bool f_is_left_splitable, f_is_right_splitable;
+
+ if ( pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance ))
+ return;
+
+ if ( pair.isBelowThreshold( threshold_brute ) ) {
+ pair.directDiameter( tree, pair_diam );
+ return;
+ }
+
+ f_is_left_splitable = pair.get_left()->isSplitable();
+ f_is_right_splitable = pair.get_right()->isSplitable();
+ assert( f_is_left_splitable || f_is_right_splitable );
+ if ( f_is_left_splitable )
+ tree.split_node( pair.get_left() );
+ if ( f_is_right_splitable )
+ tree.split_node( pair.get_right() );
+
+ if ( f_is_left_splitable
+ && f_is_right_splitable ) {
+ addPairHeap( heap,
+ pair.get_left()->get_left(),
+ pair.get_right()->get_left() );
+ addPairHeap( heap,
+ pair.get_left()->get_left(),
+ pair.get_right()->get_right() );
+ // to avoid exponential blowup
+ if ( pair.get_left() != pair.get_right() )
+ addPairHeap( heap,
+ pair.get_left()->get_right(),
+ pair.get_right()->get_left() );
+ addPairHeap( heap,
+ pair.get_left()->get_right(),
+ pair.get_right()->get_right() );
+ return;
+ }
+ if ( f_is_left_splitable ) {
+ addPairHeap( heap,
+ pair.get_left()->get_left(),
+ pair.get_right() );
+ addPairHeap( heap,
+ pair.get_left()->get_right(),
+ pair.get_right() );
+ return;
+ }
+ if ( f_is_right_splitable ) {
+ addPairHeap( heap,
+ pair.get_left(),
+ pair.get_right()->get_left() );
+ addPairHeap( heap,
+ pair.get_left(),
+ pair.get_right()->get_right() );
+ return;
+ }
+}
+
+
+GPointPair gdiam_approx_diam( gdiam_point * start, int size,
+ gdiam_real eps )
+{
+ GPointPair pair;
+
+ GTreeDiamAlg * pAlg = new GTreeDiamAlg( start, size, UDM_SIMPLE );
+ pAlg->compute_by_heap( eps );
+
+ pair = pAlg->getDiameter();
+
+ delete pAlg;
+
+ return pair;
+}
+
+
+GPointPair gdiam_approx_diam_UDM( gdiam_point * start, int size,
+ gdiam_real eps )
+{
+ GPointPair pair;
+
+ GTreeDiamAlg * pAlg = new GTreeDiamAlg( start, size, UDM_BIG_SAMPLE );
+ pAlg->compute_by_heap( eps );
+
+ pair = pAlg->getDiameter();
+
+ delete pAlg;
+
+ return pair;
+}
+
+
+gdiam_point * gdiam_convert( gdiam_real * start, int size )
+{
+ gdiam_point * p_arr;
+
+ assert( start != NULL );
+ assert( size > 0 );
+
+ p_arr = (gdiam_point *)malloc( sizeof( gdiam_point ) * size );
+ assert( p_arr != NULL );
+
+ for ( int ind = 0; ind < size; ind++ )
+ p_arr[ ind ] = (gdiam_point)&(start[ ind * 3 ]);
+
+ return p_arr;
+}
+
+
+GPointPair gdiam_approx_diam_pair( gdiam_real * start, int size,
+ gdiam_real eps )
+{
+ gdiam_point * p_arr;
+ GPointPair pair;
+
+ p_arr = gdiam_convert( start, size );
+ pair = gdiam_approx_diam( p_arr, size, eps );
+ free( p_arr );
+
+ return pair;
+}
+
+
+GPointPair gdiam_approx_diam_pair_UDM( gdiam_real * start, int size,
+ gdiam_real eps )
+{
+ gdiam_point * p_arr;
+ GPointPair pair;
+
+ p_arr = gdiam_convert( start, size );
+ pair = gdiam_approx_diam_UDM( p_arr, size, eps );
+ free( p_arr );
+
+ return pair;
+}
+
+
+gdiam_bbox gdiam_approx_const_mvbb( gdiam_point * start, int size,
+ gdiam_real eps,
+ GBBox * p_ap_bbox )
+{
+ GPointPair pair, pair_2;
+
+ GTreeDiamAlg * pAlg = new GTreeDiamAlg( start, size, UDM_SIMPLE );
+ pAlg->compute_by_heap( eps );
+
+ pair = pAlg->getDiameter();
+
+ /* Compute the resulting diameter */
+ gdiam_point_t dir, dir_2, dir_3;
+
+ dir[ 0 ] = pair.p[ 0 ] - pair.q[ 0 ];
+ dir[ 1 ] = pair.p[ 1 ] - pair.q[ 1 ];
+ dir[ 2 ] = pair.p[ 2 ] - pair.q[ 2 ];
+ pnt_normalize( dir );
+
+ pAlg->compute_by_heap_proj( eps, dir );
+
+ pair_2 = pAlg->getDiameter();
+ dir_2[ 0 ] = pair_2.p[ 0 ] - pair_2.q[ 0 ];
+ dir_2[ 1 ] = pair_2.p[ 1 ] - pair_2.q[ 1 ];
+ dir_2[ 2 ] = pair_2.p[ 2 ] - pair_2.q[ 2 ];
+
+ gdiam_real prd;
+
+ prd = pnt_dot_prod( dir, dir_2 );
+ dir_2[ 0 ] = dir_2[ 0 ] - prd * dir[ 0 ];
+ dir_2[ 1 ] = dir_2[ 1 ] - prd * dir[ 1 ];
+ dir_2[ 2 ] = dir_2[ 2 ] - prd * dir[ 2 ];
+
+ pnt_normalize( dir );
+ pnt_normalize( dir_2 );
+
+ pnt_cross_prod( dir, dir_2, dir_3 );
+ pnt_normalize( dir_3 );
+
+ gdiam_bbox bbox;
+ GBBox ap_bbox;
+
+ if ( ( pnt_length( dir_2 ) < 1e-6 )
+ && ( pnt_length( dir_3 ) < 1e-6 ) )
+ gdiam_generate_orthonormal_base( dir, dir_2, dir_3 );
+
+ if ( ( pnt_length( dir ) == 0.0 )
+ || ( pnt_length( dir_2 ) < 1e-6 )
+ || ( pnt_length( dir_3 ) < 1e-6 ) ) {
+ gdiam_generate_orthonormal_base( dir, dir_2, dir_3 );
+ pnt_dump( dir );
+ pnt_dump( dir_2 );
+ pnt_dump( dir_3 );
+
+ fflush( stdout );
+ fflush( stderr );
+ assert( false );
+ }
+
+ bbox.init( dir, dir_2, dir_3 );
+ ap_bbox.init();
+ for ( int ind = 0; ind < size; ind++ ) {
+ bbox.bound( start[ ind ] );
+ ap_bbox.bound( start[ ind ] );
+ }
+
+ delete pAlg;
+ if ( ap_bbox.volume() < bbox.volume() )
+ bbox.init( ap_bbox );
+
+ if ( p_ap_bbox != NULL )
+ *p_ap_bbox = ap_bbox;
+
+ return bbox;
+}
+
+
+gdiam_bbox gdiam_approx_const_mvbb( gdiam_real * start, int size,
+ gdiam_real eps )
+{
+ gdiam_point * p_arr;
+ gdiam_bbox gbbox;
+
+ p_arr = gdiam_convert( start, size );
+ gbbox = gdiam_approx_const_mvbb( p_arr, size, eps, NULL );
+ free( p_arr );
+
+ return gbbox;
+}
+
+/*-----------------------------------------------------------------------
+ * Code for computing best bounding box in a given drection
+\*-----------------------------------------------------------------------*/
+
+void gdiam_generate_orthonormal_base( gdiam_point in,
+ gdiam_point out1,
+ gdiam_point out2 )
+{
+ assert( pnt_length( in ) > 0.0 );
+
+ pnt_normalize( in );
+
+ // stupid cases..
+ if ( in[ 0 ] == 0.0 ) {
+ if ( in[ 1 ] == 0.0 ) {
+ pnt_init_normalize( out1, 1, 0, 0 );
+ pnt_init_normalize( out2, 0, 1, 0 );
+ return;
+ }
+ if ( in[ 2 ] == 0.0 ) {
+ pnt_init_normalize( out1, 1, 0, 0 );
+ pnt_init_normalize( out2, 0, 0, 1 );
+ return;
+ }
+ pnt_init_normalize( out1, 0, -in[ 2 ], in[ 1 ] );
+ pnt_init_normalize( out2, 1, 0, 0 );
+ return;
+ }
+ if ( ( in[ 1 ] == 0.0 ) && ( in[ 2 ] == 0.0 ) ) {
+ pnt_init_normalize( out1, 0, 1, 0 );
+ pnt_init_normalize( out2, 0, 0, 1 );
+ return;
+ }
+ if ( in[ 1 ] == 0.0 ) {
+ pnt_init_normalize( out1, -in[ 2 ], 0, in[ 0 ] );
+ pnt_init_normalize( out2, 0, 1, 0 );
+ return;
+ }
+ if ( in[ 2 ] == 0.0 ) {
+ pnt_init_normalize( out1, -in[ 1 ], in[ 0 ], 0 );
+ pnt_init_normalize( out2, 0, 0, 1 );
+ return;
+ }
+
+ // all entries in the vector are not zero.
+ pnt_init_normalize( out1, -in[ 1 ], in[ 0 ], 0 );
+ pnt_cross_prod( in, out1, out2 );
+ pnt_normalize( out2 );
+}
+
+
+class point2d
+{
+public:
+ gdiam_real x, y;
+ gdiam_point src;
+
+ void init( gdiam_point pnt,
+ gdiam_point base_x,
+ gdiam_point base_y ) {
+ src = pnt;
+ x = pnt_dot_prod( pnt, base_x );
+ y = pnt_dot_prod( pnt, base_y );
+ }
+
+ gdiam_real dist( const point2d & pnt ) {
+ return sqrt( ( x - pnt.x ) * ( x - pnt.x )
+ + ( y - pnt.y ) * ( y - pnt.y ) );
+ }
+
+ void dump() const {
+ printf( "(%g, %g)\n", x, y );
+ }
+ bool equal( const point2d & pnt ) const {
+ return ( ( x == pnt.x ) && ( y == pnt.y ) );
+ }
+ bool equal_real( const point2d & pnt ) const {
+ return ( ( fabs( x - pnt.x ) < 1e-8 )
+ && ( fabs( y - pnt.y ) < 1e-8 ) );
+ }
+};
+
+typedef point2d * point2d_ptr;
+class vec_point_2d : public std::vector<point2d_ptr> {};
+
+
+
+inline ldouble Area( const point2d & a,
+ const point2d & b,
+ const point2d & c )
+{
+ double x1, y1, x2, y2, len;
+
+ x1 = b.x - a.x;
+ y1 = b.y - a.y;
+
+ x2 = c.x - a.x;
+ y2 = c.y - a.y;
+
+ if ( ( x1 != 0.0 ) || ( y1 != 0.0 ) ) {
+ len = sqrt( x1 * x1 + y1 * y1 );
+ x1 /= len;
+ y1 /= len;
+ }
+
+ if ( ( x2 != 0.0 ) || ( y2 != 0.0 ) ) {
+ len = sqrt( x2 * x2 + y2 * y2 );
+ x2 /= len;
+ y2 /= len;
+ }
+
+ ldouble area;
+
+ area = x1 * y2 - x2 * y1;
+ return area;
+}
+
+
+inline int AreaSign( const point2d & a,
+ const point2d & b,
+ const point2d & c )
+{
+ ldouble area;//, area1, area2;
+
+ area = a.x * b.y - a.y * b.x +
+ a.y * c.x - a.x * c.y +
+ b.x * c.y - c.x * b.y;
+
+ return ( ( area < (ldouble)0.0 )? -1:
+ ( (area > (ldouble)0.0)? 1 : 0 ) );
+}
+
+inline bool Left( const point2d & a,
+ const point2d & b,
+ const point2d & c )
+{
+ return AreaSign( a, b, c ) > 0;
+}
+
+inline bool Left_colinear( const point2d & a,
+ const point2d & b,
+ const point2d & c )
+{
+ return AreaSign( a, b, c ) >= 0;
+}
+
+
+struct CompareByAngle {
+public:
+ point2d base;
+
+ bool operator()(const point2d_ptr & a, const point2d_ptr & b ) {
+ int sgn;
+ gdiam_real len1, len2;
+
+ if ( a->equal( *b ) )
+ return false;
+ assert( a != NULL );
+ assert( b != NULL );
+ if ( a->equal( base ) ) {
+ assert( false );
+ return true;
+ }
+ if ( b->equal( base ) ) {
+ assert( false );
+ return false;
+ }
+
+ sgn = AreaSign( base, *a, *b );
+ if ( sgn != 0 )
+ return (sgn > 0);
+ len1 = base.dist( *a );
+ len2 = base.dist( *b );
+
+ return (len1 > len2);
+ }
+};
+
+
+point2d_ptr get_min_point( vec_point_2d & in,
+ int & index )
+{
+ point2d_ptr min_pnt = in[ 0 ];
+
+ index = 0;
+
+ for ( int ind = 0; ind < (int)in.size(); ind++ ) {
+ if ( in[ ind ]->y < min_pnt->y ) {
+ min_pnt = in[ ind ];
+ index = ind;
+ } else
+ if ( ( in[ ind ]->y == min_pnt->y )
+ && ( in[ ind ]->x < min_pnt->x ) ) {
+ min_pnt = in[ ind ];
+ index = ind;
+ }
+ }
+
+ return min_pnt;
+}
+
+
+void dump( vec_point_2d & vec )
+{
+ for ( int ind = 0; ind < (int)vec.size(); ind++ ) {
+ printf( "-- %11d (%-11g, %-11g)\n",
+ ind, vec[ ind ]->x,
+ vec[ ind ]->y );
+ }
+ printf( "\n\n" );
+}
+
+
+static void print_pnt( point2d_ptr pnt )
+{
+ printf( "(%g, %g)\n", pnt->x, pnt->y );
+}
+
+
+void verify_convex_hull( vec_point_2d & in, vec_point_2d & ch )
+{
+ ldouble area;
+
+ for ( int ind = 0; ind < (int)( ch.size() - 2 ); ind++ )
+ assert( Left( *( ch[ ind ] ), *( ch[ ind + 1 ] ),
+ *( ch[ ind + 2 ] ) ) );
+ assert( Left( *( ch[ ch.size() - 2 ] ), *( ch[ ch.size() - 1 ] ),
+ *( ch[ 0 ] ) ) );
+ assert( Left( *( ch[ ch.size() - 1 ] ), *( ch[ 0 ] ),
+ *( ch[ 1 ] ) ) );
+
+ for ( int ind = 0; ind < (int)( in.size() - 2 ); ind++ ) {
+ for ( int jnd = 0; jnd < (int)( ch.size() - 1 ); jnd++ ) {
+ if ( ch[ jnd ]->equal_real( *( in[ ind ] ) ) )
+ continue;
+ if ( ch[ jnd + 1 ]->equal( *( in[ ind ] ) ) )
+ continue;
+
+ if ( ! Left_colinear( *( ch[ jnd ] ), *( ch[ jnd + 1 ] ),
+ *( in[ ind ] ) ) ) {
+ area = Area( *( ch[ jnd ] ), *( ch[ jnd + 1 ] ),
+ *( in[ ind ] ) );
+ if ( fabs( area ) < 1e-12 )
+ continue;
+ printf( "Failure in progress!\n\n" );
+ print_pnt( ch[ jnd ] );
+ print_pnt( ch[ jnd + 1 ] );
+ print_pnt( in[ ind ] );
+
+ printf( "ch[ jnd ]: (%g, %g)\n",
+ ch[ jnd ]->x, ch[ jnd ]->y );
+ printf( "ch[ jnd + 1 ]: (%g, %g)\n",
+ ch[ jnd + 1 ]->x, ch[ jnd + 1 ]->y );
+ printf( "ch[ ind ]: (%g, %g)\n",
+ in[ ind ]->x, in[ ind ]->y );
+ printf( "Area: %g\n", (double)Area( *( ch[ jnd ] ),
+ *( ch[ jnd + 1 ] ),
+ *( in[ ind ] ) ) );
+ printf( "jnd: %d, jnd + 1: %d, ind: %d\n",
+ jnd, jnd+1, ind );
+ fflush( stdout );
+ fflush( stderr );
+
+ assert( Left( *( ch[ jnd ] ), *( ch[ jnd + 1 ] ),
+ *( in[ ind ] ) ) );
+ }
+ }
+ if ( ch[ 0 ]->equal( *( in[ ind ] ) ) )
+ continue;
+ if ( ch[ ch.size() - 1 ]->equal( *( in[ ind ] ) ) )
+ continue;
+ assert( Left( *( ch[ ch.size() - 1 ] ),
+ *( ch[ 0 ] ),
+ *( in[ ind ] ) ) );
+ }
+ printf( "Convex hull seems to be ok!\n" );
+}
+
+
+
+static void remove_duplicate( vec_point_2d & in,
+ point2d_ptr ptr,
+ int start )
+{
+ int dest;
+
+ dest = start;
+ while ( start < (int)in.size() ) {
+ if ( in[ start ]->equal_real( *ptr ) ) {
+ start++;
+ continue;
+ }
+ in[ dest++ ] = in[ start++ ];
+ }
+ while ( (int)in.size() > dest )
+ in.pop_back();
+}
+
+
+static void remove_consecutive_dup( vec_point_2d & in )
+{
+ int dest, pos;
+
+ if ( in.size() < 1 )
+ return;
+
+ dest = pos = 0;
+
+ // copy first entry
+ in[ dest++ ] = in[ pos++ ];
+ while ( pos < ( (int)in.size() - 1 ) ) {
+ if ( in[ pos - 1 ]->equal_real( *(in[ pos ]) ) ) {
+ pos++;
+ continue;
+ }
+ in[ dest++ ] = in[ pos++ ];
+ }
+ in[ dest++ ] = in[ pos++ ];
+
+ while ( (int)in.size() > dest )
+ in.pop_back();
+}
+
+
+// Copyright 2001 softSurfer, 2012 Dan Sunday
+// This code may be freely used and modified for any purpose
+// providing that this copyright notice is included with it.
+// SoftSurfer makes no warranty for this code, and cannot be held
+// liable for any real or imagined damage resulting from its use.
+// Users of this code must verify correctness for their application.
+// isLeft(): tests if a point is Left|On|Right of an infinite line.
+// Input: three points P0, P1, and P2
+// Return: >0 for P2 left of the line through P0 and P1
+// =0 for P2 on the line
+// <0 for P2 right of the line
+// See: Algorithm 1 on Area of Triangles
+inline float
+isLeft( const point2d &P0, const point2d &P1, const point2d &P2 )
+{
+ return (P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y);
+}
+
+
+
+
+// Switched to Monotone Chain Algorithm to avoid the area computation, which
+// is problematic for strict weak ordering when doing sorting. Useful links:
+//
http://stackoverflow.com/questions/1041620/most-efficient-way-to-erase-duplicates-and-sort-a-c-vector
+// http://geomalgorithms.com/a10-_hull-1.html
+void convex_hull( vec_point_2d & in, vec_point_2d & out )
+{
+
+ CompareByAngle comp;
+ int position;
+
+ assert( in.size() > 1 );
+
+ comp.base = *( get_min_point( in, position ));
+ std::swap( in[ 0 ], in[ position ] );
+
+ remove_duplicate( in, in[ 0 ], 1 );
+
+ for ( int ind = 0; ind < (int)in.size(); ind++ ) {
+ assert( in[ ind ] != NULL );
+ }
+
+ // Copyright 2001 softSurfer, 2012 Dan Sunday
+ // This code may be freely used and modified for any purpose
+ // providing that this copyright notice is included with it.
+ // SoftSurfer makes no warranty for this code, and cannot be held
+ // liable for any real or imagined damage resulting from its use.
+ // Users of this code must verify correctness for their application.
+
+ int n = in.size();
+ sort( in.begin() + 1, in.end(), comp );
+ remove_consecutive_dup( in );
+ out.reserve(n);
+
+ // the output array out[] will be used as the stack
+ int bot=0, top=(-1); // indices for bottom and top of the stack
+ int i; // array scan index
+
+
+ // Get the indices of points with min x-coord and min|max y-coord
+ int minmin = 0, minmax;
+ float xmin = in[0]->x;
+ for (i=1; i<n; i++)
+ if (in[i]->x != xmin) break;
+ minmax = i-1;
+ if (minmax == n-1) { // degenerate case: all x-coords == xmin
+ out.push_back(in[minmin]);
+ if (in[minmax]->y != in[minmin]->y) // a nontrivial segment
+ out.push_back(in[minmax]);
+ out.push_back(in[minmin]);
+ return;
+ }
+
+ // Get the indices of points with max x-coord and min|max y-coord
+ int maxmin, maxmax = n-1;
+ float xmax = in[n-1]->x;
+ for (i=n-2; i>=0; i--)
+ if (in[i]->x != xmax) break;
+ maxmin = i+1;
+
+ // Compute the lower hull on the stack H
+ out.push_back(in[minmin]); ++top; // push minmin point onto stack
+ i = minmax;
+ while (++i <= maxmin)
+ {
+ // the lower line joins in[minmin] with in[maxmin]
+ if (isLeft( *in[minmin], *in[maxmin], *in[i]) >= 0 && i < maxmin)
+ continue; // ignore in[i] above or on the lower line
+
+ while (top > 0) // there are at least 2 points on the stack
+ {
+ // test if in[i] is left of the line at the stack top
+ if (isLeft( *out[top-1], *out[top], *in[i]) > 0)
+ break; // in[i] is a new hull vertex
+ else {
+ --top;
+ out.pop_back(); // pop top point off stack
+ }
+ }
+ ++top;
+ out.push_back(in[i]); // push in[i] onto stack
+ }
+ // Next, compute the upper hull on the stack out above the bottom hull
+ if (maxmax != maxmin) // if distinct xmax points
+ out.push_back(in[maxmax]); // push maxmax point onto stack
+ bot = top; // the bottom point of the upper hull stack
+ i = maxmin;
+ while (--i > minmax)
+ {
+ // the upper line joins in[maxmax] with in[minmax]
+ if (isLeft( *in[maxmax], *in[minmax], *in[i]) >= 0 && i > minmax)
+ continue; // ignore in[i] below or on the upper line
+
+
+ while (top > bot) // at least 2 points on the upper stack
+ {
+ // test if in[i] is left of the line at the stack top
+ if (isLeft( *out[top-1], *out[top], *in[i]) > 0) {
+ break; // in[i] is a new hull vertex
+ } else {
+ --top;
+ out.pop_back(); // pop top point off stack
+ }
+ }
+ ++top;
+ out.push_back(in[i]); // push in[i] onto stack
+ }
+ if (minmax != minmin)
+ out.push_back(in[minmin]); // push joining endpoint onto stack
+#ifndef GDIAM_QUIET
+ std::vector<point2d_ptr>::iterator it;
+ for(it = out.begin(); it != out.end(); it++){
+ std::cout << "point: " << (*it)->x << "," << (*it)->y << "\n";
+ }
+#endif
+
+}
+
+
+struct bbox_2d_info
+{
+ gdiam_point_2d_t vertices[ 4 ];
+ gdiam_real area;
+
+ void get_dir( int ind, gdiam_point_2d out ) {
+ out[ 0 ] = vertices[ ind ][ 0 ]
+ - vertices[ (ind + 1) % 4 ][ 0 ];
+ out[ 1 ] = vertices[ ind ][ 1 ]
+ - vertices[ (ind + 1) % 4 ][ 1 ];
+ }
+
+ void dump() {
+ printf( "--- bbox 2d ----------------------------\n" );
+ for ( int ind = 0; ind < 4; ind++ )
+ printf( "%d: (%g, %g)\n", ind, vertices[ ind ][0],
+ vertices[ ind ][1] );
+ for ( int ind = 0; ind < 4; ind++ ) {
+ gdiam_point_2d_t dir;
+
+ get_dir( ind, dir );
+
+ printf( "dir %d: (%g, %g)\n", ind, dir[0],
+ dir[1] );
+ }
+ printf( "\\----------------------------------\n" );
+ }
+
+};
+
+#define EPS 1e-6
+inline bool eq_real( gdiam_real a, gdiam_real b ) {
+ return ( ( ( b - EPS ) < a )
+ && ( a < (b + EPS) ) );
+}
+
+class MinAreaRectangle
+{
+private:
+ vec_point_2d ch;
+ gdiam_real * angles;
+
+public:
+ void dump_ch() {
+ for ( int ind = 0; ind < (int)ch.size(); ind++ ) {
+ printf( "%d: (%g, %g)\n", ind, ch[ ind ]->x,
+ ch[ ind ]->y );
+ }
+ }
+
+ void compute_edge_dir( int u ) {
+ int u2;
+ double ang;
+ double x1, y1, x2, y2;
+
+ u2 = (u + 1) % ch.size();
+
+ x1 = ch[ u ]->x;
+ y1 = ch[ u ]->y;
+ x2 = ch[ u2 ]->x;
+ y2 = ch[ u2 ]->y;
+
+ ang = atan2( y2 - y1, x2 - x1 );
+ if ( ang < 0 )
+ ang += 2.0 * PI;
+ angles[ u ] = ang;
+
+ // floating point nonsence. A huck to solve this...
+ if ( ( u > 0 )
+ && ( angles[ u ] < angles[ u - 1 ] )
+ && eq_real( angles[ u ], angles[ u - 1 ] ) )
+ angles[ u ] = angles[ u - 1 ];
+ }
+
+ /* Finding the first vertex whose edge is over a given angle */
+ int find_vertex( int start_vertex,
+ double trg_angle ) {
+ double prev_angle, curr_angle;
+ bool f_found;
+ int ver, prev_vertex;
+
+ f_found = false;
+
+ prev_vertex = start_vertex - 1;
+ if ( prev_vertex < 0 )
+ prev_vertex = ch.size() - 1;
+
+ prev_angle = angles[ prev_vertex ];
+ ver = start_vertex;
+ while ( ! f_found ) {
+ curr_angle = angles[ ver ];
+ if ( prev_angle <= curr_angle )
+ f_found = ( ( prev_angle < trg_angle)
+ && ( trg_angle <= curr_angle ) );
+ else
+ f_found = ( ( prev_angle < trg_angle )
+ || ( trg_angle <= curr_angle ));
+
+ if ( f_found)
+ break;
+ else
+ prev_angle = curr_angle;
+
+ ver = ( ver + 1 ) % ch.size();
+ }
+
+ return ver;
+ }
+
+
+ /* Computing the intersection point of two lines, each given by a
+ point and slope */
+ void compute_crossing( int a_ind, double a_angle,
+ int b_ind, double b_angle,
+ gdiam_real & x,
+ gdiam_real & y ) {
+ double a_slope, b_slope, x_a, x_b, y_a, y_b;
+
+ if ( eq_real( a_angle, 0.0 ) || eq_real( a_angle, PI ) ) {
+ x = ch[ b_ind ]->x;
+ y = ch[ a_ind ]->y;
+ return;
+ }
+ if ( eq_real( a_angle, PI/2.0 ) || eq_real( a_angle, PI * 1.5 ) ) {
+ x = ch[ a_ind ]->x;
+ y = ch[ b_ind ]->y;
+ return;
+ }
+
+ a_slope = tan( a_angle );
+ b_slope = tan( b_angle );
+ x_a = ch[ a_ind ]->x;
+ y_a = ch[ a_ind ]->y;
+ x_b = ch[ b_ind ]->x;
+ y_b = ch[ b_ind ]->y;
+
+ double a_const, b_const, x_out, y_out;
+ a_const = y_a - a_slope * x_a;
+ b_const = y_b - b_slope * x_b;
+ x_out = - ( a_const - b_const ) / (a_slope - b_slope);
+ y_out = a_slope * x_out + a_const;
+ x = x_out;
+ y = y_out;
+ }
+
+
+ void get_bbox( int a_ind, gdiam_real a_angle,
+ int b_ind, gdiam_real b_angle,
+ int c_ind, gdiam_real c_angle,
+ int d_ind, gdiam_real d_angle,
+ bbox_2d_info & bbox,
+ gdiam_real & area ) {
+ compute_crossing( a_ind, a_angle, b_ind, b_angle,
+ bbox.vertices[ 0 ][ 0 ],
+ bbox.vertices[ 0 ][ 1 ] );
+ compute_crossing( b_ind, b_angle, c_ind, c_angle,
+ bbox.vertices[ 1 ][ 0 ],
+ bbox.vertices[ 1 ][ 1 ] );
+ compute_crossing( c_ind, c_angle, d_ind, d_angle,
+ bbox.vertices[ 2 ][ 0 ],
+ bbox.vertices[ 2 ][ 1 ] );
+ compute_crossing( d_ind, d_angle, a_ind, a_angle,
+ bbox.vertices[ 3 ][ 0 ],
+ bbox.vertices[ 3 ][ 1 ] );
+
+ area = pnt_distance_2d( bbox.vertices[ 0 ],
+ bbox.vertices[ 1 ] )
+ * pnt_distance_2d( bbox.vertices[ 0 ],
+ bbox.vertices[ 3 ] );
+ }
+
+ /* Given one angle (bbx direction), find the other three */
+ void get_angles( int ind,
+ gdiam_real & angle_1,
+ gdiam_real & angle_2,
+ gdiam_real & angle_3 ) {
+ double angle;
+
+ angle = angles[ ind ];
+ angle_1 = angle + PI * 0.5;
+
+ if ( angle_1 >= (2.0 * PI) )
+ angle_1 -= 2.0 * PI;
+
+ angle_2 = angle + PI;
+ if ( angle_2 >= (2.0 * PI) )
+ angle_2 -= 2.0 * PI;
+
+ angle_3 = angle + 1.5 * PI;
+ if ( angle_3 >= (2.0 * PI) )
+ angle_3 -= 2.0 * PI;
+ }
+
+ void compute_min_bbox_inner( bbox_2d_info & min_bbox,
+ gdiam_real & min_area ) {
+ int u, v, s, t;
+ gdiam_real ang1, ang2, ang3, tmp_area;
+ bbox_2d_info tmp_bbox;
+
+ angles = (gdiam_real *)malloc( sizeof( gdiam_real )
+ * (int)ch.size() );
+ assert( angles != NULL );
+
+ /* Pre-computing all edge directions */
+ for (u = 0; u < (int)ch.size(); u ++)
+ compute_edge_dir( u );
+
+ /* Initializing the search */
+ get_angles( 0, ang1, ang2, ang3 );
+
+ s = find_vertex( 0, ang1 );
+ v = find_vertex( s, ang2 );
+ t = find_vertex( v, ang3 );
+
+ get_bbox( 0, angles[ 0 ],
+ s, ang1,
+ v, ang2,
+ t, ang3,
+ min_bbox, min_area );
+
+ /* Performing a double rotating calipers */
+ for ( u = 1; u < (int)ch.size(); u ++ ) {
+ get_angles( u, ang1, ang2, ang3 );
+ s = find_vertex( s, ang1 );
+ v = find_vertex( v, ang2 );
+ t = find_vertex( t, ang3 );
+ get_bbox( u, angles[ u ],
+ s, ang1,
+ v, ang2,
+ t, ang3,
+ tmp_bbox, tmp_area );
+ if ( min_area > tmp_area ) {
+ min_area = tmp_area;
+ min_bbox = tmp_bbox;
+ }
+ }
+ free( angles );
+ angles = NULL;
+ }
+
+ void compute_min_bbox( vec_point_2d & points,
+ bbox_2d_info & min_bbox,
+ gdiam_real & min_area ) {
+ convex_hull( points, ch );
+ compute_min_bbox_inner( min_bbox, min_area );
+ }
+
+};
+
+void dump_points( gdiam_point * in_arr,
+ int size )
+{
+ for ( int ind = 0; ind < size; ind++ ) {
+ printf( "%d: (%g, %g, %g)\n", ind,
+ in_arr[ ind ][ 0 ],
+ in_arr[ ind ][ 1 ],
+ in_arr[ ind ][ 2 ] );
+ }
+}
+
+#define ZERO_EPS (1e-6)
+
+static void construct_base_inner( gdiam_point pnt1,
+ gdiam_point pnt2,
+ gdiam_point pnt3 )
+{
+ pnt_normalize( pnt1 );
+ pnt_normalize( pnt2 );
+ pnt_normalize( pnt3 );
+
+ if ( ( pnt_length( pnt1 ) < ZERO_EPS )
+ && ( pnt_length( pnt2 ) < ZERO_EPS )
+ && ( pnt_length( pnt3 ) < ZERO_EPS ) ){
+ assert( false );
+ }
+
+ if ( ( pnt_length( pnt1 ) < ZERO_EPS )
+ && ( pnt_length( pnt2 ) < ZERO_EPS ) ) {
+ gdiam_generate_orthonormal_base( pnt3, pnt1, pnt2 );
+ return;
+ }
+ if ( ( pnt_length( pnt1 ) < ZERO_EPS )
+ && ( pnt_length( pnt3 ) < ZERO_EPS ) ) {
+ gdiam_generate_orthonormal_base( pnt2, pnt1, pnt3 );
+ return;
+ }
+ if ( ( pnt_length( pnt2 ) < ZERO_EPS )
+ && ( pnt_length( pnt3 ) < ZERO_EPS ) ) {
+ gdiam_generate_orthonormal_base( pnt1, pnt2, pnt3 );
+ return;
+ }
+ if ( pnt_length( pnt1 ) < ZERO_EPS ) {
+ pnt_cross_prod( pnt2, pnt3, pnt1 );
+ return;
+ }
+ if ( pnt_length( pnt2 ) < ZERO_EPS ) {
+ pnt_cross_prod( pnt1, pnt3, pnt2 );
+ return;
+ }
+ if ( pnt_length( pnt3 ) < ZERO_EPS ) {
+ pnt_cross_prod( pnt1, pnt2, pnt3 );
+ return;
+ }
+}
+
+void construct_base( gdiam_point pnt1,
+ gdiam_point pnt2,
+ gdiam_point pnt3 )
+{
+ construct_base_inner( pnt1, pnt2, pnt3 );
+ pnt_scale_and_add( pnt2, -pnt_dot_prod( pnt1, pnt2 ),
+ pnt1 );
+ pnt_normalize( pnt2 );
+
+ pnt_scale_and_add( pnt3, -pnt_dot_prod( pnt1, pnt3 ),
+ pnt1 );
+ pnt_scale_and_add( pnt3, -pnt_dot_prod( pnt2, pnt3 ),
+ pnt2 );
+ pnt_normalize( pnt3 );
+}
+
+
+class ProjPointSet
+{
+private:
+ gdiam_point_t base_x, base_y, base_proj;
+ point2d * arr;
+ vec_point_2d points, ch;
+ gdiam_bbox bbox;
+ gdiam_point * in_arr;
+ int size;
+
+public:
+ ProjPointSet() : arr(NULL), in_arr(NULL), size( 0 ) {}
+
+ ~ProjPointSet() {
+ term();
+ }
+ void init( gdiam_point dir,
+ gdiam_point * _in_arr,
+ int _size ) {
+ arr = NULL;
+
+ if ( pnt_length( dir ) == 0.0 ) {
+ dump_points( _in_arr, _size );
+ pnt_dump( dir );
+ fflush( stdout );
+ fflush( stderr );
+ assert( pnt_length( dir ) > 0.0 );
+ }
+ size = _size;
+ in_arr = _in_arr;
+ assert( size > 0 );
+
+ pnt_copy( base_proj, dir );
+ gdiam_generate_orthonormal_base( base_proj, base_x, base_y );
+
+ arr = (point2d *)malloc( sizeof( point2d ) * size );
+ assert( arr != 0 );
+
+ for ( int ind = 0; ind < size; ind++ ) {
+ arr[ ind ].init( in_arr[ ind ], base_x, base_y );
+ points.push_back( &( arr[ ind ] ) );
+ }
+ }
+
+ void compute() {
+ MinAreaRectangle mar;
+ bbox_2d_info min_bbox;
+ gdiam_real min_area;
+ double x1, y1, x2, y2;
+ gdiam_point_t out_1, out_2;
+
+ mar.compute_min_bbox( points, min_bbox, min_area );
+
+ // We take the resulting min_area rectangle and lift it back to 3d...
+ x1 = min_bbox.vertices[ 0 ][ 0 ] -
+ min_bbox.vertices[ 1 ][ 0 ];
+ y1 = min_bbox.vertices[ 0 ][ 1 ] -
+ min_bbox.vertices[ 1 ][ 1 ];
+ x2 = min_bbox.vertices[ 0 ][ 0 ] -
+ min_bbox.vertices[ 3 ][ 0 ];
+ y2 = min_bbox.vertices[ 0 ][ 1 ] -
+ min_bbox.vertices[ 3 ][ 1 ];
+
+ double len;
+
+ len = sqrt( x1 * x1 + y1 * y1 );
+ if ( len > 0.0 ) {
+ x1 /= len;
+ y1 /= len;
+ }
+
+ len = sqrt( x2 * x2 + y2 * y2 );
+ if ( len > 0.0 ) {
+ x2 /= len;
+ y2 /= len;
+ }
+
+ pnt_zero( out_1 );
+ pnt_scale_and_add( out_1, x1, base_x );
+ pnt_scale_and_add( out_1, y1, base_y );
+ pnt_normalize( out_1 );
+
+ pnt_zero( out_2 );
+ pnt_scale_and_add( out_2, x2, base_x );
+ pnt_scale_and_add( out_2, y2, base_y );
+ pnt_normalize( out_2 );
+
+ construct_base( base_proj, out_1, out_2 );
+ if ( ( ! (fabs( pnt_dot_prod( base_proj, out_1 ) ) < 1e-6 ) )
+ || ( ! (fabs( pnt_dot_prod( base_proj, out_2 ) ) < 1e-6 ) )
+ || ( ! (fabs( pnt_dot_prod( out_1, out_2 ) ) < 1e-6 ) ) ) {
+ printf( "should be all close to zero: %g, %g, %g\n",
+ pnt_dot_prod( base_proj, out_1 ),
+ pnt_dot_prod( base_proj, out_2 ),
+ pnt_dot_prod( out_1, out_2 ) );
+ pnt_dump( base_proj );
+ pnt_dump( out_1 );
+ pnt_dump( out_2 );
+
+
+ printf( "real points:\n" );
+ dump_points( in_arr, size );
+
+ printf( "points on CH:\n" );
+ dump( points );
+
+ printf( "BBox:\n" );
+ min_bbox.dump();
+
+ fflush( stdout );
+ fflush( stderr );
+ assert( fabs( pnt_dot_prod( base_proj, out_1 ) ) < 1e-6 );
+ assert( fabs( pnt_dot_prod( base_proj, out_2 ) ) < 1e-6 );
+ assert( fabs( pnt_dot_prod( out_1, out_2 ) ) < 1e-6 );
+ }
+
+ bbox.init( base_proj, out_1, out_2 );
+ for ( int ind = 0; ind < size; ind++ )
+ bbox.bound( in_arr[ ind ] );
+ }
+
+ gdiam_bbox & get_bbox() {
+ return bbox;
+ }
+
+ void term() {
+ if ( arr != NULL ) {
+ free( arr );
+ arr = NULL;
+ }
+ }
+ void init() {
+ }
+};
+
+
+gdiam_bbox gdiam_mvbb_optimize( gdiam_point * start, int size,
+ gdiam_bbox bb_out, int times )
+{
+ gdiam_bbox bb_tmp;
+
+ for ( int ind = 0; ind < times; ind++ ) {
+ ProjPointSet pps;
+
+ if ( pnt_length( bb_out.get_dir( ind % 3 ) ) == 0.0 ) {
+ printf( "Dumping!\n" );
+ bb_out.dump();
+ fflush( stdout );
+ continue;
+ }
+
+ pps.init( bb_out.get_dir( ind % 3 ), start, size );
+
+ pps.compute();
+ bb_tmp = pps.get_bbox();
+
+ if ( bb_tmp.volume() < bb_out.volume() )
+ bb_out = bb_tmp;
+ }
+
+ return bb_out;
+}
+
+
+gdiam_bbox gdiam_approx_mvbb( gdiam_point * start, int size,
+ gdiam_real UNUSED(eps) )
+{
+ gdiam_bbox bb, bb2;
+
+ bb = gdiam_approx_const_mvbb( start, size, 0.0, NULL );
+ bb2 = gdiam_mvbb_optimize( start, size, bb, 10 );
+
+ return bb2;
+}
+
+
+static int gcd2( int a, int b )
+{
+ if ( a == 0 || a == b )
+ return b;
+ if ( b == 0 )
+ return a;
+ if ( a > b )
+ return gcd2( a % b, b );
+ else
+ return gcd2( a, b % a );
+}
+
+
+static int gcd3( int a, int b, int c )
+{
+ a = abs( a );
+ b = abs( b );
+ c = abs( c );
+ if ( a == 0 )
+ return gcd2( b, c );
+ if ( b == 0 )
+ return gcd2( a, c );
+ if ( c == 0 )
+ return gcd2( a, b );
+
+ return gcd2( a, gcd2( b, c ) );
+}
+
+
+static void try_direction( gdiam_bbox & bb,
+ gdiam_bbox & bb_out,
+ gdiam_point * start,
+ int size,
+ int x_coef,
+ int y_coef,
+ int z_coef )
+{
+ gdiam_bbox bb_new;
+
+ if ( gcd3( x_coef, y_coef, z_coef ) > 1 )
+ return;
+ if ( ( x_coef == 0 ) && ( y_coef == 0 )
+ && ( z_coef == 0 ) )
+ return;
+ gdiam_point_t new_dir;
+ bb.combine( new_dir, x_coef, y_coef, z_coef );
+
+ ProjPointSet pps;
+
+ pps.init( new_dir, start, size );
+
+ pps.compute();
+
+ bb_new = pps.get_bbox();
+ bb_new = gdiam_mvbb_optimize( start, size, bb_new, 10 );
+
+ if ( bb_new.volume() < bb_out.volume() )
+ bb_out = bb_new;
+}
+
+
+gdiam_bbox gdiam_approx_mvbb_grid( gdiam_point * start, int size,
+ int grid_size )
+{
+ gdiam_bbox bb, bb_out;
+
+ bb_out = bb = gdiam_approx_const_mvbb( start, size, 0.0, NULL );
+ if ( bb.volume() == 0 ) {
+ dump_points( start, size );
+ printf( "1zero volume???\n" );
+ bb.dump();
+ }
+
+ bb_out = bb = gdiam_mvbb_optimize( start, size, bb_out, 10 );
+
+ if ( bb.volume() == 0 ) {
+ printf( "2zero volume???\n" );
+ bb.dump();
+ }
+
+ for ( int x_coef = -grid_size; x_coef <= grid_size; x_coef++ )
+ for ( int y_coef = -grid_size; y_coef <= grid_size; y_coef++ )
+ for ( int z_coef = 0; z_coef <= grid_size; z_coef++ )
+ try_direction( bb, bb_out, start, size,
+ x_coef, y_coef, z_coef );
+
+ bb_out = gdiam_mvbb_optimize( start, size, bb_out, 10 );
+
+ return bb_out;
+}
+
+
+static void register_point( gdiam_point pnt,
+ gdiam_point * tops,
+ gdiam_point * bottoms,
+ int grid_size,
+ gdiam_bbox & bb )
+{
+ gdiam_point_t pnt_bb, pnt_bottom, pnt_top;
+ int x_ind, y_ind, position;
+
+ bb.get_normalized_coordinates( pnt, pnt_bb );
+
+ // x_ind
+ x_ind = (int)( pnt_bb[ 0 ] * grid_size );
+ assert( ( -1 <= x_ind ) && ( x_ind <= grid_size ) );
+ if ( x_ind < 0 )
+ x_ind = 0;
+ if ( x_ind >= grid_size )
+ x_ind = grid_size - 1;
+
+ // y_ind
+ y_ind = (int)( pnt_bb[ 1 ] * grid_size );
+ assert( ( -1 <= y_ind ) && ( y_ind <= grid_size ) );
+ if ( y_ind < 0 )
+ y_ind = 0;
+ if ( y_ind >= grid_size )
+ y_ind = grid_size - 1;
+
+ position = x_ind + y_ind * grid_size;
+ if ( tops[ position ] == NULL ) {
+ tops[ position ] = bottoms[ position ] = pnt;
+ return;
+ }
+
+ bb.get_normalized_coordinates( tops[ position ], pnt_top );
+ if ( pnt_top[ 2 ] < pnt_bb[ 2 ] )
+ tops[ position ] = pnt;
+
+ bb.get_normalized_coordinates( bottoms[ position ], pnt_bottom );
+ if ( pnt_bottom[ 2 ] > pnt_bb[ 2 ] )
+ bottoms[ position ] = pnt;
+}
+
+
+
+// gdiam_convex_sample
+
+// We are given a point set, and (hopefully) a tight fitting
+// bounding box. We compute a sample of the given size that represents
+// the point-set. The only guarenteed is that if we use ample of size m,
+// we get an approximation of quality about 1/\sqrt{m}. Note that we pad
+// the sample if necessary to get the desired size.
+gdiam_point * gdiam_convex_sample( gdiam_point * start, int size,
+ gdiam_bbox & bb,
+ int sample_size )
+{
+ int grid_size, mem_size, grid_entries;
+ gdiam_point * bottoms, * tops, * out_arr;
+ int out_count;
+
+ assert( sample_size > 1 );
+
+ // we want the grid to be on the two minor dimensions.
+ bb.set_third_dim_longest();
+
+ grid_size = (int)sqrt((double)(sample_size / 2) );
+
+ grid_entries = grid_size * grid_size;
+ mem_size = (int)( sizeof( gdiam_point ) * grid_size * grid_size );
+
+ bottoms = (gdiam_point *)malloc( mem_size );
+ tops = (gdiam_point *)malloc( mem_size );
+ out_arr = (gdiam_point *)malloc( sizeof( gdiam_point ) * sample_size );
+
+ assert( bottoms != NULL );
+ assert( tops != NULL );
+ assert( out_arr != NULL );
+
+ for ( int ind = 0; ind < grid_entries; ind++ )
+ tops[ ind ] = bottoms[ ind ] = NULL;
+
+ // Now we stream the points registering them with the relevant
+ // shaft in the grid.
+ for ( int ind = 0; ind < size; ind++ )
+ register_point( start[ ind ], tops, bottoms,
+ grid_size, bb );
+
+ out_count = 0;
+ for ( int ind = 0; ind < grid_entries; ind++ ) {
+ if ( tops[ ind ] == NULL )
+ continue;
+
+ out_arr[ out_count++ ] = tops[ ind ];
+ if ( tops[ ind ] != bottoms[ ind ] )
+ out_arr[ out_count++ ] = bottoms[ ind ];
+ }
+
+ // We dont have enough entries in our sample - lets randomly pick
+ // points.
+ while ( out_count < sample_size )
+ out_arr[ out_count++ ] = start[ rand() % size ];
+
+ free( tops );
+ free( bottoms );
+
+ return out_arr;
+}
+
+
+gdiam_bbox gdiam_approx_mvbb_grid_sample( gdiam_point * start, int size,
+ int grid_size, int sample_size )
+{
+ gdiam_bbox bb, bb_new;
+ gdiam_point * sample;
+
+ if ( sample_size >= size )
+ return gdiam_approx_mvbb_grid( start, size, grid_size );
+
+ bb = gdiam_approx_const_mvbb( start, size, 0.0, NULL );
+
+ sample = gdiam_convex_sample( start, size, bb, sample_size );
+
+ bb_new = gdiam_approx_mvbb_grid( sample, sample_size, grid_size );
+
+ for ( int ind = 0; ind < size; ind++ )
+ bb_new.bound( start[ ind ] );
+
+ return bb_new;
+}
+
+
+
+gdiam_bbox gdiam_approx_mvbb_grid_sample( gdiam_real * start, int size,
+ int grid_size, int sample_size )
+{
+ gdiam_point * p_arr;
+ gdiam_bbox bb;
+
+ p_arr = gdiam_convert( start, size );
+ bb = gdiam_approx_mvbb_grid_sample( p_arr, size, grid_size, sample_size );
+ free( p_arr );
+
+ return bb;
+}
+
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic pop /* end ignoring warnings */
+#elif defined(__clang__)
+# pragma clang diagnostic pop /* end ignoring warnings */
+#endif
+
+/* gdiam.C - End of File ------------------------------------------*/
+
+
+// Local Variables:
+// tab-width: 8
+// mode: C++
+// c-basic-offset: 4
+// indent-tabs-mode: t
+// c-file-style: "stroustrup"
+// End:
+// ex: shiftwidth=4 tabstop=8
Copied: brlcad/trunk/src/librt/gdiam/gdiam.hpp (from rev 77314,
brlcad/trunk/src/other/libgdiam/gdiam.hpp)
===================================================================
--- brlcad/trunk/src/librt/gdiam/gdiam.hpp (rev 0)
+++ brlcad/trunk/src/librt/gdiam/gdiam.hpp 2020-10-01 21:41:50 UTC (rev
77315)
@@ -0,0 +1,727 @@
+/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
+ * gdiam.hpp -
+ * Implmeents an algorithm for computing a diameter,
+ *
+ * Copyright 2001 Sariel Har-Peled ([email protected])
+ * https://sarielhp.org/research/papers/00/diameter/diam_prog.html
+ *
+ * Note: gdiam is available under multiple licenses - we make
+ * use of it under the MIT license:
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
copy of
+ * this software and associated documentation files (the "Software"), to deal
in
+ * the Software without restriction, including without limitation the rights to
+ * use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies
+ * of the Software, and to permit persons to whom the Software is furnished to
do
+ * so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE
+ * SOFTWARE.
+ *
+ * Code is based on the paper:
+ * A Practical Approach for Computing the Diameter of a Point-Set.
+ * Sariel Har-Peled
+ \*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/
+
+#ifndef __GDIAM__H
+#define __GDIAM__H
+
+/* for g++ to quell warnings */
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic push /* start new diagnostic pragma */
+# pragma GCC diagnostic ignored "-Wfloat-equal"
+#elif defined(__clang__)
+# pragma clang diagnostic push /* start new diagnostic pragma */
+# pragma clang diagnostic ignored "-Wfloat-equal"
+#endif
+
+#define GDIAM_FLOAT_ZERO(val) (((val) > -1.0e-37) && ((val) < 1.0e-37))
+
+#define GDIAM_DIM 3
+typedef double gdiam_real;
+typedef gdiam_real gdiam_point_t[ GDIAM_DIM ];
+typedef gdiam_real gdiam_point_2d_t[ 2 ];
+typedef gdiam_real * gdiam_point_2d;
+typedef gdiam_real * gdiam_point;
+typedef const gdiam_real * gdiam_point_cnt;
+
+#ifndef __MINMAX_DEFINED
+#define __MINMAX_DEFINED
+
+#include<stdio.h>
+#include<cmath>
+#include<algorithm>
+#include<assert.h>
+#define min std::min
+#define max std::max
+
+#endif /* MIN_MAX */
+
+template <class T>
+inline void gdiam_exchange( T & a, T & b )
+{
+ T tmp = a;
+
+ a = b;
+ b = tmp;
+}
+
+inline gdiam_real pnt_length( const gdiam_point pnt )
+{
+ return sqrt( pnt[ 0 ] * pnt[ 0 ] + pnt[ 1 ] * pnt[ 1 ]
+ + pnt[ 2 ] * pnt[ 2 ] );
+}
+
+inline void pnt_normalize( gdiam_point pnt )
+{
+ gdiam_real len = pnt_length( pnt );
+ if ( GDIAM_FLOAT_ZERO(len) )
+ return;
+
+ pnt[ 0 ] /= len;
+ pnt[ 1 ] /= len;
+ pnt[ 2 ] /= len;
+}
+
+inline void pnt_copy( gdiam_point_t dest,
+ gdiam_point_t src )
+{
+ dest[ 0 ] = src[ 0 ];
+ dest[ 1 ] = src[ 1 ];
+ dest[ 2 ] = src[ 2 ];
+}
+inline void pnt_zero( gdiam_point dst ) {
+ dst[ 0 ] = dst[ 1 ] = dst[ 2 ] = 0;
+}
+inline void pnt_dump( gdiam_point_cnt pnt ) {
+ printf( "(%g, %g, %g)\n", pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+}
+
+
+inline gdiam_real pnt_dot_prod( gdiam_point_cnt a,
+ gdiam_point_cnt b )
+{
+ return a[ 0 ] * b[ 0 ]
+ + a[ 1 ] * b[ 1 ]
+ + a[ 2 ] * b[ 2 ];
+}
+
+inline void pnt_cross_prod( const gdiam_point a,
+ const gdiam_point b,
+ const gdiam_point out )
+{
+ out[ 0 ] = a[ 1 ] * b[ 2 ] - a[ 2 ] * b[ 1 ];
+ out[ 1 ] = - ( a[ 0 ] * b[ 2 ] - a[ 2 ] * b[ 0 ] );
+ out[ 2 ] = a[ 0 ] * b[ 1 ] - a[ 1 ] * b[ 0 ];
+}
+
+
+inline gdiam_real pnt_distance_2d( gdiam_point_2d p,
+ gdiam_point_2d q )
+{
+ gdiam_real valx = (p[ 0 ] - q[ 0 ]);
+ gdiam_real valy = (p[ 1 ] - q[ 1 ]);
+
+ return sqrt( valx * valx + valy * valy );
+}
+
+
+inline gdiam_real pnt_distance( gdiam_point p, gdiam_point q )
+{
+ gdiam_real valx = (p[ 0 ] - q[ 0 ]);
+ gdiam_real valy = (p[ 1 ] - q[ 1 ]);
+ gdiam_real valz = (p[ 2 ] - q[ 2 ]);
+
+ return sqrt( valx * valx + valy * valy + valz * valz );
+}
+
+inline gdiam_real pnt_distance( gdiam_point p, gdiam_point q,
+ gdiam_point_cnt dir )
+{
+ gdiam_real valx = (p[ 0 ] - q[ 0 ]);
+ gdiam_real valy = (p[ 1 ] - q[ 1 ]);
+ gdiam_real valz = (p[ 2 ] - q[ 2 ]);
+
+ gdiam_real len, proj_len;
+
+ len = sqrt( valx * valx + valy * valy + valz * valz );
+
+ proj_len = dir[ 0 ] * valx + dir[ 1 ] * valy + dir[ 2 ] * valz;
+
+ return sqrt( len * len - proj_len * proj_len );
+}
+
+inline void pnt_init( gdiam_point pnt,
+ gdiam_real x,
+ gdiam_real y,
+ gdiam_real z )
+{
+ pnt[ 0 ] = x;
+ pnt[ 1 ] = y;
+ pnt[ 2 ] = z;
+}
+
+
+inline void pnt_init_normalize( gdiam_point pnt,
+ gdiam_real x,
+ gdiam_real y,
+ gdiam_real z )
+{
+ pnt[ 0 ] = x;
+ pnt[ 1 ] = y;
+ pnt[ 2 ] = z;
+
+ pnt_normalize( pnt );
+}
+
+inline bool pnt_isEqual( const gdiam_point p,
+ const gdiam_point q )
+{
+ // Assuming here the GDIAM_DIM == 3 !!!!
+ return ( ( GDIAM_FLOAT_ZERO(p[0] - q[0]) )
+ && ( GDIAM_FLOAT_ZERO(p[0] - q[0]) )
+ && ( GDIAM_FLOAT_ZERO(p[0] - q[0]) ) );
+}
+
+inline void pnt_scale_and_add( gdiam_point dest,
+ gdiam_real coef,
+ gdiam_point_cnt vec ) {
+ dest[ 0 ] += coef * vec[ 0 ];
+ dest[ 1 ] += coef * vec[ 1 ];
+ dest[ 2 ] += coef * vec[ 2 ];
+}
+
+class GPointPair
+{
+ public:
+ gdiam_real distance;
+ gdiam_point p, q;
+
+ GPointPair() : p(), q() {
+ distance = 0.0;
+ }
+
+ void init( const gdiam_point _p,
+ const gdiam_point _q ) {
+ p = _p;
+ q = _q;
+ distance = pnt_distance( p, q );
+ }
+ void init( const gdiam_point _p,
+ const gdiam_point _q,
+ const gdiam_point proj ) {
+ p = _p;
+ q = _q;
+ distance = pnt_distance( p, q, proj );
+ }
+
+ void init( const gdiam_point pnt ) {
+ distance = 0;
+ p = q = pnt;
+ }
+
+ void update_diam_simple( const gdiam_point _p, const gdiam_point _q
) {
+ gdiam_real new_dist;
+
+ new_dist = pnt_distance( _p, _q );
+ if ( new_dist <= distance )
+ return;
+
+ //printf( "new_dist: %g\n", new_dist );
+ distance = new_dist;
+ p = _p;
+ q = _q;
+ }
+ void update_diam_simple( const gdiam_point _p, const gdiam_point _q,
+ const gdiam_point dir ) {
+ gdiam_real new_dist;
+
+ new_dist = pnt_distance( _p, _q, dir );
+ if ( new_dist <= distance )
+ return;
+
+ distance = new_dist;
+ p = _p;
+ q = _q;
+ }
+
+ void update_diam( const gdiam_point _p, const gdiam_point _q ) {
+ //update_diam_simple( p, _p );
+ //update_diam_simple( p, _q );
+ //update_diam_simple( q, _p );
+ //update_diam_simple( q, _q );
+ update_diam_simple( _p, _q );
+ }
+ void update_diam( const gdiam_point _p, const gdiam_point _q,
+ const gdiam_point dir ) {
+ //update_diam_simple( p, _p );
+ //update_diam_simple( p, _q );
+ //update_diam_simple( q, _p );
+ //update_diam_simple( q, _q );
+ update_diam_simple( _p, _q, dir );
+ }
+
+
+ void update_diam( GPointPair & pp ) {
+ //update_diam_simple( p, pp.p );
+ //update_diam_simple( p, pp.q );
+ //update_diam_simple( q, pp.p );
+ //update_diam_simple( q, pp.q );
+ update_diam_simple( pp.p, pp.q );
+ }
+};
+
+
+
+
+class GBBox {
+ private:
+ gdiam_real min_coords[ GDIAM_DIM ];
+ gdiam_real max_coords[ GDIAM_DIM ];
+
+ public:
+ void init( const GBBox & a,
+ const GBBox & b ) {
+ for (int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ min_coords[ ind ] = min( a.min_coords[ ind ],
+ b.min_coords[ ind ] );
+ max_coords[ ind ] = max( a.max_coords[ ind ],
+ b.max_coords[ ind ] );
+ }
+ }
+ void dump() const {
+ gdiam_real prod, diff;
+
+ prod = 1.0;
+ printf( "__________________________________________\n" );
+ for (int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ printf( "%d: [%g...%g]\n", ind, min_coords[ ind ],
+ max_coords[ ind ] );
+ diff = max_coords[ ind ] - min_coords[ ind ];
+ prod *= diff;
+ }
+ printf( "volume = %g\n", prod );
+ printf( "\\__________________________________________\n" );
+ }
+
+ void center( gdiam_point out ) const {
+ for (int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ out[ ind ] = ( min_coords[ ind ] + max_coords[ ind ] ) / 2.0;
+ }
+ }
+
+ gdiam_real volume() const {
+ gdiam_real prod, val;
+
+ prod = 1;
+ for ( int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ val = length_dim( ind );
+ prod *= val;
+ }
+
+ return prod;
+ }
+ void init() {
+ for (int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ min_coords[ ind ] = 1e20;
+ max_coords[ ind ] = -1e20;
+ }
+ }
+ /*
+ void dump() const {
+ printf( "---(" );
+ for (int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ printf( "[%g, %g] ",
+ min_coords[ ind ],
+ max_coords[ ind ] );
+ }
+ printf( ")\n" );
+ }*/
+
+ const gdiam_real & min_coord( int coord ) const {
+ return min_coords[ coord ];
+ }
+ const gdiam_real & max_coord( int coord ) const {
+ return max_coords[ coord ];
+ }
+
+ void bound( const gdiam_point pnt ) {
+ //cout << "bounding: " << pnt << "\n";
+ for ( int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ if ( pnt[ ind ] < min_coords[ ind ] )
+ min_coords[ ind ] = pnt[ ind ];
+ if ( pnt[ ind ] > max_coords[ ind ] )
+ max_coords[ ind ] = pnt[ ind ];
+ }
+ }
+
+ gdiam_real length_dim( int dim ) const {
+ return max_coords[ dim ] - min_coords[ dim ];
+ }
+ int getLongestDim() const {
+ int dim = 0;
+ gdiam_real len = length_dim( 0 );
+
+ for ( int ind = 1; ind < GDIAM_DIM; ind++ )
+ if ( length_dim( ind ) > len ) {
+ len = length_dim( ind );
+ dim = ind;
+ }
+
+ return dim;
+ }
+ gdiam_real getLongestEdge() const {
+ return length_dim( getLongestDim() );
+ }
+
+ gdiam_real get_diam() const
+ {
+ gdiam_real sum, val;
+
+ sum = 0;
+ for ( int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ val = length_dim( ind );
+ sum += val * val;
+ }
+
+ return sqrt( sum );
+ }
+
+ // in the following we assume that the length of dir is ONE!!!
+ // Note that the following is an overestaime - the diameter of
+ // projection of a cube, is bounded by the length of projections
+ // of its edges....
+ gdiam_real get_diam_proj( gdiam_point dir ) const
+ {
+ gdiam_real sum, coord;
+ gdiam_real prod, val;
+
+ //printf( "get_diam_proj: " );
+ sum = 0;
+ for ( int ind = 0; ind < GDIAM_DIM; ind++ ) {
+ coord = length_dim( ind );
+ //printf( "coord[%d]: %g\n",ind, coord );
+ prod = coord * dir [ ind ];
+ val = coord * coord - prod * prod;
+ assert( val >= 0 );
+ sum += sqrt( val );
+ }
+ //printf( "sum: %g, %g\n", sum, get_diam() );
+
+ // sum = squard diameter of the bounding box
+ // prod = length of projection of the diameter of cube on the
+ // direction.
+
+ return get_diam();
+ }
+
+ gdiam_real get_min_coord( int dim ) const {
+ return min_coords[ dim ];
+ }
+ gdiam_real get_max_coord( int dim ) const {
+ return max_coords[ dim ];
+ }
+};
+
+
+// gdiam_bbox is the famous - arbitrary orieted bounding box
+class gdiam_bbox
+{
+ private:
+ gdiam_point_t dir_1, dir_2, dir_3;
+
+ gdiam_real low_1, high_1;
+ gdiam_real low_2, high_2;
+ gdiam_real low_3, high_3;
+ bool f_init;
+
+ public:
+ gdiam_bbox() {
+ f_init = false;
+ }
+
+ gdiam_real volume() const {
+ gdiam_real len1, len2, len3;
+
+ len1 = ( high_1 - low_1 );
+ len2 = ( high_2 - low_2 );
+ len3 = ( high_3 - low_3 );
+
+ return len1 * len2 * len3;
+ }
+
+
+ void set_third_dim_longest() {
+ gdiam_point_t pnt_tmp;
+
+ if ( ( high_1 - low_1 ) > ( high_3 - low_3 ) ) {
+ gdiam_exchange( high_1, high_3 );
+ gdiam_exchange( low_1, low_3 );
+
+ memcpy( pnt_tmp, dir_1, sizeof( dir_1 ) );
+ //pnt_tmp = dir_1;
+ memcpy( dir_1, dir_3, sizeof( dir_3) );
+ //dir_1 = dir_3;
+ memcpy( dir_3, pnt_tmp, sizeof( dir_3 ) );
+ //dir_3 = pnt_tmp;
+ }
+ if ( ( high_2 - low_2 ) > ( high_3 - low_3 ) ) {
+ gdiam_exchange( high_2, high_3 );
+ gdiam_exchange( low_2, low_3 );
+
+ pnt_copy( pnt_tmp, dir_2 );
+ pnt_copy( dir_2, dir_3 );
+ pnt_copy( dir_3, pnt_tmp );
+ }
+ }
+
+ gdiam_point get_dir( int ind ) {
+ if ( ind == 0 )
+ return dir_1;
+ if ( ind == 1 )
+ return dir_2;
+ if ( ind == 2 )
+ return dir_3;
+
+ assert( false );
+
+ return NULL;
+ }
+
+ void combine( gdiam_point out,
+ double a_coef, double b_coef,
+ double c_coef ) {
+ for ( int ind = 0; ind < GDIAM_DIM; ind++ )
+ out[ ind ] = a_coef * dir_1[ ind ]
+ + b_coef * dir_2[ ind ]
+ + c_coef * dir_3[ ind ];
+ }
+
+ void get_vertex( double sel1, double sel2, double sel3, double *p_x,
double *p_y, double *p_z ) {
+ gdiam_real coef1, coef2, coef3;
+ gdiam_point_t pnt;
+
+ coef1 = low_1 + sel1 * ( high_1 - low_1 );
+ coef2 = low_2 + sel2 * ( high_2 - low_2 );
+ coef3 = low_3 + sel3 * ( high_3 - low_3 );
+
+ pnt_zero( pnt );
+ pnt_scale_and_add( pnt, coef1, dir_1 );
+ pnt_scale_and_add( pnt, coef2, dir_2 );
+ pnt_scale_and_add( pnt, coef3, dir_3 );
+
+ (*p_x) = pnt[0];
+ (*p_y) = pnt[1];
+ (*p_z) = pnt[2];
+
+ }
+
+
+ void dump_vertex( double sel1, double sel2, double sel3 ) const {
+ gdiam_real coef1, coef2, coef3;
+
+ //printf( "selection: (%g, %g, %g)\n", sel1, sel2, sel3 );
+ coef1 = low_1 + sel1 * ( high_1 - low_1 );
+ coef2 = low_2 + sel2 * ( high_2 - low_2 );
+ coef3 = low_3 + sel3 * ( high_3 - low_3 );
+
+ //printf( "coeficients: (%g, %g, %g)\n",
+ // coef1, coef2, coef3 );
+
+ gdiam_point_t pnt;
+ pnt_zero( pnt );
+
+ //printf( "starting...\n" );
+ //pnt_dump( pnt );
+ pnt_scale_and_add( pnt, coef1, dir_1 );
+ //pnt_dump( pnt );
+ pnt_scale_and_add( pnt, coef2, dir_2 );
+ //pnt_dump( pnt );
+ pnt_scale_and_add( pnt, coef3, dir_3 );
+ //pnt_dump( pnt );
+
+ pnt_dump( pnt );
+ }
+
+ void dump() const {
+ printf( "-----------------------------------------------\n" );
+ dump_vertex( 0, 0, 0 );
+ dump_vertex( 0, 0, 1 );
+ dump_vertex( 0, 1, 0 );
+ dump_vertex( 0, 1, 1 );
+ dump_vertex( 1, 0, 0 );
+ dump_vertex( 1, 0, 1 );
+ dump_vertex( 1, 1, 0 );
+ dump_vertex( 1, 1, 1 );
+ printf( "volume: %g\n", volume() );
+ printf( "Directions:\n" );
+ pnt_dump( dir_1 );
+ pnt_dump( dir_2 );
+ pnt_dump( dir_3 );
+ printf( "prods: %g, %g, %g\n",
+ pnt_dot_prod( dir_1, dir_2 ),
+ pnt_dot_prod( dir_1, dir_3 ),
+ pnt_dot_prod( dir_2, dir_3 ) );
+
+ printf( "--------------------------------------------------\n" );
+ //printf( "range_1: %g... %g\n", low_1, high_1 );
+ //printf( "range_2: %g... %g\n", low_2, high_2 );
+ //printf( "range_3: %g... %g\n", low_3, high_3 );
+ //printf( "prd: %g\n", pnt_dot_prod( dir_1, dir_2 ) );
+ //printf( "prd: %g\n", pnt_dot_prod( dir_2, dir_3 ) );
+ //printf( "prd: %g\n", pnt_dot_prod( dir_1, dir_3 ) );
+ }
+
+
+ void init( const GBBox & bb ) {
+ pnt_init( dir_1, 1, 0, 0 );
+ pnt_init( dir_2, 0, 1, 0 );
+ pnt_init( dir_3, 0, 0, 1 );
+
+ low_1 = bb.min_coord( 0 );
+ high_1 = bb.max_coord( 0 );
+
+ low_2 = bb.min_coord( 1 );
+ high_2 = bb.max_coord( 1 );
+
+ low_3 = bb.min_coord( 2 );
+ high_3 = bb.max_coord( 2 );
+ f_init = true;
+ }
+
+ void init( const gdiam_point _dir_1,
+ const gdiam_point _dir_2,
+ const gdiam_point _dir_3 ) {
+ memset( this, 0, sizeof( gdiam_bbox ) );
+
+ pnt_copy( dir_1, _dir_1 );
+ pnt_copy( dir_2, _dir_2 );
+ pnt_copy( dir_3, _dir_3 );
+ pnt_normalize( dir_1 );
+ pnt_normalize( dir_2 );
+ pnt_normalize( dir_3 );
+
+ if ( ( ! (fabs( pnt_dot_prod( dir_1, dir_2 ) ) < 1e-6 ) )
+ || ( ! (fabs( pnt_dot_prod( dir_1, dir_3 ) ) < 1e-6 ) )
+ || ( ! (fabs( pnt_dot_prod( dir_2, dir_3 ) ) < 1e-6 ) ) ) {
+ printf( "should be all close to zero: %g, %g, %g\n",
+ pnt_dot_prod( dir_1, dir_2 ),
+ pnt_dot_prod( dir_1, dir_3 ),
+ pnt_dot_prod( dir_2, dir_3 ) );
+ pnt_dump( _dir_1 );
+ pnt_dump( _dir_2 );
+ pnt_dump( _dir_3 );
+ fflush( stdout );
+ fflush( stderr );
+ assert( fabs( pnt_dot_prod( dir_1, dir_2 ) ) < 1e-6 );
+ assert( fabs( pnt_dot_prod( dir_1, dir_3 ) ) < 1e-6 );
+ assert( fabs( pnt_dot_prod( dir_2, dir_3 ) ) < 1e-6 );
+ }
+
+ // The following reduce the error by slightly improve the
@@ Diff output truncated at 100000 characters. @@
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
_______________________________________________
BRL-CAD Source Commits mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
