Changes merged, delete brep-debug (svn branch delete)

svn:revision:77232 svn:branch:brep-debug svn:account:starseeker
Add libged symbol for libtclcad's use of ged_dplot
2020-09-25 16:52:15 +00:00 · 2020-09-25 16:29:19 +00:00 · 2020-09-25 16:28:34 +00:00 · 2020-09-25 15:57:21 +00:00 · 2020-09-25 15:35:15 +00:00 · 2020-09-25 15:35:15 +00:00
39 changed files with 5647 additions and 124 deletions
@@ -13,6 +13,10 @@ changes made.  See document footer for additional details.
 --- 2020-xx-xx  Release 7.32.2                                     ---
 ----------------------------------------------------------------------

+* added subtract subcommand to MGED analyze command - Cliff Yapp
+* added intersect subcommand to MGED analyze command - Cliff Yapp
+* added summarize subcommand to MGED analyze command - Cliff Yapp
+* added quiet (-q) flag to MGED kill command - Cliff Yapp
 * added MGED stat command for reporting of object info - Cliff Yapp
 * corrected examples in gdiff man page - Cliff Yapp
 * fixed help message for -F option on gdiff command - Cliff Yapp
@@ -36,6 +36,7 @@ include(${BRLCAD_SOURCE_DIR}/misc/CMake/DocBook.cmake)

 add_subdirectory(articles)
 add_subdirectory(books)
+add_subdirectory(devguides)
 add_subdirectory(lessons)
 add_subdirectory(presentations)
 add_subdirectory(specifications)
@@ -0,0 +1,62 @@
+# Style sheet for XSLT transformation to HTML pages
+if(BRLCAD_EXTRADOCS_HTML)
+  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/../resources/brlcad/brlcad-man-xhtml-stylesheet.xsl.in
+    ${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/brlcad-man-xhtml-stylesheet.xsl)
+endif(BRLCAD_EXTRADOCS_HTML)
+
+if(BRLCAD_EXTRADOCS_PHP)
+  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/../resources/brlcad/brlcad-man-xhtml-stylesheet.xsl.in
+    ${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/brlcad-man-xhtml-stylesheet.xsl)
+endif(BRLCAD_EXTRADOCS_PHP)
+
+# Style sheet for XSLT transformation to manual pages
+if(BRLCAD_EXTRADOCS_MAN)
+  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/../resources/brlcad/brlcad-man-stylesheet.xsl.in
+    ${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/brlcad-man-stylesheet.xsl)
+endif(BRLCAD_EXTRADOCS_MAN)
+
+# Files for PDF
+if(BRLCAD_EXTRADOCS_PDF)
+  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/../resources/brlcad/brlcad-man-fo-stylesheet.xsl.in
+    ${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/brlcad-man-fo-stylesheet.xsl)
+endif(BRLCAD_EXTRADOCS_PDF)
+
+# For HTML, MAN and FO (FO is an intermediate file used in the
+# XML->PDF transformation) we use variables to hold the full
+# stylesheet path. In the case we need to further
+# customize FO stylesheets we can have separate CMake logic in
+# appropriate directories to handle the customization (e.g., the
+# BRL-CAD manuals in books/en/CMakeLists.txt).
+set(XSL_PHP_STYLESHEET "${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/wordpress.xsl")
+set(XSL_HTML_STYLESHEET "${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/brlcad-man-xhtml-stylesheet.xsl")
+set(XSL_MAN_STYLESHEET "${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/brlcad-man-stylesheet.xsl")
+set(XSL_FO_STYLESHEET "${CMAKE_CURRENT_BINARY_DIR}/../resources/brlcad/brlcad-man-fo-stylesheet.xsl")
+
+set(devguides_EN
+  bool_eval_development.xml
+  )
+
+set(devguides_EN_IMAGES
+  images/ccx_overlap_event.png
+  images/ssx_transverse_event.png
+  images/ssx_tangent_event.png
+  images/compare_endpoint_style.png
+  images/lcurves_with_shaded_context.png
+  images/ssx_10_vs_13.png
+  images/curve_traversal_directions.png
+  images/axis_X.png
+  images/intermediate_linked_curves.png
+  images/ssx_overlap_event.png
+  images/evaluation_overview.png
+  )
+
+ADD_DOC(devguides_EN html/devguides)
+ADD_DOC(devguides_EN_IMAGES html/devguides/images)
+ADD_DOCBOOK("HTML;PHP;PDF" devguides_EN devguides/en devguides_EN_IMAGES_cp)
+
+# Local Variables:
+# tab-width: 8
+# mode: cmake
+# indent-tabs-mode: t
+# End:
+# ex: shiftwidth=2 tabstop=8
@@ -10,9 +10,7 @@
  <refnamediv xml:id="name">
    <refname>analyze</refname>
    <refpurpose>
-      Displays the rotation and fallback angles, surface area, and
-      plane equation for each face of the ARB specified on the command line. The total
-      surface area, volume, and the length of each edge are also displayed.
+      Report and/or compute analytical information about the listed objects.
    </refpurpose>
  </refnamediv>

@@ -20,39 +18,243 @@
  <refsynopsisdiv xml:id="synopsis">
    <cmdsynopsis sepchar=" ">
      <command>analyze</command>
-      <arg choice="req" rep="norepeat"><replaceable>arb_name</replaceable></arg>
+      <arg choice="opt" rep="repeat"><replaceable>options</replaceable></arg>
+      <arg choice="opt" rep="norepeat"><replaceable>subcommand</replaceable></arg>
    </cmdsynopsis>
+    <variablelist>
+      
+      <varlistentry><term></term><listitem>
+      <cmdsynopsis sepchar=" ">
+	<command>summarize</command>
+	<arg choice="req" rep="repeat">obj</arg>
+      </cmdsynopsis>
+      </listitem></varlistentry>
+      
+      <varlistentry><term></term><listitem>
+      <cmdsynopsis sepchar=" ">
+	<command>intersect</command>
+	<arg choice="opt" rep="repeat"><replaceable>options</replaceable></arg>
+	<arg choice="req" rep="repeat">obj1 obj2 [...]</arg>
+      </cmdsynopsis>
+      </listitem></varlistentry>
+      
+      <varlistentry><term></term><listitem>
+      <cmdsynopsis sepchar=" ">
+	<command>subtract</command>
+	<arg choice="opt" rep="repeat"><replaceable>options</replaceable></arg>
+	<arg choice="req" rep="repeat">obj1 obj2 [...]</arg>
+      </cmdsynopsis>
+      </listitem></varlistentry>
+    </variablelist>
  </refsynopsisdiv>

  <refsection xml:id="description"><title>DESCRIPTION</title>

    <para>
-      The <command>analyze</command> command displays the rotation and fallback
-      angles, surface area, and plane equation for each face of the ARB specified
-      on the command line. The total surface area, volume, and the length of each
-      edge are also displayed. If executed while editing an ARB, the arb_name may
-      be omitted, and the ARB being edited will be analyzed.
+      The <command>analyze</command> command provides an array of reporting and
+      shape generation capabilities.  The classic behavior (and the behavior users
+      will see if calling without specifying a subcommand) is a report summarizing
+      geometric information about the specified solid(s).  This behavior may also
+      be explicitly invoked by the <command>summarize</command> subcommand.
    </para>
+
+    <para>
+      The <command>intersect</command> will take the first object specified and
+      intersect it with any subsequently specified objects - in effect, building
+      up a new object with the geoemtry common to all objects specified.
+    </para>
+
+    <para>
+      The <command>subtract</command> will take the first object specified and
+      subtract from it any subsequently specified objects.
+    </para>
+    
  </refsection>

  <refsection xml:id="examples"><title>EXAMPLES</title>

-    <para>
-      The example shows how to display information about an ARB's rotation and
-      fallback angles, surface area, and plane equation for its face.
-      Edge information is also displayed.
-    </para>
    <example><title>Analyze a particular <emphasis>ARB</emphasis></title>
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>analyze box.s</userinput>
+      </para>
+      <literallayout>
+box.s:  ARB8
+	1 (30.0283, -5.21153, -16.3791)
+	2 (30.0283, 21.5812, -16.3791)
+	3 (30.0283, 21.5812, 10.4137)
+	4 (30.0283, -5.21153, 10.4137)
+	5 (3.23558, -5.21153, -16.3791)
+	6 (3.23558, 21.5812, -16.3791)
+	7 (3.23558, 21.5812, 10.4137)
+	8 (3.23558, -5.21153, 10.4137)
+
+------+---------------------------+-------------------------------------------------+--------------+
+| FACE |     ROT          FB       |                 PLANE EQUATION                  | SURFACE AREA |
+------+---------------------------+-------------------------------------------------+--------------+
+| 1234 |   0.00000000   0.00000000 |  1.00000000  0.00000000  0.00000000 30.02833557 | 717.85172729 |
+| 5678 | 180.00000000  -0.00000000 | -1.00000000 -0.00000000 -0.00000000 -3.23558044 | 717.85172729 |
+| 1584 | 270.00000000  -0.00000000 | -0.00000000 -1.00000000 -0.00000000  5.21152973 | 717.85172729 |
+| 2376 |  90.00000000  -0.00000000 | -0.00000000  1.00000000 -0.00000000 21.58122540 | 717.85172729 |
+| 1265 |   0.00000000 -90.00000000 | -0.00000000  0.00000000 -1.00000000 16.37908936 | 717.85172729 |
+| 4378 |   0.00000000  90.00000000 |  0.00000000 -0.00000000  1.00000000 10.41366577 | 717.85172729 |
+------+---------------------------+-------------------------------------------------+--------------+
+
+  +--------------------+--------------------+--------------------+--------------------+
+  | EDGE        LENGTH | EDGE        LENGTH | EDGE        LENGTH | EDGE        LENGTH |
+  +--------------------+--------------------+--------------------+--------------------+
+  | 12     26.79275513 | 23     26.79275513 | 34     26.79275513 | 14     26.79275513 |
+  | 15     26.79275513 | 56     26.79275513 | 26     26.79275513 | 67     26.79275513 |
+  | 78     26.79275513 | 58     26.79275513 | 48     26.79275513 | 37     26.79275513 |
+  +--------------------+--------------------+--------------------+--------------------+
+
+    +-------------------------------+
+    | Volume       = 19233.22554681 |
+    | Surface Area =  4307.11036376 |
+    | Gallons      =     0.00508088 |
+    +-------------------------------+
+      </literallayout>
+    </example>
+
+     <example><title>Analyze a particular <emphasis>TOR</emphasis></title>
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>analyze tor</userinput>
+      </para>
+      <literallayout>
+tor:  torus (TOR)
+	V (4.91624, -32.8022, 31.7118), r1=25.4 (A), r2=5.08 (H)
+	N=(0, 1, 0)
+	A=(0, -0, 1)
+	B=(1, 0, -0)
+	vector to inner edge = (0, -0, 20.32)
+	vector to outer edge = (0, -0, 30.48)
+
+    Centroid: (4.91624, -32.8022, 31.7118)
+
+    +-------------------------------+
+    | Volume       = 12938.70529707 |
+    | Surface Area =  5093.97853992 |
+    | Gallons      =     0.00341804 |
+    +-------------------------------+
+      </literallayout>
+    </example>
+    
+    <example><title>Points inside a sphere</title>
+    
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>l sph.s</userinput>
+      </para>
+      <literallayout>
+sph.s:  ellipsoid (ELL)
+	V (0, 0, 0)
+	A (10, 0, 0) mag=10
+	B (0, 10, 0) mag=10
+	C (0, 0, 10) mag=10
+	A direction cosines=(0, 90, 90)
+	A rotation angle=0, fallback angle=0
+	B direction cosines=(90, 0, 90)
+	B rotation angle=90, fallback angle=0
+	C direction cosines=(90, 90, 0)
+	C rotation angle=0, fallback angle=90
+      </literallayout>
+      
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>l pnts.s</userinput>
+      </para>
+      <literallayout>
+pnts.s:  Point Cloud (PNTS)
+Total number of points: 6
+Default scale: 0.000000
+point#, (point)
+1, 	 (2.000000 2.000000 2.000000)
+2, 	 (0.000000 0.000000 -5.000000)
+3, 	 (0.000000 0.000000 10.000000)
+4, 	 (0.000000 0.000000 20.000000)
+5, 	 (10.000000 0.000000 10.000000)
+6, 	 (0.000000 10.000000 20.000000)
+      </literallayout>
+      
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>analyze intersect -o inside.s pnts.s sph.s</userinput>
+      </para>
+      <literallayout>
+3
+      </literallayout>
+      
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>l inside.s</userinput>
+      </para>
+      <literallayout>
+inside.s:  Point Cloud (PNTS)
+Total number of points: 3
+Default scale: 0.000000
+point#, (point)
+1, 	 (0.000000 0.000000 10.000000)
+2, 	 (0.000000 0.000000 -5.000000)
+3, 	 (2.000000 2.000000 2.000000)
+      </literallayout>
+
+    </example>
+
+
+    <example><title>Points outside a sphere</title>
+    
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>l sph.s</userinput>
+      </para>
+      <literallayout>
+sph.s:  ellipsoid (ELL)
+	V (0, 0, 0)
+	A (10, 0, 0) mag=10
+	B (0, 10, 0) mag=10
+	C (0, 0, 10) mag=10
+	A direction cosines=(0, 90, 90)
+	A rotation angle=0, fallback angle=0
+	B direction cosines=(90, 0, 90)
+	B rotation angle=90, fallback angle=0
+	C direction cosines=(90, 90, 0)
+	C rotation angle=0, fallback angle=90
+      </literallayout>
+      
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>l pnts.s</userinput>
+      </para>
+      <literallayout>
+pnts.s:  Point Cloud (PNTS)
+Total number of points: 6
+Default scale: 0.000000
+point#, (point)
+1, 	 (2.000000 2.000000 2.000000)
+2, 	 (0.000000 0.000000 -5.000000)
+3, 	 (0.000000 0.000000 10.000000)
+4, 	 (0.000000 0.000000 20.000000)
+5, 	 (10.000000 0.000000 10.000000)
+6, 	 (0.000000 10.000000 20.000000)
+      </literallayout>
+      
+      <para>
+	<prompt>mged&gt;</prompt> <userinput>analyze subtract -o outside.s pnts.s sph</userinput>
+      </para>
+      <literallayout>
+3
+      </literallayout>

      <para>
-	<prompt>mged&gt;</prompt> <userinput>analyze arb_name</userinput>
-      </para>
-      <para>
-	Rotation and fallback angles, surface area, and plane equation for
-	each face of the ARB specified on the command line are displayed along
-	with edge information.
+	<prompt>mged&gt;</prompt> <userinput>l outside.s</userinput>
      </para>
+      <literallayout>
+outside.s:  Point Cloud (PNTS)
+Total number of points: 3
+Default scale: 0.000000
+point#, (point)
+1, 	 (0.000000 10.000000 20.000000)
+2, 	 (10.000000 0.000000 10.000000)
+3, 	 (0.000000 0.000000 20.000000)
+      </literallayout>
+
+      
    </example>
+
+    
  </refsection>

  <refsection xml:id="author"><title>AUTHOR</title>
@@ -1,61 +1,79 @@
 <refentry xmlns="http://docbook.org/ns/docbook" version="5.0" xml:id="kill">

-<refmeta>
-  <refentrytitle>KILL</refentrytitle>
-  <manvolnum>nged</manvolnum>
-  <refmiscinfo class="source">BRL-CAD</refmiscinfo>
-  <refmiscinfo class="manual">BRL-CAD User Commands</refmiscinfo>
-</refmeta>
+  <refmeta>
+    <refentrytitle>KILL</refentrytitle>
+    <manvolnum>nged</manvolnum>
+    <refmiscinfo class="source">BRL-CAD</refmiscinfo>
+    <refmiscinfo class="manual">BRL-CAD User Commands</refmiscinfo>
+  </refmeta>

-<refnamediv xml:id="name">
-  <refname>kill</refname>
-  <refpurpose>Deletes specified <emphasis>objects</emphasis> from the current database.
-   </refpurpose>
-</refnamediv>
+  <refnamediv xml:id="name">
+    <refname>kill</refname>
+    <refpurpose>Deletes specified <emphasis>objects</emphasis> from the current database.
+    </refpurpose>
+  </refnamediv>

-<!-- body begins here -->
-<refsynopsisdiv xml:id="synopsis">
-  <cmdsynopsis sepchar=" ">
-    <command>kill</command>
-     <arg choice="opt" rep="norepeat">-f</arg>
-     <arg choice="req" rep="norepeat"><replaceable>objects</replaceable></arg>
-  </cmdsynopsis>
-</refsynopsisdiv>
+  <!-- body begins here -->
+  <refsynopsisdiv xml:id="synopsis">
+    <cmdsynopsis sepchar=" ">
+      <command>kill</command>
+      <arg choice="opt" rep="norepeat">-f</arg>
+      <arg choice="req" rep="norepeat"><replaceable>objects</replaceable></arg>
+    </cmdsynopsis>
+  </refsynopsisdiv>

-<refsection xml:id="description"><title>DESCRIPTION</title>
+  <refsection xml:id="description"><title>DESCRIPTION</title>

  <para>
-    Deletes the specified <emphasis>objects</emphasis> from the current database.
-	This command affects only the <emphasis>objects</emphasis> actually listed on the command 	line. If a combination is killed, its members are not affected. If the <emphasis>-f	</emphasis>option is specified, then kill will not complain if some, or all, of the 	<emphasis>objects</emphasis> specified do not actually exist in the database. Note that 	the <emphasis>objects</emphasis> are killed immediately. There is no need for a <command>	write file</command> command in MGED, and there is no <command>undo</command> command. 	<emphasis remap="B" role="bold">Use this command with caution.</emphasis> Other commands 	that remove objects from the database are <command>killall</command> and <command>killtree	</command>.
+    Deletes the specified <emphasis>objects</emphasis> from the current database.  This command
+    affects only the <emphasis>objects</emphasis> actually listed on the command line. If a
+    combination is killed, its members are not affected. If the <emphasis>-f</emphasis>option
+    is specified, then kill will not complain if some, or all, of the <emphasis>objects</emphasis>
+    specified do not actually exist in the database. Note that 	the <emphasis>objects</emphasis>
+    are killed immediately. There is no need for a <command>write file</command> command in MGED,
+    and there is no <command>undo</command> command. <emphasis remap="B" role="bold">Use this
+    command with caution.</emphasis> Other commands that remove objects from the database are
+    <command>killall</command> and <command>killtree</command>.
  </para>
-</refsection>
+  <para>
+    The <option>-q</option> option may be added to suppress database object lookup messages
+    reporting failure, which will otherwise let the user know an invalid object was specified.
+  </para>
+  </refsection>

-<refsection xml:id="examples"><title>EXAMPLES</title>
-
-  <para>The example shows the use of the <command>kill</command> command to delete specified 	objects from the database.
+  <refsection xml:id="examples"><title>EXAMPLES</title>
+  <para>
+    The example shows the use of the <command>kill</command> command to delete specified objects
+    from the database.
  </para>
  <example><title>Delete objects from the database.</title>

-    <variablelist>
-      <varlistentry>
-	   <term><prompt>mged&gt;</prompt> <userinput>kill group1 region2 shapeb</userinput></term>
-	   <listitem>
-	     <para>Deletes <emphasis>group1, region2</emphasis>, and <emphasis>shapeb</emphasis> 				from the database.
-	      </para>
-	   </listitem>
-      </varlistentry>
-    </variablelist>
+  <variablelist>
+    <varlistentry>
+      <term><prompt>mged&gt;</prompt> <userinput>kill group1 region2 shapeb</userinput></term>
+      <listitem>
+	<para>
+	  Deletes <emphasis>group1, region2</emphasis>, and <emphasis>shapeb</emphasis>	from
+	  the database.
+	</para>
+      </listitem>
+    </varlistentry>
+  </variablelist>
  </example>

-</refsection>
+  </refsection>

-<refsection xml:id="author"><title>AUTHOR</title><para>BRL-CAD Team</para></refsection>
+  <refsection xml:id="author"><title>AUTHOR</title>
+  <para>
+    BRL-CAD Team
+  </para>
+  </refsection>

-<refsection xml:id="bug_reports"><title>BUG REPORTS</title>
+  <refsection xml:id="bug_reports"><title>BUG REPORTS</title>

  <para>
    Reports of bugs or problems should be submitted via electronic
    mail to <email>devs@brlcad.org</email>
  </para>
-</refsection>
+  </refsection>
 </refentry>
@@ -860,6 +860,11 @@ GED_EXPORT int ged_pnts(struct ged *gedp, int argc, const char *argv[]);
 */
 GED_EXPORT extern int ged_stat(struct ged *gedp, int argc, const char *argv[]);

+
+
+/* Debugging command for brep plotting */
+GED_EXPORT extern int ged_dplot(struct ged *gedp, int argc, const char *argv[]);
+
 /** @} */


@@ -33,6 +33,7 @@ set(LIBBREP_SOURCES
  cdt/mesh.cpp
  cdt/ovlps_simple.cpp
  cdt/tri_isect.cpp
+  debug_plot.cpp
  intersect.cpp
  tools/tools.cpp
  opennurbs_ext.cpp
@@ -21,6 +21,8 @@
 *
 * Evaluate NURBS booleans (union, intersection and difference).
 *
+ * Additional documentation can be found in the "NURBS Boolean Evaluation
+ * Development Guide" docbook article (bool_eval_development.html).
 */

 #include "common.h"
@@ -31,6 +33,7 @@
 #include <queue>
 #include <set>
 #include <map>
+#include <sstream>

 #include "bio.h"

@@ -43,20 +46,15 @@
 #include "brep/ray.h"
 #include "brep/util.h"

+#include "debug_plot.h"
 #include "brep_except.h"
+#include "brep_defines.h"

+DebugPlot *dplot = NULL;

 // Whether to output the debug messages about b-rep booleans.
 #define DEBUG_BREP_BOOLEAN 0

-// tol value used in ON_Intersect()s. We use a smaller tolerance than the
-// default one 0.001.
-#define INTERSECTION_TOL 1e-4
-
-// tol value used in ON_3dVector::IsParallelTo(). We use a smaller tolerance
-// than the default one ON_PI/180.
-#define ANGLE_TOL ON_PI/1800.0
-

 struct IntersectPoint {
    ON_3dPoint m_pt;	// 3D intersection point
@@ -104,7 +102,7 @@ struct SSICurve {
 };


-HIDDEN void
+void
 append_to_polycurve(ON_Curve *curve, ON_PolyCurve &polycurve);
 // We link the SSICurves that share an endpoint, and form this new structure,
 // which has many similar behaviors as ON_Curve, e.g. PointAt(), Reverse().
@@ -344,7 +342,7 @@ curve_t_compare(const IntersectPoint *p1, const IntersectPoint *p2)
 }


-HIDDEN void
+void
 append_to_polycurve(ON_Curve *curve, ON_PolyCurve &polycurve)
 {
    // use this function rather than ON_PolyCurve::Append() to avoid
@@ -3445,7 +3443,7 @@ get_face_intersection_curves(
 	std::set<int> *unused = i < face_count1 ? &unused1 : &unused2;
 	std::set<int> *intact = i < face_count1 ? &finalform1 : &finalform2;
 	int curr_index = i < face_count1 ? i : i - face_count1;
-	if (face.BoundingBox().MinimumDistanceTo(brep->BoundingBox()) > ON_ZERO_TOLERANCE) {
+	if (face.BoundingBox().MinimumDistanceTo(brep->BoundingBox()) > INTERSECTION_TOL) {
 	    switch (operation) {
 		case BOOLEAN_UNION:
 		    intact->insert(curr_index);
@@ -3483,8 +3481,8 @@ get_face_intersection_curves(
 		if (unused2.find(j) == unused2.end() &&  finalform2.find(j) == finalform2.end()) {
 		    // If the two faces don't interact according to their bounding boxes,
 		    // they won't be a source of events - otherwise, they must be checked.
-		    fastf_t disjoint = brep1->m_F[i].BoundingBox().MinimumDistanceTo(brep2->m_F[j].BoundingBox());
-		    if (!(disjoint > ON_ZERO_TOLERANCE)) {
+		    fastf_t face_dist = brep1->m_F[i].BoundingBox().MinimumDistanceTo(brep2->m_F[j].BoundingBox());
+		    if (face_dist <= INTERSECTION_TOL) {
 			intersection_candidates.insert(std::pair<int, int>(i, j));
 		    }
 		}
@@ -3563,6 +3561,10 @@ get_face_intersection_curves(
 		if (results <= 0) {
 		    continue;
 		}
+
+		dplot->SSX(events, brep1, brep1->m_F[i].m_si, brep2, brep2->m_F[j].m_si);
+		dplot->WriteLog();
+
 		ON_SimpleArray<ON_Curve *> face1_curves, face2_curves;
 		for (int k = 0; k < events.Count(); k++) {
 		    if (events[k].m_type == ON_SSX_EVENT::ssx_tangent ||
@@ -3590,6 +3592,8 @@ get_face_intersection_curves(
 			}
 		    }
 		}
+		dplot->ClippedFaceCurves(surf1, surf2, face1_curves, face2_curves);
+		dplot->WriteLog();

 		if (DEBUG_BREP_BOOLEAN) {
 		    // Look for coplanar faces
@@ -3603,7 +3607,6 @@ get_face_intersection_curves(
 		}

 	    }
-
 	}
    }

@@ -3694,6 +3697,7 @@ categorize_trimmed_faces(
 	    }

 	    splitted[j]->m_rev = false;
+	    splitted[j]->m_belong_to_final = TrimmedFace::NOT_BELONG;
 	    switch (face_location) {
 		case INSIDE_BREP:
 		    if (operation == BOOLEAN_INTERSECT ||
@@ -3798,6 +3802,8 @@ get_evaluated_faces(const ON_Brep *brep1, const ON_Brep *brep2, op_type operatio
    for (int i = 0; i < original_faces.Count(); i++) {
 	TrimmedFace *first = original_faces[i];
 	ON_ClassArray<LinkedCurve> linked_curves = link_curves(curves_array[i]);
+	dplot->LinkedCurves(first->m_face->SurfaceOf(), linked_curves);
+	dplot->WriteLog();

 	ON_SimpleArray<TrimmedFace *> splitted = split_trimmed_face(first, linked_curves);
 	trimmed_faces.Append(splitted);
@@ -3826,6 +3832,9 @@ get_evaluated_faces(const ON_Brep *brep1, const ON_Brep *brep2, op_type operatio

    categorize_trimmed_faces(trimmed_faces, brep1, brep2, surf_tree1, surf_tree2, operation);

+    dplot->SplitFaces(trimmed_faces);
+    dplot->WriteLog();
+
    for (int i = 0; i < surf_tree1.Count(); i++) {
 	delete surf_tree1[i];
    }
@@ -3914,6 +3923,14 @@ standardize_loop_orientations(ON_Brep *brep)
 int
 ON_Boolean(ON_Brep *evaluated_brep, const ON_Brep *brep1, const ON_Brep *brep2, op_type operation)
 {
+    static int calls = 0;
+    ++calls;
+    std::ostringstream prefix;
+    prefix << "bool" << calls;
+    dplot = new DebugPlot(prefix.str().c_str());
+    dplot->Surfaces(brep1, brep2);
+    dplot->WriteLog();
+
    ON_ClassArray<ON_SimpleArray<TrimmedFace *> > trimmed_faces;
    try {
 	/* Deal with the trivial cases up front */
@@ -3934,14 +3951,17 @@ ON_Boolean(ON_Brep *evaluated_brep, const ON_Brep *brep1, const ON_Brep *brep2,
 	    }
 	    evaluated_brep->ShrinkSurfaces();
 	    evaluated_brep->Compact();
+	    dplot->WriteLog();
 	    return 0;
 	}
 	trimmed_faces = get_evaluated_faces(brep1, brep2, operation);
    } catch (InvalidBooleanOperation &e) {
 	bu_log("%s", e.what());
+	dplot->WriteLog();
 	return -1;
    } catch (GeometryGenerationError &e) {
 	bu_log("%s", e.what());
+	dplot->WriteLog();
 	return -1;
    }

@@ -3999,6 +4019,10 @@ ON_Boolean(ON_Brep *evaluated_brep, const ON_Brep *brep1, const ON_Brep *brep2,
 	bu_log("%s", ON_String(ws).Array());
    }

+    dplot->WriteLog();
+    delete dplot;
+    dplot = NULL;
+
    return 0;
 }

@@ -0,0 +1,71 @@
+/*                  B R E P _ D E F I N E S . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2014 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file brep_defines.h
+ *
+ * Private defines.
+ *
+ */
+
+// The maximal depth for subdivision - trade-off between accuracy and
+// performance.
+#define NR_MAX_DEPTH 8
+#define MAX_PCI_DEPTH NR_MAX_DEPTH
+#define MAX_PSI_DEPTH NR_MAX_DEPTH
+#define MAX_CCI_DEPTH NR_MAX_DEPTH
+#define MAX_CSI_DEPTH NR_MAX_DEPTH
+#define MAX_SSI_DEPTH NR_MAX_DEPTH
+
+// Used to prevent an infinite loop in the unlikely event that we
+// can't provide a good starting point for the Newton-Raphson
+// Iteration.
+#define NR_MAX_ITERATIONS 100
+#define PCI_MAX_ITERATIONS NR_MAX_ITERATIONS
+#define PSI_MAX_ITERATIONS NR_MAX_ITERATIONS
+#define CCI_MAX_ITERATIONS NR_MAX_ITERATIONS
+#define CSI_MAX_ITERATIONS NR_MAX_ITERATIONS
+#define SSI_MAX_ITERATIONS NR_MAX_ITERATIONS
+
+// We make the default tolerance for PSI the same as that of curve and
+// surface intersections defined by openNURBS (see opennurbs_curve.h
+// and opennurbs_surface.h).
+#define NR_DEFAULT_TOLERANCE 0.001
+#define PCI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
+#define PSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
+#define CCI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
+#define CSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
+#define SSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
+
+// tol value used in ON_Intersect()s. We use a smaller tolerance than the
+// default one 0.001.
+#define INTERSECTION_TOL 1e-4
+
+// tol value used in ON_3dVector::IsParallelTo(). We use a smaller tolerance
+// than the default one ON_PI/180.
+#define ANGLE_TOL ON_PI/1800.0
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -0,0 +1,990 @@
+/*                  D E B U G _ P L O T . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2014 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file debug_plot.cpp
+ *
+ * DebugPlot implementation. Currently borrows code from librt to
+ * handle the creation of vlists for brep geometry and conversion of
+ * those vlists to unix plot files.
+ */
+
+#include <iostream>
+#include <sstream>
+#include <string>
+#include "bu/log.h"
+#include "bn.h"
+#include "raytrace.h"
+#include "vmath.h"
+#include "debug_plot.h"
+#include "brep_except.h"
+
+static unsigned char surface1_color[] = {0, 0, 62};
+static unsigned char surface2_color[] = {62, 0, 0};
+static unsigned char surface1_highlight_color[] = {56, 56, 255};
+static unsigned char surface2_highlight_color[] = {255, 56, 56};
+
+static unsigned char tangent_color[] = {255, 255, 255};
+static unsigned char transverse_color[] = {255, 255, 0};
+static unsigned char overlap_color[] = {0, 255, 0};
+
+static unsigned char accepted_outerloop_color[] = {0, 255, 0};
+static unsigned char accepted_innerloop_color[] = {255, 0, 0};
+static unsigned char unknown_outerloop_color[] = {158, 158, 0};
+static unsigned char unknown_innerloop_color[] = {158, 158, 0};
+static unsigned char rejected_outerloop_color[] = {0, 62, 0};
+static unsigned char rejected_innerloop_color[] = {62, 0, 0};
+
+DebugPlot::DebugPlot(const char *basename) :
+    prefix(basename),
+    have_surfaces(false),
+    brep1_surf_count(0),
+    brep2_surf_count(0),
+    linked_curve_count(0)
+{
+    BU_LIST_INIT(&vlist_free_list);
+}
+
+DebugPlot::~DebugPlot()
+{
+    struct bn_vlist *vp;
+    while (BU_LIST_WHILE(vp, bn_vlist, &vlist_free_list)) {
+	BU_LIST_DEQUEUE(&(vp->l));
+	bu_free((char *)vp, "bn_vlist");
+    }
+}
+
+int
+DebugPlot::SurfacePairs(void)
+{
+    return (int)intersecting_surfaces.size();
+}
+
+int
+DebugPlot::IntersectingIsocurves(int ssx_idx)
+{
+    int max_isocsx_idx = (int)ssx_isocsx_events.size() - 1;
+    if (ssx_idx < 0 || ssx_idx > max_isocsx_idx) {
+	std::cerr << "DebugPlot::IntersectingIsocurves passed invalid ssx index.\n";
+	return 0;
+    }
+    return (int)ssx_isocsx_events[ssx_idx].size();
+}
+
+int
+DebugPlot::LinkedCurves(void)
+{
+    return linked_curve_count;
+}
+
+HIDDEN void
+rt_vlist_to_uplot(FILE *fp, const struct bu_list *vhead)
+{
+    struct bn_vlist *vp;
+
+    for (BU_LIST_FOR(vp, bn_vlist, vhead)) {
+	 int i;
+	 int nused = vp->nused;
+	 const int *cmd = vp->cmd;
+	 point_t *pt = vp->pt;
+
+	for (i = 0; i < nused; i++, cmd++, pt++) {
+	    switch (*cmd) {
+		case BN_VLIST_POLY_START:
+		case BN_VLIST_TRI_START:
+		    break;
+		case BN_VLIST_POLY_MOVE:
+		case BN_VLIST_LINE_MOVE:
+		case BN_VLIST_TRI_MOVE:
+		    pdv_3move(fp, *pt);
+		    break;
+		case BN_VLIST_POLY_DRAW:
+		case BN_VLIST_POLY_END:
+		case BN_VLIST_LINE_DRAW:
+		case BN_VLIST_TRI_DRAW:
+		case BN_VLIST_TRI_END:
+		    pdv_3cont(fp, *pt);
+		    break;
+		default:
+		    bu_log("rt_vlist_to_uplot: unknown vlist cmd x%x\n",
+			   *cmd);
+	    }
+	}
+    }
+}
+
+HIDDEN void
+write_plot_to_file(
+    const char *filename,
+    const struct bu_list *vhead,
+    const unsigned char *color)
+{
+    FILE *fp = fopen(filename, "w");
+
+    if (!color) {
+	unsigned char clr[] = {255, 0, 0};
+	color = clr;
+    }
+
+    pl_linmod(fp, "solid");
+    pl_color(fp, color[0], color[1], color[2]);
+    rt_vlist_to_uplot(fp, vhead);
+
+    fclose(fp);
+}
+
+void
+DebugPlot::WriteLog()
+{
+    std::ostringstream filename;
+
+    // First write out empty plots of different colors, we use this to
+    // "erase" unwanted overlays when running the dplot command.
+    struct bu_list vhead;
+    BU_LIST_INIT(&vhead);
+    point_t origin = {0.0, 0.0, 0.0};
+    BN_ADD_VLIST(&vlist_free_list, &vhead, origin, BN_VLIST_LINE_MOVE);
+
+    filename << prefix << "_empty0.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, tangent_color);
+    filename.str("");
+    filename << prefix << "_empty1.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, transverse_color);
+    filename.str("");
+    filename << prefix << "_empty2.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, overlap_color);
+    filename.str("");
+    filename << prefix << "_empty3.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, surface1_color);
+    filename.str("");
+    filename << prefix << "_empty4.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, surface2_color);
+    filename.str("");
+    filename << prefix << "_empty5.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, surface1_highlight_color);
+    filename.str("");
+    filename << prefix << "_empty6.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, surface2_highlight_color);
+    filename.str("");
+    filename << prefix << "_empty7.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, accepted_outerloop_color);
+    filename.str("");
+    filename << prefix << "_empty8.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, accepted_innerloop_color);
+    filename.str("");
+    filename << prefix << "_empty9.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, unknown_outerloop_color);
+    filename.str("");
+    filename << prefix << "_empty10.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, unknown_innerloop_color);
+    filename.str("");
+    filename << prefix << "_empty11.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, rejected_outerloop_color);
+    filename.str("");
+    filename << prefix << "_empty12.plot3";
+    write_plot_to_file(filename.str().c_str(), &vhead, rejected_innerloop_color);
+
+    // create dplot log file
+    filename.str("");
+    filename << prefix << ".dplot";
+    FILE *fp = fopen(filename.str().c_str(), "w");
+
+    // write out surface-surface intersections
+    fprintf(fp, "surfaces %d %d\n", brep1_surf_count, brep2_surf_count);
+    for (size_t i = 0; i < intersecting_surfaces.size(); ++i) {
+	std::pair<int, int> intersecting = intersecting_surfaces[i];
+
+	int events, b1_isocurves, intersecting_isocurves, b1_clipped, b2_clipped;
+
+	try {
+	    events = ssx_events.at(i);
+	} catch (std::out_of_range &e) {
+	    events = 0;
+	}
+	try {
+	    std::pair<int, int> ccount = ssx_clipped_curves.at(i);
+	    b1_clipped = ccount.first;
+	    b2_clipped = ccount.second;
+	} catch (std::out_of_range &e) {
+	    b1_clipped = b2_clipped = 0;
+	}
+	try {
+	    b1_isocurves = ssx_isocsx_brep1_curves.at(i);
+	} catch (std::out_of_range &e) {
+	    b1_isocurves = 0;
+	}
+	try {
+	    intersecting_isocurves = (int)ssx_isocsx_events.at(i).size();
+	} catch (std::out_of_range &e) {
+	    intersecting_isocurves = 0;
+	}
+
+	// b1si b2si finalevents b1ccurves b2ccurves b1_isocurve_xs total_isocurve_xs isocsx0_event0 ...
+	fprintf(fp, "ssx %d %d %d %d %d %d %d", intersecting.first,
+		intersecting.second, events, b1_clipped, b2_clipped,
+		b1_isocurves, intersecting_isocurves);
+
+	if (ssx_isocsx_events.size() > i) {
+	    for (size_t j = 0; j < ssx_isocsx_events[i].size(); ++j) {
+		fprintf(fp, " %d", ssx_isocsx_events[i][j]);
+	    }
+	}
+	fprintf(fp, "\n");
+    }
+
+    // write out linked curve count
+    if (linked_curve_count > 0) {
+	fprintf(fp, "linkedcurves %d\n", linked_curve_count);
+    }
+
+    // write out split faces
+    size_t split_faces = split_face_outerloop_curves.size();
+    if (split_faces > 0) {
+	fprintf(fp, "splitfaces %d\n", (int)split_faces);
+    }
+    for (size_t i = 0; i < split_faces; ++i) {
+	fprintf(fp, "splitface %d %d %d\n", (int)i,
+		split_face_outerloop_curves[i],
+		split_face_innerloop_curves[i]);
+    }
+    fclose(fp);
+    BN_FREE_VLIST(&vlist_free_list, &vhead);
+}
+
+HIDDEN double
+find_next_t(const ON_Curve* crv, double start_t, double step, double max_dist)
+{
+    ON_Interval dom = crv->Domain();
+    ON_3dPoint prev_pt = crv->PointAt(dom.ParameterAt(start_t));
+    ON_3dPoint next_pt;
+
+    // ensure that (start + step) < 1.0
+    if (start_t + step > 1.0) {
+	step = 1.0 - start_t - BN_TOL_DIST;
+    }
+    
+    // reduce step until next point is within tolerance
+    while (step > BN_TOL_DIST) {
+	next_pt = crv->PointAt(dom.ParameterAt(start_t + step));
+
+	if (prev_pt.DistanceTo(next_pt) <= max_dist) {
+	    return start_t + step;
+	}
+	step /= 2.0;
+    }
+    // if we couldn't find a point within tolerance, give up and jump
+    // to end of domain
+    return 1.0;
+}
+
+void
+DebugPlot::Plot3DCurve(
+    const ON_Curve *crv,
+    const char *filename,
+    unsigned char *color,
+    struct bu_list *vlist /* = NULL */)
+{
+    struct bu_list vhead_tmp;
+    BU_LIST_INIT(&vhead_tmp);
+
+    struct bu_list *vhead = &vhead_tmp;
+    if (vlist) {
+	vhead = vlist;
+    }
+
+    ON_Interval crv_dom = crv->Domain();
+
+    // Insert first point.
+    point_t pt1;
+    ON_3dPoint p;
+    p = crv->PointAt(crv_dom.ParameterAt(0.0));
+    VMOVE(pt1, p);
+    BN_ADD_VLIST(&vlist_free_list, vhead, pt1, BN_VLIST_LINE_MOVE);
+
+    /* Dynamic sampling approach - start with an initial guess
+     * for the next point of one tenth of the domain length
+     * further down the domain from the previous value.  Set a
+     * maximum physical distance between points of 100 times
+     * the model tolerance.  Reduce the increment until the
+     * tolerance is satisfied, then add the point and use it
+     * as the starting point for the next calculation until
+     * the whole domain is finished.  Perhaps it would be more
+     * ideal to base the tolerance on some fraction of the
+     * curve bounding box dimensions?
+     */
+    double t = 0.0;
+    while (t < 1.0) {
+	t = find_next_t(crv, t, 0.1, BN_TOL_DIST * 100);
+	p = crv->PointAt(crv_dom.ParameterAt(t));
+	VMOVE(pt1, p);
+
+	BN_ADD_VLIST(&vlist_free_list, vhead, pt1, BN_VLIST_LINE_DRAW);
+    }
+
+    if (!vlist) {
+	write_plot_to_file(filename, vhead, color);
+	BN_FREE_VLIST(&vlist_free_list, vhead);
+    }
+}
+
+void
+DebugPlot::Plot3DCurveFrom2D(
+    const ON_Surface *surf,
+    const ON_Curve *crv,
+    const char *filename,
+    unsigned char *color,
+    bool decorate /* = false */)
+{
+    struct bu_list vhead;
+    BU_LIST_INIT(&vhead);
+
+    ON_Interval crv_dom = crv->Domain();
+
+    ON_3dPoint p, uv;
+
+    // Insert first point.
+    point_t pt1, first_pt, last_pt, prev_pt;
+    ON_3dVector normal;
+    uv = crv->PointAt(crv_dom.ParameterAt(0.0));
+    surf->EvNormal(uv.x, uv.y, p, normal);
+    VMOVE(first_pt, p);
+
+    uv = crv->PointAt(crv_dom.ParameterAt(1.0));
+    surf->EvNormal(uv.x, uv.y, p, normal);
+    VMOVE(last_pt, p);
+
+    bool closed = false;
+    if (VNEAR_EQUAL(first_pt, last_pt, BN_TOL_DIST)) {
+	closed = true;
+    }
+
+    VMOVE(pt1, first_pt);
+    BN_ADD_VLIST(&vlist_free_list, &vhead, pt1, BN_VLIST_LINE_MOVE);
+
+    /* Dynamic sampling approach - start with an initial guess
+     * for the next point of one tenth of the domain length
+     * further down the domain from the previous value.  Set a
+     * maximum physical distance between points of 100 times
+     * the model tolerance.  Reduce the increment until the
+     * tolerance is satisfied, then add the point and use it
+     * as the starting point for the next calculation until
+     * the whole domain is finished.  Perhaps it would be more
+     * ideal to base the tolerance on some fraction of the
+     * curve bounding box dimensions?
+     */
+    double t = 0.0;
+    bool first = true;
+    double min_mag = BN_TOL_DIST * 10.0;
+    double mag_tan = min_mag;
+    vect_t tangent = {0.0, 0.0, min_mag};
+    vect_t perp, barb;
+    while (t < 1.0) {
+	t = find_next_t(crv, t, 0.1, BN_TOL_DIST * 100);
+
+	uv = crv->PointAt(crv_dom.ParameterAt(t));
+	surf->EvNormal(uv.x, uv.y, p, normal);
+	VMOVE(prev_pt, pt1);
+	VMOVE(pt1, p);
+
+	if (first || !VNEAR_EQUAL(pt1, prev_pt, BN_TOL_DIST)) {
+	    VSUB2(tangent, pt1, prev_pt);
+	}
+
+	if (decorate && first) {
+	    first = false;
+
+	    mag_tan = DIST_PNT_PNT(prev_pt, pt1);
+	    mag_tan = FMAX(mag_tan, min_mag);
+
+	    VUNITIZE(tangent);
+	    VCROSS(perp, tangent, normal);
+
+	    if (!closed) {
+		VSCALE(tangent, tangent, mag_tan);
+		VUNITIZE(perp);
+		VSCALE(perp, perp, mag_tan);
+
+		VADD3(barb, prev_pt, tangent, perp);
+		BN_ADD_VLIST(&vlist_free_list, &vhead, barb, BN_VLIST_LINE_DRAW);
+		BN_ADD_VLIST(&vlist_free_list, &vhead, prev_pt, BN_VLIST_LINE_MOVE);
+
+		VSCALE(perp, perp, -1.0);
+		VADD3(barb, prev_pt, tangent, perp);
+		BN_ADD_VLIST(&vlist_free_list, &vhead, barb, BN_VLIST_LINE_DRAW);
+		BN_ADD_VLIST(&vlist_free_list, &vhead, prev_pt, BN_VLIST_LINE_MOVE);
+	    }
+	}
+	BN_ADD_VLIST(&vlist_free_list, &vhead, pt1, BN_VLIST_LINE_DRAW);
+    }
+    if (decorate) {
+	VUNITIZE(tangent);
+	VSCALE(tangent, tangent, -mag_tan);
+
+	VCROSS(perp, tangent, normal);
+	VUNITIZE(perp);
+	VSCALE(perp, perp, mag_tan);
+
+	VADD2(barb, pt1, perp);
+	if (!closed) {
+	    VADD2(barb, barb, tangent);
+	}
+	BN_ADD_VLIST(&vlist_free_list, &vhead, barb, BN_VLIST_LINE_DRAW);
+	BN_ADD_VLIST(&vlist_free_list, &vhead, pt1, BN_VLIST_LINE_MOVE);
+
+	VSCALE(perp, perp, -1.0);
+	VADD2(barb, pt1, perp);
+	if (!closed) {
+	    VADD2(barb, barb, tangent);
+	}
+	BN_ADD_VLIST(&vlist_free_list, &vhead, barb, BN_VLIST_LINE_DRAW);
+    }
+
+    write_plot_to_file(filename, &vhead, color);
+    BN_FREE_VLIST(&vlist_free_list, &vhead);
+}
+
+void
+DebugPlot::PlotBoundaryIsocurves(
+    struct bu_list *vlist,
+    const ON_Surface &surf,
+    int knot_dir)
+{
+    int surf_dir = 1 - knot_dir;
+    int knot_count = surf.SpanCount(surf_dir) + 1;
+
+    double *surf_knots = new double[knot_count];
+    surf.GetSpanVector(surf_dir, surf_knots);
+
+    // knots that can be boundaries of Bezier patches
+    ON_SimpleArray<double> surf_bknots;
+    surf_bknots.Append(surf_knots[0]);
+    for (int i = 1; i < knot_count; i++) {
+	if (surf_knots[i] > *(surf_bknots.Last())) {
+	    surf_bknots.Append(surf_knots[i]);
+	}
+    }
+    delete[] surf_knots;
+
+    if (surf.IsClosed(surf_dir)) {
+	surf_bknots.Remove();
+    }
+
+    for (int i = 0; i < surf_bknots.Count(); i++) {
+	ON_Curve *surf_boundary_iso = surf.IsoCurve(knot_dir, surf_bknots[i]);
+	Plot3DCurve(surf_boundary_iso, NULL, NULL, vlist);
+	delete surf_boundary_iso;
+    }
+}
+
+void
+DebugPlot::PlotSurface(
+    const ON_Surface &surf,
+    const char *filename,
+    unsigned char *color)
+{
+    struct bu_list vhead;
+    BU_LIST_INIT(&vhead);
+
+    PlotBoundaryIsocurves(&vhead, surf, 0);
+    PlotBoundaryIsocurves(&vhead, surf, 1);
+
+    write_plot_to_file(filename, &vhead, color);
+
+    BN_FREE_VLIST(&vlist_free_list, &vhead);
+}
+
+
+void
+DebugPlot::Surfaces(const ON_Brep *brep1, const ON_Brep *brep2)
+{
+    if (!brep1 || !brep2) {
+	std::cerr << "error: dplot_surfaces: NULL args\n";
+	return;
+    }
+
+    brep1_surf_count = brep1->m_S.Count();
+    for (int i = 0; i < brep1->m_S.Count(); i++) {
+	ON_Surface *surf = brep1->m_S[i];
+	std::ostringstream filename;
+	filename << prefix << "_brep1_surface" << i << ".plot3";
+	PlotSurface(*surf, filename.str().c_str(), surface1_color);
+    }
+
+    brep2_surf_count = brep2->m_S.Count();
+    for (int i = 0; i < brep2->m_S.Count(); i++) {
+	ON_Surface *surf = brep2->m_S[i];
+	std::ostringstream filename;
+	filename << prefix << "_brep2_surface" << i << ".plot3";
+	PlotSurface(*surf, filename.str().c_str(), surface2_color);
+    }
+    have_surfaces = true;
+}
+
+int get_subcurve_inside_faces(const ON_Brep *brep1, const ON_Brep *brep2, int face_i1, int face_i2, ON_SSX_EVENT *event);
+
+void
+DebugPlot::SSX(
+    const ON_ClassArray<ON_SSX_EVENT> &events,
+    const ON_Brep *brep1, int brep1_surf,
+    const ON_Brep *brep2, int brep2_surf)
+{
+    ON_Surface *surf;
+    std::ostringstream filename;
+
+    // create highlighted plot of brep1 surface if it doesn't exist
+    filename << prefix << "_highlight_brep1_surface" << brep1_surf << ".plot3";
+    if (!bu_file_exists(filename.str().c_str(), NULL)) {
+	surf = brep1->m_S[brep1_surf];
+	PlotSurface(*surf, filename.str().c_str(), surface1_highlight_color);
+    }
+
+    // create highlighted plot of brep2 surface if it doesn't exist
+    filename.str("");
+    filename << prefix << "_highlight_brep2_surface" << brep2_surf << ".plot3";
+    if (!bu_file_exists(filename.str().c_str(), NULL)) {
+	surf = brep2->m_S[brep2_surf];
+	PlotSurface(*surf, filename.str().c_str(), surface2_highlight_color);
+    }
+
+    // create plot of the intersections between these surfaces
+    surf = brep1->m_S[brep1_surf];
+    size_t ssx_idx = intersecting_surfaces.size();
+    int plot_count = 0;
+    for (int i = 0; i < events.Count(); ++i) {
+	filename.str("");
+	filename << prefix << "_ssx" << ssx_idx << "_event" << plot_count <<
+	    ".plot3";
+
+	if (events[i].m_type == ON_SSX_EVENT::ssx_tangent) {
+	    Plot3DCurveFrom2D(surf, events[i].m_curveA,
+		    filename.str().c_str(), tangent_color, true);
+	    ++plot_count;
+	} else if (events[i].m_type == ON_SSX_EVENT::ssx_transverse) {
+	    Plot3DCurveFrom2D(surf, events[i].m_curveA,
+		    filename.str().c_str(), transverse_color, true);
+	    ++plot_count;
+	} else if (events[i].m_type == ON_SSX_EVENT::ssx_overlap) {
+	    Plot3DCurveFrom2D(surf, events[i].m_curveA,
+		    filename.str().c_str(), overlap_color, true);
+	    ++plot_count;
+	}
+    }
+    // stash surface indices and event count
+    std::pair<int, int> ssx_pair(brep1_surf, brep2_surf);
+    intersecting_surfaces.push_back(ssx_pair);
+    ssx_events.push_back(plot_count);
+}
+
+void
+DebugPlot::IsoCSX(
+    const ON_SimpleArray<ON_X_EVENT> &events,
+    const ON_Curve *isocurve,
+    bool is_brep1_iso) // is the isocurve from brep1?
+{
+    size_t ssx_idx = intersecting_surfaces.size();
+
+    // create plot of the intersections between the curve and surface
+    while (ssx_isocsx_events.size() < (ssx_idx + 1)) {
+	ssx_isocsx_events.push_back(std::vector<int>());
+    }
+    size_t isocsx_idx = ssx_isocsx_events[ssx_idx].size();
+    int plot_count = 0;
+    for (int i = 0; i < events.Count(); ++i) {
+	if (events[i].m_type == ON_X_EVENT::csx_overlap) {
+	    std::ostringstream filename;
+	    filename << prefix << "_ssx" << ssx_idx << "_isocsx" << isocsx_idx
+		<< "_event" << plot_count++ << ".plot3";
+
+	    try {
+		ON_Curve *event_curve = sub_curve(isocurve, events[i].m_a[0],
+			events[i].m_a[1]);
+		Plot3DCurve(event_curve, filename.str().c_str(), overlap_color);
+	    } catch (InvalidInterval &e) {
+		std::cerr << "error: IsoCSX event contains degenerate interval\n";
+	    }
+	}
+    }
+    if (plot_count) {
+	// create highlighted plot of isocurve if it doesn't already exist
+	std::ostringstream filename;
+	filename << prefix << "_highlight_ssx" << ssx_idx << "_isocurve" <<
+	    isocsx_idx << ".plot3";
+	if (!bu_file_exists(filename.str().c_str(), NULL)) {
+	    if (is_brep1_iso) {
+		Plot3DCurve(isocurve, filename.str().c_str(),
+			surface1_highlight_color);
+	    } else {
+		Plot3DCurve(isocurve, filename.str().c_str(),
+			surface2_highlight_color);
+	    }
+	}
+
+	// remember event count for this isocsx
+	ssx_isocsx_events[ssx_idx].push_back(plot_count);
+
+	// remember how many events are for brep1 isocurve and brep2 surface,
+	if (is_brep1_iso) {
+	    while (ssx_isocsx_brep1_curves.size() < (ssx_idx + 1)) {
+		ssx_isocsx_brep1_curves.push_back(0);
+	    }
+	    ++ssx_isocsx_brep1_curves[ssx_idx];
+	}
+    }
+}
+
+void
+DebugPlot::ClippedFaceCurves(
+    const ON_Surface *surf1,
+    const ON_Surface *surf2,
+    const ON_SimpleArray<ON_Curve *> &face1_curves,
+    const ON_SimpleArray<ON_Curve *> &face2_curves)
+{
+    // plot clipped tangent/transverse/overlap curves
+    size_t ssx_idx = intersecting_surfaces.size() - 1;
+    for (int i = 0; i < face1_curves.Count(); ++i) {
+	std::ostringstream filename;
+	filename << prefix << "_ssx" << ssx_idx << "_brep1face_clipped_curve" << i << ".plot3";
+	Plot3DCurveFrom2D(surf1, face1_curves[i], filename.str().c_str(),
+		surface1_highlight_color, true);
+    }
+    for (int i = 0; i < face2_curves.Count(); ++i) {
+	std::ostringstream filename;
+	filename << prefix << "_ssx" << ssx_idx << "_brep2face_clipped_curve" << i << ".plot3";
+	Plot3DCurveFrom2D(surf2, face2_curves[i], filename.str().c_str(),
+		surface2_highlight_color, true);
+    }
+
+    while (ssx_clipped_curves.size() < (ssx_idx + 1)) {
+	ssx_clipped_curves.push_back(std::pair<int, int>(0, 0));
+    }
+    std::pair<int, int> counts(face1_curves.Count(), face2_curves.Count());
+    ssx_clipped_curves[ssx_idx] = counts;
+}
+
+struct TrimmedFace {
+    // curve segments in the face's outer loop
+    ON_SimpleArray<ON_Curve *> m_outerloop;
+    // several inner loops, each has some curves
+    std::vector<ON_SimpleArray<ON_Curve *> > m_innerloop;
+    const ON_BrepFace *m_face;
+    enum {
+	UNKNOWN = -1,
+	NOT_BELONG = 0,
+	BELONG = 1
+    } m_belong_to_final;
+    bool m_rev;
+
+    // Default constructor
+    TrimmedFace()
+    {
+	m_face = NULL;
+	m_belong_to_final = UNKNOWN;
+	m_rev = false;
+    }
+
+    // Destructor
+    ~TrimmedFace()
+    {
+	// Delete the curve segments if it's not belong to the result.
+	if (m_belong_to_final != BELONG) {
+	    for (int i = 0; i < m_outerloop.Count(); i++) {
+		if (m_outerloop[i]) {
+		    delete m_outerloop[i];
+		    m_outerloop[i] = NULL;
+		}
+	    }
+	    for (unsigned int i = 0; i < m_innerloop.size(); i++) {
+		for (int j = 0; j < m_innerloop[i].Count(); j++) {
+		    if (m_innerloop[i][j]) {
+			delete m_innerloop[i][j];
+			m_innerloop[i][j] = NULL;
+		    }
+		}
+	    }
+	}
+    }
+
+    TrimmedFace *Duplicate() const
+    {
+	TrimmedFace *out = new TrimmedFace();
+	out->m_face = m_face;
+	for (int i = 0; i < m_outerloop.Count(); i++) {
+	    if (m_outerloop[i]) {
+		out->m_outerloop.Append(m_outerloop[i]->Duplicate());
+	    }
+	}
+	out->m_innerloop = m_innerloop;
+	for (unsigned int i = 0; i < m_innerloop.size(); i++) {
+	    for (int j = 0; j < m_innerloop[i].Count(); j++) {
+		if (m_innerloop[i][j]) {
+		    out->m_innerloop[i][j] = m_innerloop[i][j]->Duplicate();
+		}
+	    }
+	}
+	return out;
+    }
+};
+
+void
+DebugPlot::SplitFaces(
+    const ON_ClassArray<ON_SimpleArray<TrimmedFace *> > &split_faces)
+{
+    for (int i = 0; i < split_faces.Count(); ++i) {
+	for (int j = 0; j < split_faces[i].Count(); ++j) {
+	    TrimmedFace *face = split_faces[i][j];
+
+	    unsigned char *outerloop_color = unknown_outerloop_color;
+	    unsigned char *innerloop_color = unknown_innerloop_color;
+	    switch (face->m_belong_to_final) {
+		case TrimmedFace::NOT_BELONG:
+		    outerloop_color = rejected_outerloop_color;
+		    innerloop_color = rejected_innerloop_color;
+		    break;
+		case TrimmedFace::BELONG:
+		    outerloop_color = accepted_outerloop_color;
+		    innerloop_color = accepted_innerloop_color;
+		    break;
+		default:
+		    outerloop_color = unknown_outerloop_color;
+		    innerloop_color = unknown_innerloop_color;
+	    }
+
+	    int split_face_count = split_face_outerloop_curves.size();
+	    for (int k = 0; k < face->m_outerloop.Count(); ++k) {
+		std::ostringstream filename;
+		filename << prefix << "_split_face" << split_face_count <<
+		    "_outerloop_curve" << k << ".plot3";
+
+		Plot3DCurveFrom2D(face->m_face->SurfaceOf(),
+			face->m_outerloop[k], filename.str().c_str(),
+			outerloop_color);
+	    }
+	    split_face_outerloop_curves.push_back(face->m_outerloop.Count());
+
+	    int innerloop_count = 0;
+	    for (size_t k = 0; k < face->m_innerloop.size(); ++k) {
+		for (int l = 0; l < face->m_innerloop[k].Count(); ++l) {
+		    std::ostringstream filename;
+		    filename << prefix << "_split_face" << split_face_count <<
+			"_innerloop_curve" << innerloop_count++ << ".plot3";
+
+		    Plot3DCurveFrom2D(face->m_face->SurfaceOf(),
+			    face->m_innerloop[k][l], filename.str().c_str(),
+			    innerloop_color);
+		}
+	    }
+	    split_face_innerloop_curves.push_back(innerloop_count);
+	}
+    }
+}
+
+void append_to_polycurve(ON_Curve *curve, ON_PolyCurve &polycurve);
+
+struct SSICurve {
+    ON_Curve *m_curve;
+
+    SSICurve()
+    {
+	m_curve = NULL;
+    }
+
+    SSICurve(ON_Curve *curve)
+    {
+	m_curve = curve;
+    }
+
+    SSICurve *Duplicate() const
+    {
+	SSICurve *out = new SSICurve();
+	if (out != NULL) {
+	    *out = *this;
+	    out->m_curve = m_curve->Duplicate();
+	}
+	return out;
+    }
+};
+
+struct LinkedCurve {
+private:
+    ON_Curve *m_curve;	// an explicit storage of the whole curve
+public:
+    // The curves contained in this LinkedCurve, including
+    // the information needed by the connectivity graph
+    ON_SimpleArray<SSICurve> m_ssi_curves;
+
+    // Default constructor
+    LinkedCurve()
+    {
+	m_curve = NULL;
+    }
+
+    ~LinkedCurve()
+    {
+	if (m_curve) {
+	    delete m_curve;
+	}
+	m_curve = NULL;
+    }
+
+    LinkedCurve &operator= (const LinkedCurve &_lc)
+    {
+	m_curve = _lc.m_curve ? _lc.m_curve->Duplicate() : NULL;
+	m_ssi_curves = _lc.m_ssi_curves;
+	return *this;
+    }
+
+    ON_3dPoint PointAtStart() const
+    {
+	if (m_ssi_curves.Count()) {
+	    return m_ssi_curves[0].m_curve->PointAtStart();
+	} else {
+	    return ON_3dPoint::UnsetPoint;
+	}
+    }
+
+    ON_3dPoint PointAtEnd() const
+    {
+	if (m_ssi_curves.Count()) {
+	    return m_ssi_curves.Last()->m_curve->PointAtEnd();
+	} else {
+	    return ON_3dPoint::UnsetPoint;
+	}
+    }
+
+    bool IsClosed() const
+    {
+	if (m_ssi_curves.Count() == 0) {
+	    return false;
+	}
+	return PointAtStart().DistanceTo(PointAtEnd()) < ON_ZERO_TOLERANCE;
+    }
+
+    bool IsValid() const
+    {
+	// Check whether the curve has "gaps".
+	for (int i = 1; i < m_ssi_curves.Count(); i++) {
+	    if (m_ssi_curves[i].m_curve->PointAtStart().DistanceTo(m_ssi_curves[i - 1].m_curve->PointAtEnd()) >= ON_ZERO_TOLERANCE) {
+		bu_log("The LinkedCurve is not valid.\n");
+		return false;
+	    }
+	}
+	return true;
+    }
+
+    bool Reverse()
+    {
+	ON_SimpleArray<SSICurve> new_array;
+	for (int i = m_ssi_curves.Count() - 1; i >= 0; i--) {
+	    if (!m_ssi_curves[i].m_curve->Reverse()) {
+		return false;
+	    }
+	    new_array.Append(m_ssi_curves[i]);
+	}
+	m_ssi_curves = new_array;
+	return true;
+    }
+
+    void Append(const LinkedCurve &lc)
+    {
+	m_ssi_curves.Append(lc.m_ssi_curves.Count(), lc.m_ssi_curves.Array());
+    }
+
+    void Append(const SSICurve &sc)
+    {
+	m_ssi_curves.Append(sc);
+    }
+
+    void AppendCurvesToArray(ON_SimpleArray<ON_Curve *> &arr) const
+    {
+	for (int i = 0; i < m_ssi_curves.Count(); i++) {
+	    arr.Append(m_ssi_curves[i].m_curve->Duplicate());
+	}
+    }
+
+    const ON_Curve *Curve()
+    {
+	if (m_curve != NULL) {
+	    return m_curve;
+	}
+	if (m_ssi_curves.Count() == 0 || !IsValid()) {
+	    return NULL;
+	}
+	ON_PolyCurve *polycurve = new ON_PolyCurve;
+	for (int i = 0; i < m_ssi_curves.Count(); i++) {
+	    append_to_polycurve(m_ssi_curves[i].m_curve->Duplicate(), *polycurve);
+	}
+	m_curve = polycurve;
+	return m_curve;
+    }
+
+    const ON_3dPoint PointAt(double t)
+    {
+	const ON_Curve *c = Curve();
+	if (c == NULL) {
+	    return ON_3dPoint::UnsetPoint;
+	}
+	return c->PointAt(t);
+    }
+
+    const ON_Interval Domain()
+    {
+	const ON_Curve *c = Curve();
+	if (c == NULL) {
+	    return ON_Interval::EmptyInterval;
+	}
+	return c->Domain();
+    }
+
+    ON_Curve *SubCurve(double t1, double t2)
+    {
+	const ON_Curve *c = Curve();
+	if (c == NULL) {
+	    return NULL;
+	}
+	try {
+	    return sub_curve(c, t1, t2);
+	} catch (InvalidInterval &e) {
+	    bu_log("%s", e.what());
+	    return NULL;
+	}
+    }
+};
+
+void
+DebugPlot::LinkedCurves(
+    const ON_Surface *surf,
+    ON_ClassArray<LinkedCurve> &linked_curves)
+{
+    for (int i = 0; i < linked_curves.Count(); ++i) {
+	const ON_Curve *linked_curve = linked_curves[i].Curve();
+
+	std::ostringstream filename;
+	filename << prefix << "_linked_curve" << linked_curve_count++ << ".plot3";
+
+	Plot3DCurveFrom2D(surf, linked_curve, filename.str().c_str(), transverse_color);
+    }
+}
+
+// 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
@@ -0,0 +1,141 @@
+/*                    D E B U G _ P L O T . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2014 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file debug_plot.h
+ *
+ * Class used to write out unix plot files representing brep geometry
+ * and intersection events from various stages of brep boolean
+ * evaluation.
+ *
+ * Also outputs a parsable log file with implicit information about
+ * the output plot files.
+ */
+
+#ifndef DEBUG_PLOT_H
+#define DEBUG_PLOT_H
+#include "bu/list.h"
+#include "raytrace.h"
+
+struct TrimmedFace;
+struct LinkedCurve;
+
+class DebugPlot {
+public:
+    // prefix all dplot output files with 'basename'
+    DebugPlot(const char *basename);
+
+    ~DebugPlot();
+
+    // write plot files for all surfaces of two breps
+    void Surfaces(const ON_Brep *brep1, const ON_Brep *brep2);
+
+    // record a surface surface intersection
+    void SSX(
+	const ON_ClassArray<ON_SSX_EVENT> &events,
+	const ON_Brep *brep1, int brep1_surf,
+	const ON_Brep *brep2, int brep2_surf);
+
+    // record surface surface isocurve intersections
+    void IsoCSX(
+	const ON_SimpleArray<ON_X_EVENT> &events,
+	const ON_Curve *isocurve,
+	bool is_brep1_iso);
+
+    // record clipped surface-surface intersection curves
+    void ClippedFaceCurves(
+	const ON_Surface *surf1,
+	const ON_Surface *surf2,
+	const ON_SimpleArray<ON_Curve *> &face1_curves,
+	const ON_SimpleArray<ON_Curve *> &face2_curves);
+
+    // write out the log file that the dplot command references to
+    // navigate the generated plot files
+    void WriteLog();
+
+    // get the number of intersecting surface pairs recorded
+    int SurfacePairs();
+
+    // get the number of isocurves that intersect either surface in
+    // given surface pair
+    int IntersectingIsocurves(int ssx_idx);
+
+    // get the number of linked curves recorded
+    int LinkedCurves();
+
+    void SplitFaces(
+	const ON_ClassArray<ON_SimpleArray<TrimmedFace *> > &split_faces);
+
+    void LinkedCurves(
+	const ON_Surface *surf,
+	ON_ClassArray<LinkedCurve> &linked_curves);
+
+    void
+    Plot3DCurveFrom2D(
+	const ON_Surface *surf,
+	const ON_Curve *crv,
+	const char *filename,
+	unsigned char *color,
+	bool decorate = false
+	);
+
+    void
+    Plot3DCurve(
+	const ON_Curve *crv,
+	const char *filename,
+	unsigned char *color,
+	struct bu_list *vlist = NULL);
+
+private:
+    struct bu_list vlist_free_list;
+    std::string prefix;
+    bool have_surfaces;
+    int brep1_surf_count;
+    int brep2_surf_count;
+    int linked_curve_count;
+    std::vector< std::pair<int, int> > intersecting_surfaces; // ssx surface index pairs
+    std::vector<int> ssx_events; // num final events of each ssx
+    std::vector< std::vector<int> > ssx_isocsx_events; // num events for each isocsx of each ssx
+    std::vector<int> ssx_isocsx_brep1_curves; // num ssx isocsx events using brep1 isocurves
+    std::vector< std::pair<int, int> > ssx_clipped_curves; // num clipped intersection curves
+    std::vector<int> split_face_outerloop_curves;
+    std::vector<int> split_face_innerloop_curves;
+
+    void PlotSurface(
+	const ON_Surface &surf,
+	const char *filename,
+	unsigned char *color);
+
+    void
+    PlotBoundaryIsocurves(
+	struct bu_list *vlist,
+	const ON_Surface &surf,
+	int knot_dir);
+};
+
+#endif
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -21,6 +21,8 @@
 *
 * Implementation of intersection routines openNURBS left out.
 *
+ * Additional documentation can be found in the "NURBS Boolean Evaluation
+ * Development Guide" docbook article (bool_eval_development.html).
 */

 #include "common.h"
@@ -38,42 +40,15 @@
 #include "brep/defines.h"
 #include "brep/intersect.h"
 #include "brep/util.h"
+#include "debug_plot.h"
 #include "brep_except.h"
+#include "brep_defines.h"

+extern DebugPlot *dplot;

 // Whether to output the debug messages about b-rep intersections.
 static int DEBUG_BREP_INTERSECT = 0;

-// The maximal depth for subdivision - trade-off between accuracy and
-// performance.
-#define NR_MAX_DEPTH 8
-#define MAX_PCI_DEPTH NR_MAX_DEPTH
-#define MAX_PSI_DEPTH NR_MAX_DEPTH
-#define MAX_CCI_DEPTH NR_MAX_DEPTH
-#define MAX_CSI_DEPTH NR_MAX_DEPTH
-#define MAX_SSI_DEPTH NR_MAX_DEPTH
-
-
-// We make the default tolerance for PSI the same as that of curve and
-// surface intersections defined by openNURBS (see opennurbs_curve.h
-// and opennurbs_surface.h).
-#define NR_DEFAULT_TOLERANCE 0.001
-#define PCI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
-#define PSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
-#define CCI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
-#define CSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
-#define SSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
-
-// Used to prevent an infinite loop in the unlikely event that we
-// can't provide a good starting point for the Newton-Raphson
-// Iteration.
-#define NR_MAX_ITERATIONS 100
-#define PCI_MAX_ITERATIONS NR_MAX_ITERATIONS
-#define PSI_MAX_ITERATIONS NR_MAX_ITERATIONS
-#define CCI_MAX_ITERATIONS NR_MAX_ITERATIONS
-#define CSI_MAX_ITERATIONS NR_MAX_ITERATIONS
-#define SSI_MAX_ITERATIONS NR_MAX_ITERATIONS
-
 class XEventProxy {
 public:
    XEventProxy(ON_X_EVENT::TYPE type);
@@ -2556,10 +2531,29 @@ newton_ssi(double &uA, double &vA, double &uB, double &vB, const ON_Surface *sur
    return (pointA.DistanceTo(pointB) < isect_tol) && !std::isnan(uA) && !std::isnan(vA) && !std::isnan(uB) & !std::isnan(vB);
 }

+// if curve end is greater than tol distance from pt, append a linear
+// segment to the curve so it extends to the pt
+HIDDEN void
+extend_curve_end_to_pt(ON_Curve *&curve, ON_3dPoint pt, double tol)
+{
+    ON_NurbsCurve *nc = curve->NurbsCurve();
+    if (nc->PointAtEnd().DistanceTo(pt) > tol) {
+	ON_LineCurve line_bridge(nc->PointAtEnd(), pt);
+
+	ON_NurbsCurve bridge;
+	if (line_bridge.GetNurbForm(bridge)) {
+	    nc->Append(bridge);
+	}
+	delete curve;
+	curve = nc;
+    }
+}

 HIDDEN ON_Curve *
 link_curves(ON_Curve *&c1, ON_Curve *&c2)
 {
+    extend_curve_end_to_pt(c1, c2->PointAtStart(), ON_ZERO_TOLERANCE);
+
    ON_NurbsCurve *nc1 = c1->NurbsCurve();
    ON_NurbsCurve *nc2 = c2->NurbsCurve();
    if (nc1 && nc2) {
@@ -2578,6 +2572,16 @@ link_curves(ON_Curve *&c1, ON_Curve *&c2)
    return NULL;
 }

+// if curve start and end are within tolerance, append a linear
+// segment to the curve to close it completely so it passes
+// IsClosed() tests
+HIDDEN void
+fill_gap_if_closed(ON_Curve *&curve, double tol)
+{
+    if (curve->PointAtStart().DistanceTo(curve->PointAtEnd()) <= tol) {
+	extend_curve_end_to_pt(curve, curve->PointAtStart(), ON_ZERO_TOLERANCE);
+    }
+}

 struct OverlapSegment {
    ON_Curve *m_curve3d, *m_curveA, *m_curveB;
@@ -2969,17 +2973,10 @@ is_curve_on_overlap_boundary(
    bool ret = false;

    if (iso.overlap_t[0] < iso.overlap_t[1]) {
-	bool at_first_knot = iso.src.knot.IsFirst();
-	bool at_last_knot = iso.src.knot.IsLast();
-	bool closed_at_iso = surf1->IsClosed(1 - iso.src.knot.dir) != 0;
+	int location = isocurve_surface_overlap_location(iso, surf1, surf2,
+		surf2_tree, isect_tol, isect_tol1);

-	ret = true;
-	if (closed_at_iso || (!at_first_knot && !at_last_knot)) {
-	    int location = isocurve_surface_overlap_location(iso, surf1,
-							     surf2, surf2_tree, isect_tol, isect_tol1);
-
-	    ret = (location == ON_OVERLAP_BOUNDARY);
-	}
+	ret = (location == ON_OVERLAP_BOUNDARY);
    }
    return ret;
 }
@@ -3513,6 +3510,7 @@ set_ssx_event_from_curves(
 HIDDEN void
 find_overlap_boundary_curves(
    ON_SimpleArray<OverlapSegment *> &overlaps,
+    ON_ClassArray<ON_SSX_EVENT> &csx_events,
    ON_3dPointArray &isocurve_3d,
    ON_2dPointArray &isocurveA_2d,
    ON_2dPointArray &isocurveB_2d,
@@ -3554,6 +3552,9 @@ find_overlap_boundary_curves(
 	    ON_Intersect(surf1_isocurve, surf2, events, isect_tol,
 			 overlap_tol, 0, 0, 0, &overlap2d);

+	    dplot->IsoCSX(events, surf1_isocurve, is_surfA_iso);
+	    dplot->WriteLog();
+
 	    append_csx_event_points(isocurve_3d, isocurve2_2d, events);

 	    for (int k = 0; k < events.Count(); k++) {
@@ -3575,6 +3576,36 @@ find_overlap_boundary_curves(
 		    if (curve_on_overlap_boundary) {
 			append_overlap_segments(overlaps, iso, overlap2d[k],
 						surf1);
+		    } else {
+			// the intersection isn't part of a
+			// surface-surface overlap, but we want it all
+			// the same
+			ON_SimpleArray<OverlapSegment *> tmp_overlaps;
+
+			append_overlap_segments(tmp_overlaps, iso, overlap2d[k], surf1);
+
+			for (int l = 0; l < tmp_overlaps.Count(); ++l) {
+			    ON_SSX_EVENT ssx_event;
+
+			    int ret = set_ssx_event_from_curves(ssx_event,
+				    tmp_overlaps[l]->m_curve3d,
+				    tmp_overlaps[l]->m_curveA,
+				    tmp_overlaps[l]->m_curveB, surfA, surfB);
+			    if (ret != 0) {
+				bu_log("warning: reverse failed.");
+			    }
+			    csx_events.Append(ssx_event);
+
+			    // set the curves to NULL so they aren't
+			    // deleted by destructors
+			    tmp_overlaps[l]->m_curve3d = NULL;
+			    tmp_overlaps[l]->m_curveA = NULL;
+			    tmp_overlaps[l]->m_curveB = NULL;
+
+			    ssx_event.m_curve3d = NULL;
+			    ssx_event.m_curveA = NULL;
+			    ssx_event.m_curveB = NULL;
+			}
 		    }
 		}
 	    }
@@ -3735,12 +3766,25 @@ ON_Intersect(const ON_Surface *surfA,
    ON_2dPointArray curve_uvA, curve_uvB, tmp_curve_uvA, tmp_curve_uvB;

    ON_SimpleArray<OverlapSegment *> overlaps;
-    find_overlap_boundary_curves(overlaps, tmp_curvept, tmp_curve_uvA,
+    ON_ClassArray<ON_SSX_EVENT> csx_events;
+
+    find_overlap_boundary_curves(overlaps, csx_events, tmp_curvept, tmp_curve_uvA,
 				 tmp_curve_uvB, surfA, surfB, treeA, treeB, isect_tol, isect_tolA,
 				 isect_tolB, overlap_tol);
+
    split_overlaps_at_intersections(overlaps, surfA, surfB, treeA, treeB,
 				    isect_tol, isect_tolA, isect_tolB);

+    // add csx_events
+    for (int i = 0; i < csx_events.Count(); ++i) {
+	x.Append(csx_events[i]);
+
+	// set the curves to NULL so they aren't deleted by destructor
+	csx_events[i].m_curve3d = NULL;
+	csx_events[i].m_curveA = NULL;
+	csx_events[i].m_curveB = NULL;
+    }
+
    // find the neighbors for every overlap segment
    ON_SimpleArray<bool> start_linked(overlaps.Count()), end_linked(overlaps.Count());
    for (int i = 0; i < overlaps.Count(); i++) {
@@ -3875,6 +3919,9 @@ ON_Intersect(const ON_Surface *surfA,
 		delete overlaps[j];
 		overlaps[j] = NULL;
 	    }
+	    fill_gap_if_closed(overlaps[i]->m_curve3d, isect_tol);
+	    fill_gap_if_closed(overlaps[i]->m_curveA, isect_tolA);
+	    fill_gap_if_closed(overlaps[i]->m_curveB, isect_tolB);
 	}
    }

@@ -14,14 +14,8 @@ set(distcheck_files
  vsp.stp
  )

-BRLCAD_ADDEXEC(brep_cdt_mesh brep_cdt_mesh.cpp "libbrep;libbg;libbn;libbu" TEST)
-
-BRLCAD_ADDEXEC(test_point_intersect test_point_intersect.cpp "libbrep" NO_STRICT NO_INSTALL)
-
-BRLCAD_ADDEXEC(test_curve_intersect test_curve_intersect.cpp "libbrep" NO_STRICT NO_INSTALL)
-
-CMAKEFILES(${distcheck_files})
-CMAKEFILES(CMakeLists.txt)
+BRLCAD_ADDEXEC(test_brep_ppx ppx.cpp "libbrep"  NO_INSTALL)
+add_test(NAME brep_ppx COMMAND test_brep_ppx)

 # Local Variables:
 # tab-width: 8
@@ -0,0 +1,230 @@
+/*                         P P X . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2014 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file ppx.cpp
+ *
+ * Point-point intersection tests.
+ *
+ */
+#include "common.h"
+#include "bu.h"
+#include "brep.h"
+#include <limits>
+#include <sstream>
+#include <string>
+#include "../brep_defines.h"
+
+// Here we determine an upper limit on the magnitude of coordinate
+// values for test points.
+//
+// We need to be able to calculate point to point distance as:
+//     sqrt((x2 - x1)^2 + (y2 - y1)^2 + (z2 - z1)^2)
+//
+// To prevent overflow, we must restrict the coordinate values such
+// that:
+//     (x2 - x1)^2 + (y2 - y1)^2 + (z2 - z1)^2 <= DBL_MAX
+//
+// If we want to choose an upper limit M on the magnitude of a
+// coordinate value to enforce this constraint, the restriction
+// simplifies to:
+//     coord1, coord2 in [-M, M] s.t.
+//     (coord2 - coord1)^2 <= DBL_MAX / 3.0
+//  => (coord2 - coord1) <= sqrt(DBL_MAX / 3.0)
+//
+// Knowing that the maximum value of (coord2 - coord1) is 2M, we can
+// solve for M:
+//     2M <= sqrt(DBL_MAX / 3.0)
+//  => M <= sqrt(DBL_MAX / 3.0) / 2.0
+//
+// We choose to use a value slightly less than the true DBL_MAX for
+// the calculation to avoid the possibility that floating point error
+// causes [3.0 * sqrt(DBL_MAX / 3.0)^2 <= DBL_MAX] to be false.
+static const double NEAR_DBL_MAX =
+    std::numeric_limits<double>::max() -
+    std::numeric_limits<double>::round_error();
+
+static const double COORD_MAG_MAX = sqrt(NEAR_DBL_MAX / 3.0) / 2.0;
+
+class PPX_Input {
+public:
+    PPX_Input()
+	: tol(0.0)
+    {
+    }
+    ON_3dPoint ptA;
+    ON_3dPoint ptB;
+    ON_ClassArray<ON_PX_EVENT> events;
+    bool default_tol;
+    double tol;
+};
+
+class PPX_Output {
+public:
+    bool status;
+    ON_ClassArray<ON_PX_EVENT> events;
+
+    bool Equals(std::string &log, const PPX_Output &other);
+};
+
+static const double EVENT_EQUAL_TOL =
+    std::numeric_limits<double>::round_error();
+
+static bool
+point_events_equal(ON_PX_EVENT a, ON_PX_EVENT b)
+{
+    if (a.m_type != b.m_type) {
+	return false;
+    }
+    if (a.m_A.DistanceTo(b.m_A) >= EVENT_EQUAL_TOL) {
+	return false;
+    }
+    if (a.m_B.DistanceTo(b.m_B) >= EVENT_EQUAL_TOL) {
+	return false;
+    }
+    if (a.m_b.DistanceTo(b.m_b) >= EVENT_EQUAL_TOL) {
+	return false;
+    }
+    if (a.m_Mid.DistanceTo(b.m_Mid) >= EVENT_EQUAL_TOL) {
+	return false;
+    }
+    if (!ON_NearZero(a.m_radius - b.m_radius, EVENT_EQUAL_TOL)) {
+	return false;
+    }
+    return true;
+}
+
+bool
+PPX_Output::Equals(std::string &log, const PPX_Output &other)
+{
+    bool ret = true;
+    std::stringstream out;
+
+    if (status != other.status) {
+	out << "return status: " << std::boolalpha << status << " vs " <<
+	    other.status << "\n";
+	ret = false;
+    }
+
+    if (events.Count() != other.events.Count()) {
+	out << "event count: " << events.Count() << " vs " <<
+	    other.events.Count() << "\n";
+	ret = false;
+    } else {
+	for (int i = 0; i < events.Count(); ++i) {
+	    if (!point_events_equal(other.events[i], events[i])) {
+		out << "event arrays don't match\n";
+		ret = false;
+		break;
+	    }
+	}
+    }
+
+    log = out.str();
+
+    return ret;
+}
+
+void
+test_intersection(PPX_Input in, PPX_Output expected_out)
+{
+    PPX_Output out;
+
+    out.events = in.events;
+
+    if (in.default_tol) {
+	out.status = ON_Intersect(in.ptA, in.ptB, out.events);
+    } else {
+	out.status = ON_Intersect(in.ptA, in.ptB, out.events, in.tol);
+    }
+
+    std::string err_msg;
+    if (!expected_out.Equals(err_msg, out)) {
+	bu_exit(1, "Unexpected intersection result. Expected vs actual:\n%s",
+		err_msg.c_str());
+    }
+}
+
+static ON_PX_EVENT
+ppx_point_event(ON_3dPoint ptA, ON_3dPoint ptB)
+{
+    ON_PX_EVENT event;
+    event.m_type = ON_PX_EVENT::ppx_point;
+    event.m_A = ptA;
+    event.m_B = ptB;
+    event.m_b = ON_2dPoint(0.0, 0.0);
+    event.m_Mid = (ptA + ptB) * 0.5;
+    event.m_radius = ptA.DistanceTo(ptB) * 0.5;
+    return event;
+}
+
+static void
+test_equal_points(ON_3dPoint pt)
+{
+    ON_ClassArray<ON_PX_EVENT> events;
+    events.Append(ppx_point_event(pt, pt));
+
+    PPX_Output expected_output;
+    expected_output.status = true;
+    expected_output.events = events;
+
+    PPX_Input input;
+    input.ptA = input.ptB = pt;
+    input.default_tol = true;
+
+    test_intersection(input, expected_output);
+}
+
+static void
+do_equal_point_tests(void)
+{
+    std::vector<double> coord_vals;
+    coord_vals.push_back(-COORD_MAG_MAX);
+    coord_vals.push_back(-INTERSECTION_TOL);
+    coord_vals.push_back(0.0);
+    coord_vals.push_back(INTERSECTION_TOL);
+    coord_vals.push_back(COORD_MAG_MAX);
+
+    for (size_t i = 0; i < coord_vals.size(); ++i) {
+	for (size_t j = 0; j < coord_vals.size(); ++j) {
+	    for (size_t k = 0; k < coord_vals.size(); ++k) {
+		ON_3dPoint test_pt(coord_vals[i], coord_vals[j],
+			coord_vals[k]);
+
+		test_equal_points(test_pt);
+	    }
+	}
+    }
+}
+
+int
+main(void)
+{
+    do_equal_point_tests();
+
+    return 0;
+}
+
+// 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
@@ -1,12 +1,32 @@
+
+set(DPLOT_READER dplot)
+
+PERPLEX_TARGET(
+  ${DPLOT_READER}_scanner
+  ${DPLOT_READER}_scanner.perplex
+  OUT_SRC_FILE ${DPLOT_READER}_scanner.c
+  OUT_HDR_FILE ${DPLOT_READER}_scanner.h)
+LEMON_TARGET(
+  ${DPLOT_READER}_parser
+  ${DPLOT_READER}_parser.lemon
+  OUT_SRC_FILE ${DPLOT_READER}_parser.c
+  OUT_HDR_FILE ${DPLOT_READER}_parser.h)
+ADD_PERPLEX_LEMON_DEPENDENCY(${DPLOT_READER}_scanner ${DPLOT_READER}_parser)
+
 include_directories(
  ${CMAKE_CURRENT_SOURCE_DIR}
  ${BRLCAD_BINARY_DIR}/include
  ${BRLCAD_SOURCE_DIR}/include
  ${BU_INCLUDE_DIRS}
  ${GED_INCLUDE_DIRS}
+  ${PERPLEX_${DPLOT_READER}_scanner_INCLUDE_DIR}
+  ${LEMON_${DPLOT_READER}_parser_INCLUDE_DIR}
  )

 set(BREP_SRCS
+  ${PERPLEX_${DPLOT_READER}_scanner_SRC}
+  ${LEMON_${DPLOT_READER}_parser_SRC}
+  dplot.c
  brep.cpp
  conversion.cpp
  csg.cpp
@@ -1480,9 +1480,11 @@ ged_brep_core(struct ged *gedp, int argc, const char *argv[])
 extern "C" {
    struct ged_cmd_impl brep_cmd_impl = { "brep", ged_brep_core, GED_CMD_DEFAULT };
    const struct ged_cmd brep_cmd = { &brep_cmd_impl };
-    const struct ged_cmd *brep_cmds[] = { &brep_cmd,  NULL };
+    struct ged_cmd_impl dplot_cmd_impl = { "dplot", ged_dplot_core, GED_CMD_DEFAULT };
+    const struct ged_cmd dplot_cmd = { &dplot_cmd_impl };
+    const struct ged_cmd *brep_cmds[] = { &brep_cmd, &dplot_cmd, NULL };

-    static const struct ged_plugin pinfo = { GED_API,  brep_cmds, 1 };
+    static const struct ged_plugin pinfo = { GED_API,  brep_cmds, 2 };

    COMPILER_DLLEXPORT const struct ged_plugin *ged_plugin_info()
    {
@@ -0,0 +1,799 @@
+/*                         D P L O T . C
+ * BRL-CAD
+ *
+ * Copyright (c) 2014-2020 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file dplot.c
+ *
+ * Brief description
+ *
+ */
+
+#include "common.h"
+
+#include <stdlib.h>
+#include <string.h>
+
+#include "bu/color.h"
+#include "bu/opt.h"
+#include "raytrace.h"
+#include "rt/geom.h"
+#include "wdb.h"
+
+#include "../ged_private.h"
+#include "./dplot_reader.h"
+
+enum {
+    DPLOT_INITIAL,
+    DPLOT_SSX_FIRST,
+    DPLOT_SSX,
+    DPLOT_SSX_EVENTS,
+    DPLOT_ISOCSX_FIRST,
+    DPLOT_ISOCSX,
+    DPLOT_ISOCSX_EVENTS,
+    DPLOT_FACE_CURVES,
+    DPLOT_LINKED_CURVES,
+    DPLOT_SPLIT_FACES
+};
+
+struct dplot_info {
+    struct dplot_data fdata;
+    int mode;
+    struct ged *gedp;
+    FILE *logfile;
+    char *prefix;
+    int ssx_idx;
+    int isocsx_idx;
+    int brep1_surf_idx;
+    int brep2_surf_idx;
+    int brep1_surf_count;
+    int brep2_surf_count;
+    int event_idx;
+    int event_count;
+    int brep1_isocsx_count;
+    int isocsx_count;
+};
+
+#define CLEANUP \
+    fclose(info.logfile); \
+    bu_free(info.prefix, "prefix"); \
+    info.prefix = NULL; \
+    if (info.fdata.ssx) bu_free(info.fdata.ssx, "ssx array");
+
+#define RETURN_MORE \
+    CLEANUP \
+    return GED_MORE;
+
+#define RETURN_ERROR \
+    CLEANUP \
+    info.mode = DPLOT_INITIAL; \
+    return GED_ERROR;
+
+HIDDEN int
+dplot_overlay(
+	struct ged *gedp,
+	const char *prefix,
+	const char *infix,
+	int idx,
+	const char *name)
+{
+    const char *cmd_av[] = {"overlay", "[filename]", "1.0", "[name]"};
+    int ret, cmd_ac = sizeof(cmd_av) / sizeof(char *);
+    struct bu_vls overlay_name = BU_VLS_INIT_ZERO;
+
+    bu_vls_printf(&overlay_name, "%s%s%d.plot3", prefix, infix, idx);
+    cmd_av[1] = cmd_av[3] = bu_vls_cstr(&overlay_name);
+    if (name) {
+	cmd_av[3] = name;
+    }
+    ret = ged_overlay(gedp, cmd_ac, cmd_av);
+    bu_vls_free(&overlay_name);
+
+    if (ret != GED_OK) {
+	bu_vls_printf(gedp->ged_result_str, "error overlaying plot\n");
+	return GED_ERROR;
+    }
+    return GED_OK;
+}
+
+HIDDEN int
+dplot_erase_overlay(
+	struct dplot_info *info,
+	const char *name)
+{
+    const int NUM_EMPTY_PLOTS = 13;
+    int i;
+
+    /* We can't actually erase the old plot without its real name,
+     * which is unknown. Instead, we'll write a plot with the same
+     * base name and color, which will overwrite the old one. We
+     * don't actually know the color either, so we resort to writing
+     * an empty plot with the given name using every color we created
+     * plots with.
+     */
+    for (i = 0; i < NUM_EMPTY_PLOTS; ++i) {
+	int ret = dplot_overlay(info->gedp, info->prefix, "_empty", i, name);
+	if (ret != GED_OK) {
+	    return ret;
+	}
+    }
+    return GED_OK;
+}
+
+HIDDEN int
+dplot_ssx(
+	struct dplot_info *info)
+{
+    int i, ret;
+
+    /* draw surfaces, skipping intersecting surfaces if in SSI mode */
+    /* TODO: need to name these overlays so I can selectively erase them */
+    for (i = 0; i < info->brep1_surf_count; ++i) {
+	struct bu_vls name = BU_VLS_INIT_ZERO;
+	bu_vls_printf(&name, "%s_brep1_surface%d", info->prefix, i);
+	dplot_erase_overlay(info, bu_vls_cstr(&name));
+	bu_vls_free(&name);
+    }
+    for (i = 0; i < info->brep2_surf_count; ++i) {
+	struct bu_vls name = BU_VLS_INIT_ZERO;
+	bu_vls_printf(&name, "%s_brep2_surface%d", info->prefix, i);
+	dplot_erase_overlay(info, bu_vls_cstr(&name));
+	bu_vls_free(&name);
+    }
+    if (info->mode == DPLOT_SSX_FIRST || info->mode == DPLOT_SSX || info->mode == DPLOT_SSX_EVENTS) {
+	for (i = 0; i < info->brep1_surf_count; ++i) {
+	    if (info->mode != DPLOT_SSX_FIRST && i == info->brep1_surf_idx) {
+		continue;
+	    }
+	    ret = dplot_overlay(info->gedp, info->prefix, "_brep1_surface", i, NULL);
+	    if (ret != GED_OK) {
+		return GED_ERROR;
+	    }
+	}
+	for (i = 0; i < info->brep2_surf_count; ++i) {
+	    if (info->mode != DPLOT_SSX_FIRST && i == info->brep2_surf_idx) {
+		continue;
+	    }
+	    ret = dplot_overlay(info->gedp, info->prefix, "_brep2_surface", i, NULL);
+	    if (ret != GED_OK) {
+		return GED_ERROR;
+	    }
+	}
+    }
+
+    /* draw highlighted surface-surface intersection pair */
+    if (info->mode == DPLOT_SSX || info->mode == DPLOT_SSX_EVENTS) {
+	ret = dplot_overlay(info->gedp, info->prefix, "_highlight_brep1_surface",
+		info->brep1_surf_idx, "dplot_ssx1");
+	if (ret != GED_OK) {
+	    return GED_ERROR;
+	}
+	ret = dplot_overlay(info->gedp, info->prefix, "_highlight_brep2_surface",
+		info->brep2_surf_idx, "dplot_ssx2");
+	if (ret != GED_OK) {
+	    return GED_ERROR;
+	}
+	if (info->mode == DPLOT_SSX) {
+	    /* advance past the completed pair */
+	    ++info->ssx_idx;
+	}
+    }
+    if (info->mode == DPLOT_SSX_FIRST) {
+	info->mode = DPLOT_SSX;
+    }
+    if (info->mode == DPLOT_SSX && info->ssx_idx < info->fdata.ssx_count) {
+	bu_vls_printf(info->gedp->ged_result_str, "Press [Enter] to show surface-"
+		"surface intersection %d", info->ssx_idx);
+	return GED_MORE;
+    }
+    return GED_OK;
+}
+
+void
+dplot_print_event_legend(struct dplot_info *info)
+{
+    bu_vls_printf(info->gedp->ged_result_str, "yellow = transverse\n");
+    bu_vls_printf(info->gedp->ged_result_str, "white  = tangent\n");
+    bu_vls_printf(info->gedp->ged_result_str, "green  = overlap\n");
+}
+
+HIDDEN int
+dplot_ssx_events(
+	struct dplot_info *info)
+{
+    int ret;
+
+    /* erase old event plots */
+    ret = dplot_erase_overlay(info, "curr_event");
+    if (ret != GED_OK) {
+	return ret;
+    }
+
+    if (info->mode != DPLOT_SSX_EVENTS) {
+	return GED_OK;
+    }
+
+    if (info->event_count > 0) {
+	/* convert event ssx_idx to string */
+	struct bu_vls infix = BU_VLS_INIT_ZERO;
+	bu_vls_printf(&infix, "_ssx%d_event", info->ssx_idx);
+
+	/* plot overlay */
+	ret = dplot_overlay(info->gedp, info->prefix, bu_vls_cstr(&infix),
+		info->event_idx, "curr_event");
+	bu_vls_free(&infix);
+
+	if (ret != GED_OK) {
+	    return ret;
+	}
+	if (info->event_idx == 0) {
+	    dplot_print_event_legend(info);
+	}
+    }
+    /* advance to next event, or return to initial state */
+    if (++info->event_idx < info->event_count) {
+	bu_vls_printf(info->gedp->ged_result_str, "Press [Enter] to show next event\n");
+	return GED_MORE;
+    }
+    info->mode = DPLOT_INITIAL;
+    return GED_OK;
+}
+
+HIDDEN int
+dplot_isocsx(
+	struct dplot_info *info)
+{
+    if (info->mode != DPLOT_ISOCSX &&
+	    info->mode != DPLOT_ISOCSX_FIRST &&
+	    info->mode != DPLOT_ISOCSX_EVENTS)
+    {
+	return GED_OK;
+    }
+
+    if (info->fdata.ssx[info->ssx_idx].isocsx_events == NULL) {
+	bu_vls_printf(info->gedp->ged_result_str, "The isocurves of neither "
+		"surface intersected the opposing surface in surface-surface"
+		" intersection %d.\n", info->ssx_idx);
+	info->mode = DPLOT_INITIAL;
+	return GED_OK;
+    }
+
+    dplot_overlay(info->gedp, info->prefix, "_brep1_surface",
+	    info->brep1_surf_idx, "isocsx_b1");
+    dplot_overlay(info->gedp, info->prefix, "_brep2_surface",
+	    info->brep2_surf_idx, "isocsx_b2");
+
+    if (info->mode == DPLOT_ISOCSX) {
+	struct bu_vls infix = BU_VLS_INIT_ZERO;
+	/* plot surface and the isocurve that intersects it */
+	if (info->isocsx_idx < info->brep1_isocsx_count) {
+	    dplot_overlay(info->gedp, info->prefix, "_highlight_brep2_surface",
+		    info->brep2_surf_idx, "isocsx_b2");
+	} else {
+	    dplot_overlay(info->gedp, info->prefix, "_highlight_brep1_surface",
+		    info->brep1_surf_idx, "isocsx_b1");
+	}
+	bu_vls_printf(&infix, "_highlight_ssx%d_isocurve", info->ssx_idx);
+	dplot_overlay(info->gedp, info->prefix, bu_vls_cstr(&infix), info->isocsx_idx, "isocsx_isocurve");
+	bu_vls_free(&infix);
+    }
+
+    if (info->mode == DPLOT_ISOCSX_FIRST ||
+	    info->mode == DPLOT_ISOCSX)
+    {
+	if (info->mode == DPLOT_ISOCSX_FIRST ||
+		++info->isocsx_idx < info->isocsx_count)
+	{
+	    bu_vls_printf(info->gedp->ged_result_str, "Press [Enter] to show "
+		    "isocurve-surface intersection %d", info->isocsx_idx);
+	    info->mode = DPLOT_ISOCSX;
+	    return GED_MORE;
+	} else {
+	    info->mode = DPLOT_INITIAL;
+	}
+    }
+    return GED_OK;
+}
+
+HIDDEN int
+dplot_isocsx_events(struct dplot_info *info)
+{
+    int ret;
+
+    if (info->mode != DPLOT_ISOCSX_EVENTS) {
+	return GED_OK;
+    }
+    ret = dplot_erase_overlay(info, "curr_event");
+    if (ret != GED_OK) {
+	return ret;
+    }
+    if (info->event_count > 0) {
+	/* convert event ssx_idx to string */
+	struct bu_vls infix = BU_VLS_INIT_ZERO;
+	bu_vls_printf(&infix, "_ssx%d_isocsx%d_event", info->ssx_idx,
+		info->isocsx_idx);
+
+	/* plot overlay */
+	ret = dplot_overlay(info->gedp, info->prefix, bu_vls_cstr(&infix),
+		info->event_idx, "curr_event");
+	bu_vls_free(&infix);
+
+	if (ret != GED_OK) {
+	    bu_vls_printf(info->gedp->ged_result_str,
+		    "error overlaying plot\n");
+	    return ret;
+	}
+	if (info->event_idx == 0) {
+	    dplot_print_event_legend(info);
+	}
+    }
+    /* advance to next event, or return to initial state */
+    if (++info->event_idx < info->event_count) {
+	bu_vls_printf(info->gedp->ged_result_str,
+		"Press [Enter] to show next event\n");
+	return GED_MORE;
+    }
+
+    info->mode = DPLOT_INITIAL;
+    return GED_OK;
+}
+
+HIDDEN int
+dplot_face_curves(struct dplot_info *info)
+{
+    int f1_curves, f2_curves;
+    if (info->mode != DPLOT_FACE_CURVES) {
+	return GED_OK;
+    }
+
+    f1_curves = info->fdata.ssx[info->ssx_idx].face1_clipped_curves;
+    f2_curves = info->fdata.ssx[info->ssx_idx].face2_clipped_curves;
+    info->event_count = f1_curves + f2_curves;
+
+    if (info->event_count == 0) {
+	bu_vls_printf(info->gedp->ged_result_str, "No clipped curves for ssx"
+		" pair %d.\n", info->ssx_idx);
+	return GED_OK;
+    }
+
+    if (info->event_idx < info->event_count) {
+	struct bu_vls prefix;
+
+	dplot_overlay(info->gedp, info->prefix, "_brep1_surface",
+		info->brep1_surf_idx, "face_b1");
+	dplot_overlay(info->gedp, info->prefix, "_brep2_surface",
+		info->brep2_surf_idx, "face_b2");
+
+	BU_VLS_INIT(&prefix);
+	bu_vls_printf(&prefix, "%s_ssx%d", info->prefix, info->ssx_idx);
+	dplot_erase_overlay(info, "clipped_fcurve");
+	if (info->event_idx < f1_curves) {
+	    bu_vls_printf(&prefix, "_brep1face_clipped_curve");
+	    dplot_overlay(info->gedp, bu_vls_cstr(&prefix), "",
+		    info->event_idx, "clipped_fcurve");
+	} else {
+	    bu_vls_printf(&prefix, "_brep2face_clipped_curve");
+	    dplot_overlay(info->gedp, bu_vls_cstr(&prefix), "",
+		    info->event_idx - f1_curves, "clipped_fcurve");
+	}
+	++info->event_idx;
+	if (info->event_idx < info->event_count) {
+	    bu_vls_printf(info->gedp->ged_result_str, "Press [Enter] to show the"
+		    " next curve.");
+	    return GED_MORE;
+	}
+    }
+
+    info->mode = DPLOT_INITIAL;
+    return GED_OK;
+}
+
+HIDDEN int
+dplot_split_faces(
+	struct dplot_info *info)
+{
+    int i, j;
+    struct bu_vls name = BU_VLS_INIT_ZERO;
+    struct bu_vls short_name = BU_VLS_INIT_ZERO;
+    struct split_face split_face;
+
+    if (info->mode != DPLOT_SPLIT_FACES) {
+	return GED_OK;
+    }
+
+    if (info->event_idx >= info->fdata.split_face_count) {
+	for (i = 0; i < info->fdata.split_face_count; ++i) {
+	    split_face = info->fdata.face[i];
+
+	    bu_vls_trunc(&name, 0);
+	    bu_vls_printf(&name, "_split_face%d_outerloop_curve", i);
+	    for (j = 0; j < split_face.outerloop_curves; ++j) {
+		bu_vls_trunc(&short_name, 0);
+		bu_vls_printf(&short_name, "sf%do%d", i, j);
+		dplot_overlay(info->gedp, info->prefix, bu_vls_cstr(&name), j,
+			bu_vls_cstr(&short_name));
+	    }
+
+	    bu_vls_trunc(&name, 0);
+	    bu_vls_printf(&name, "_split_face%d_innerloop_curve", i);
+	    for (j = 0; j < split_face.innerloop_curves; ++j) {
+		bu_vls_trunc(&short_name, 0);
+		bu_vls_printf(&short_name, "sf%di%d", i, j);
+		dplot_overlay(info->gedp, info->prefix, bu_vls_cstr(&name), j,
+			bu_vls_cstr(&short_name));
+	    }
+	}
+    } else {
+	if (info->event_idx > 0) {
+	    /* erase curves of previous split face */
+	    split_face = info->fdata.face[info->event_idx - 1];
+	    for (i = 0; i < split_face.outerloop_curves; ++i) {
+		bu_vls_trunc(&short_name, 0);
+		bu_vls_printf(&short_name, "sfo%d", i);
+		dplot_erase_overlay(info, bu_vls_cstr(&short_name));
+	    }
+	    for (i = 0; i < split_face.innerloop_curves; ++i) {
+		bu_vls_trunc(&short_name, 0);
+		bu_vls_printf(&short_name, "sfi%d", i);
+		dplot_erase_overlay(info, bu_vls_cstr(&short_name));
+	    }
+	}
+
+	split_face = info->fdata.face[info->event_idx];
+	bu_vls_printf(&name, "_split_face%d_outerloop_curve", info->event_idx);
+	for (i = 0; i < info->fdata.face[info->event_idx].outerloop_curves; ++i) {
+	    bu_vls_trunc(&short_name, 0);
+	    bu_vls_printf(&short_name, "sfo%d", i);
+
+	    dplot_erase_overlay(info, bu_vls_cstr(&short_name));
+	    dplot_overlay(info->gedp, info->prefix, bu_vls_cstr(&name), i,
+		    bu_vls_cstr(&short_name));
+	}
+
+	bu_vls_trunc(&name, 0);
+	bu_vls_printf(&name, "_split_face%d_innerloop_curve", info->event_idx);
+	for (i = 0; i < info->fdata.face[info->event_idx].innerloop_curves; ++i) {
+	    bu_vls_trunc(&short_name, 0);
+	    bu_vls_printf(&short_name, "sfi%d", i);
+
+	    dplot_erase_overlay(info, bu_vls_cstr(&short_name));
+	    dplot_overlay(info->gedp, info->prefix, bu_vls_cstr(&name), i,
+		    bu_vls_cstr(&short_name));
+	}
+
+	bu_vls_printf(info->gedp->ged_result_str, "Press [Enter] to show "
+		"split face %d", ++info->event_idx);
+	return GED_MORE;
+    }
+    return GED_OK;
+}
+
+HIDDEN int
+dplot_linked_curves(
+	struct dplot_info *info)
+{
+    int i;
+    if (info->mode != DPLOT_LINKED_CURVES) {
+	return GED_OK;
+    }
+
+    if (info->event_idx >= info->fdata.linked_curve_count) {
+	for (i = 0; i < info->fdata.linked_curve_count; ++i) {
+	    dplot_overlay(info->gedp, info->prefix, "_linked_curve", i, NULL);
+	}
+    } else {
+	dplot_overlay(info->gedp, info->prefix, "_linked_curve",
+		info->event_idx, "linked_curve");
+	bu_vls_printf(info->gedp->ged_result_str, "Press [Enter] to show "
+		"linked curve %d", ++info->event_idx);
+	return GED_MORE;
+    }
+    return GED_OK;
+}
+
+HIDDEN void *
+dplot_malloc(size_t s) {
+    return bu_malloc(s, "dplot_malloc");
+}
+
+HIDDEN void
+dplot_free(void *p) {
+    bu_free(p, "dplot_free");
+}
+
+HIDDEN void
+dplot_load_file_data(struct dplot_info *info)
+{
+    int i, j;
+    struct ssx *curr_ssx;
+    struct isocsx *curr_isocsx;
+    int token_id;
+    perplex_t scanner;
+    void *parser;
+
+    /* initialize scanner and parser */
+    parser = ParseAlloc(dplot_malloc);
+    scanner = perplexFileScanner(info->logfile);
+
+    info->fdata.brep1_surface_count = info->fdata.brep2_surface_count = 0;
+    info->fdata.ssx_count = 0;
+    BU_LIST_INIT(&info->fdata.ssx_list);
+    BU_LIST_INIT(&info->fdata.isocsx_list);
+    perplexSetExtra(scanner, (void *)&info->fdata);
+
+    /* parse */
+    while ((token_id = yylex(scanner)) != YYEOF) {
+	Parse(parser, token_id, info->fdata.token_data, &info->fdata);
+    }
+    Parse(parser, 0, info->fdata.token_data, &info->fdata);
+
+    /* clean up */
+    ParseFree(parser, dplot_free);
+    perplexFree(scanner);
+
+    /* move ssx to dynamic array for easy access */
+    info->fdata.ssx = NULL;
+    if (info->fdata.ssx_count > 0) {
+
+	info->fdata.ssx = (struct ssx *)bu_malloc(
+		sizeof(struct ssx) * info->fdata.ssx_count, "ssx array");
+
+	i = info->fdata.ssx_count - 1;
+	while (BU_LIST_WHILE(curr_ssx, ssx, &info->fdata.ssx_list)) {
+	    BU_LIST_DEQUEUE(&curr_ssx->l);
+
+	    curr_ssx->isocsx_events = NULL;
+	    if (curr_ssx->intersecting_isocurves > 0) {
+		curr_ssx->isocsx_events = (int *)bu_malloc(sizeof(int) *
+			curr_ssx->intersecting_isocurves, "isocsx array");
+
+		j = curr_ssx->intersecting_isocurves - 1;
+		while (BU_LIST_WHILE(curr_isocsx, isocsx, &curr_ssx->isocsx_list)) {
+		    BU_LIST_DEQUEUE(&curr_isocsx->l);
+		    curr_ssx->isocsx_events[j--] = curr_isocsx->events;
+		    BU_PUT(curr_isocsx, struct isocsx);
+		}
+	    }
+	    info->fdata.ssx[i--] = *curr_ssx;
+	    BU_PUT(curr_ssx, struct ssx);
+	}
+    }
+}
+
+int
+ged_dplot_core(struct ged *gedp, int argc, const char *argv[])
+{
+    static struct dplot_info info;
+    int ret;
+    const char *filename, *cmd;
+    char *dot;
+
+    info.gedp = gedp;
+    bu_vls_trunc(gedp->ged_result_str, 0);
+
+    if (argc < 3) {
+	bu_vls_printf(gedp->ged_result_str, "usage: %s logfile cmd\n",
+		argv[0]);
+	bu_vls_printf(gedp->ged_result_str, "  where cmd is one of:\n");
+	bu_vls_printf(gedp->ged_result_str,
+		"      ssx     (show intersecting surface pairs)\n");
+	bu_vls_printf(gedp->ged_result_str,
+		"      ssx N   (show intersections of ssx pair N)\n");
+	bu_vls_printf(gedp->ged_result_str,
+		"   isocsx N   (show intersecting isocurve-surface pairs of ssx pair N)\n");
+	bu_vls_printf(gedp->ged_result_str,
+		"   isocsx N M (show intersections of ssx pair N, isocsx pair M)\n");
+	bu_vls_printf(gedp->ged_result_str,
+		"  fcurves N   (show clipped face curves of ssx pair N)\n");
+	bu_vls_printf(gedp->ged_result_str,
+		"  lcurves     (show linked ssx curves used to split faces)\n");
+	bu_vls_printf(gedp->ged_result_str,
+		"    faces     (show split faces used to construct result)\n");
+	return GED_HELP;
+    }
+    filename = argv[1];
+    cmd = argv[2];
+
+    if (info.mode == DPLOT_INITIAL) {
+	if (BU_STR_EQUAL(cmd, "ssx") && argc == 3) {
+	    info.mode = DPLOT_SSX_FIRST;
+	    info.ssx_idx = 0;
+	} else if (BU_STR_EQUAL(cmd, "ssx") && argc == 4) {
+	    /* parse surface pair index */
+	    const char *idx_str = argv[3];
+	    ret = bu_sscanf(idx_str, "%d", &info.ssx_idx);
+	    if (ret != 1) {
+		bu_vls_printf(gedp->ged_result_str, "%s is not a valid "
+			"surface pair (must be a non-negative integer)\n", idx_str);
+		return GED_ERROR;
+	    }
+	    info.mode = DPLOT_SSX_EVENTS;
+	    info.event_idx = 0;
+	} else if (BU_STR_EQUAL(cmd, "isocsx") && argc == 4) {
+	    const char *idx_str = argv[3];
+	    ret = bu_sscanf(idx_str, "%d", &info.ssx_idx);
+	    if (ret != 1) {
+		bu_vls_printf(gedp->ged_result_str, "%s is not a valid "
+			"surface pair (must be a non-negative integer)\n", idx_str);
+		return GED_ERROR;
+	    }
+	    info.mode = DPLOT_ISOCSX_FIRST;
+	    info.isocsx_idx = 0;
+	} else if (BU_STR_EQUAL(cmd, "isocsx") && argc == 5) {
+	    /* parse surface pair index */
+	    const char *idx_str = argv[3];
+	    ret = bu_sscanf(idx_str, "%d", &info.ssx_idx);
+	    if (ret != 1) {
+		bu_vls_printf(gedp->ged_result_str, "%s is not a valid "
+			"surface pair (must be a non-negative integer)\n", idx_str);
+		return GED_ERROR;
+	    }
+	    idx_str = argv[4];
+	    ret = bu_sscanf(idx_str, "%d", &info.isocsx_idx);
+	    if (ret != 1) {
+		bu_vls_printf(gedp->ged_result_str, "%s is not a valid "
+			"isocurve-surface pair (must be a non-negative integer)\n", idx_str);
+		return GED_ERROR;
+	    }
+	    info.mode = DPLOT_ISOCSX_EVENTS;
+	    info.event_idx = 0;
+	} else if (BU_STR_EQUAL(cmd, "fcurves") && argc == 4) {
+	    const char *idx_str = argv[3];
+	    ret = bu_sscanf(idx_str, "%d", &info.ssx_idx);
+	    if (ret != 1) {
+		bu_vls_printf(gedp->ged_result_str, "%s is not a valid "
+			"surface pair (must be a non-negative integer)\n", idx_str);
+		return GED_ERROR;
+	    }
+	    info.mode = DPLOT_FACE_CURVES;
+	    info.event_idx = 0;
+	} else if (BU_STR_EQUAL(cmd, "lcurves") && argc == 3) {
+	    info.mode = DPLOT_LINKED_CURVES;
+	    info.event_idx = 0;
+	} else if (BU_STR_EQUAL(cmd, "faces") && argc == 3) {
+	    info.mode = DPLOT_SPLIT_FACES;
+	    info.event_idx = 0;
+	} else {
+	    bu_vls_printf(gedp->ged_result_str, "%s is not a recognized "
+		    "command or was given the wrong number of arguments\n",
+		    cmd);
+	    return GED_ERROR;
+	}
+    }
+
+    /* open dplot log file */
+    info.logfile = fopen(filename, "r");
+    if (!info.logfile) {
+	bu_vls_printf(gedp->ged_result_str, "couldn't open log file \"%s\"\n", filename);
+	return GED_ERROR;
+    }
+
+    /* filename before '.' is assumed to be the prefix for all
+     * plot-file names
+     */
+    info.prefix = bu_strdup(filename);
+    dot = strchr(info.prefix, '.');
+    if (dot) {
+	*dot = '\0';
+    }
+
+    dplot_load_file_data(&info);
+
+    if (info.mode == DPLOT_SSX_FIRST	||
+	    info.mode == DPLOT_SSX		||
+	    info.mode == DPLOT_SSX_EVENTS	||
+	    info.mode == DPLOT_ISOCSX_FIRST	||
+	    info.mode == DPLOT_ISOCSX	||
+	    info.mode == DPLOT_ISOCSX_EVENTS)
+    {
+	if (info.fdata.ssx_count == 0) {
+	    bu_vls_printf(info.gedp->ged_result_str, "no surface surface"
+		    "intersections");
+	    RETURN_ERROR;
+	} else if (info.ssx_idx < 0 ||
+		info.ssx_idx > (info.fdata.ssx_count - 1))
+	{
+	    bu_vls_printf(info.gedp->ged_result_str, "no surface pair %d (valid"
+		    " range is [0, %d])\n", info.ssx_idx, info.fdata.ssx_count - 1);
+	    RETURN_ERROR;
+	}
+    }
+    if (info.fdata.ssx_count > 0) {
+	info.brep1_surf_idx = info.fdata.ssx[info.ssx_idx].brep1_surface;
+	info.brep2_surf_idx = info.fdata.ssx[info.ssx_idx].brep2_surface;
+	info.event_count = info.fdata.ssx[info.ssx_idx].final_curve_events;
+	info.brep1_isocsx_count =
+	    info.fdata.ssx[info.ssx_idx].intersecting_brep1_isocurves;
+	info.isocsx_count =
+	    info.fdata.ssx[info.ssx_idx].intersecting_isocurves;
+    }
+    info.brep1_surf_count = info.fdata.brep1_surface_count;
+    info.brep2_surf_count = info.fdata.brep2_surface_count;
+
+    if (info.mode == DPLOT_ISOCSX_EVENTS) {
+	int *isocsx_events = info.fdata.ssx[info.ssx_idx].isocsx_events;
+
+	info.event_count = 0;
+	if (isocsx_events) {
+	    info.event_count = isocsx_events[info.event_idx];
+	}
+    }
+
+    ret = dplot_ssx(&info);
+    if (ret == GED_ERROR) {
+	RETURN_ERROR;
+    } else if (ret == GED_MORE) {
+	RETURN_MORE;
+    }
+
+    ret = dplot_ssx_events(&info);
+    if (ret == GED_ERROR) {
+	RETURN_ERROR;
+    } else if (ret == GED_MORE) {
+	RETURN_MORE;
+    }
+
+    ret = dplot_isocsx(&info);
+    if (ret == GED_ERROR) {
+	RETURN_ERROR;
+    } else if (ret == GED_MORE) {
+	RETURN_MORE;
+    }
+
+    ret = dplot_isocsx_events(&info);
+    if (ret == GED_ERROR) {
+	RETURN_ERROR;
+    } else if (ret == GED_MORE) {
+	RETURN_MORE;
+    }
+
+    ret = dplot_face_curves(&info);
+    if (ret == GED_ERROR) {
+	RETURN_ERROR;
+    } else if (ret == GED_MORE) {
+	RETURN_MORE;
+    }
+
+    ret = dplot_split_faces(&info);
+    if (ret == GED_ERROR) {
+	RETURN_ERROR;
+    } else if (ret == GED_MORE) {
+	RETURN_MORE;
+    }
+
+    ret = dplot_linked_curves(&info);
+    if (ret == GED_ERROR) {
+	RETURN_ERROR;
+    } else if (ret == GED_MORE) {
+	RETURN_MORE;
+    }
+
+    info.mode = DPLOT_INITIAL;
+    CLEANUP;
+
+    return GED_OK;
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -0,0 +1,103 @@
+%include {
+#include "bu.h"
+#include <assert.h>
+#include "dplot_reader.h"
+}
+
+%destructor statement_list {
+    $$.n = data->ssx_count;
+}
+
+%stack_overflow {
+    fprintf(stderr, "Parser stack overflow. Last token was %d\n", yypMinor->yy0.n);
+}
+
+%syntax_error {
+    fprintf(stderr, "Syntax error. Last token (type %d) was %d. Ignoring "
+	    "input tokens until valid input is found.\n", yymajor, yyminor.yy0.n);
+}
+
+%parse_failure {
+    fprintf(stderr, "Parsing failed. Parser has reset.\n");
+}
+
+%token_type {token_t}
+%extra_argument {struct dplot_data *data}
+
+%token_prefix {TOKEN_}
+
+%type count {int}
+%type index {int}
+
+start_symbol ::= statement_list.
+
+statement_list ::= /* empty */.
+statement_list ::= statement_list statement.
+
+statement ::= surface_counts EOL.
+statement ::= ssx_info EOL.
+statement ::= linked_curve_count EOL.
+statement ::= split_face_count EOL.
+statement ::= split_face_info EOL.
+
+surface_counts ::= SURFACES count(C1) count(C2). {
+    data->brep1_surface_count = C1;
+    data->brep2_surface_count = C2;
+}
+
+ssx_info ::= SSX index(SURF_IDX1) index(SURF_IDX2) count(EVENTS) count(CRV_COUNT1) count(CRV_COUNT2) count(B1_ISO) count(ISO) isocsx_event_counts. {
+    struct ssx *ssx;
+    BU_GET(ssx, struct ssx);
+
+    ssx->brep1_surface = SURF_IDX1;
+    ssx->brep2_surface = SURF_IDX2;
+    ssx->final_curve_events = EVENTS;
+    ssx->face1_clipped_curves = CRV_COUNT1;
+    ssx->face2_clipped_curves = CRV_COUNT2;
+    ssx->intersecting_brep1_isocurves = B1_ISO;
+    ssx->intersecting_isocurves = ISO;
+    BU_LIST_INIT(&ssx->isocsx_list);
+
+    /* move the list of event counts from the tmp location to this ssx
+     * instance
+     */
+    BU_LIST_APPEND_LIST(&ssx->isocsx_list, &data->isocsx_list);
+
+    BU_LIST_PUSH(&data->ssx_list, &ssx->l);
+    ++data->ssx_count;
+}
+
+isocsx_event_counts ::= /* empty */.
+
+isocsx_event_counts ::= isocsx_event_counts count(C). {
+    struct isocsx *isocsx;
+    BU_GET(isocsx, struct isocsx);
+
+    isocsx->events = C;
+    BU_LIST_INIT(&isocsx->l);
+
+    BU_LIST_PUSH(&data->isocsx_list, &isocsx->l);
+}
+
+linked_curve_count ::= LINKED_CURVES count(COUNT). {
+    data->linked_curve_count = COUNT;
+}
+
+split_face_count ::= SPLIT_FACES count(SPLIT_FACES). {
+    data->split_face_count = SPLIT_FACES;
+    data->face = (struct split_face *)bu_malloc(SPLIT_FACES *
+	sizeof(struct split_face), "split face array");
+}
+
+split_face_info ::= SPLIT_FACE index(I) count(OUTER_COUNT) count(INNER_COUNT). {
+    data->face[I].outerloop_curves = OUTER_COUNT;
+    data->face[I].innerloop_curves = INNER_COUNT;
+}
+
+count(C) ::= NONNEGATIVE_INT(N). {
+    C = N.n;
+}
+
+index(I) ::= NONNEGATIVE_INT(N). {
+    I = N.n;
+}
@@ -0,0 +1,87 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "dplot_reader.h"
+#include "bu.h"
+
+HIDDEN void *
+dplot_malloc(size_t s) {
+    return bu_malloc(s, "dplot_malloc");
+}
+
+HIDDEN void
+dplot_free(void *p) {
+    bu_free(p, "dplot_free");
+}
+
+int
+main(int argc, char *argv[])
+{
+    int i, token_id;
+    FILE *input_file;
+    perplex_t scanner;
+    void *parser;
+    struct ssx *curr;
+    struct dplot_data data;
+
+    if (argc != 2) {
+	bu_exit(1, "usage: %s input\n", argv[0]);
+    }
+
+    /* initialize scanner */
+    input_file = fopen(argv[1], "r");
+    if (!input_file) {
+	bu_exit(2, "error: couldn't open \"%s\" for reading.\n", argv[1]);
+    }
+    scanner = perplexFileScanner(input_file);
+
+    data.ssx_count = 0;
+    BU_LIST_INIT(&data.ssx_list);
+    perplexSetExtra(scanner, (void *)&data);
+
+    /* initialize parser */
+    parser = ParseAlloc(dplot_malloc);
+
+    /* parse */
+    while ((token_id = yylex(scanner)) != YYEOF) {
+	Parse(parser, token_id, data.token_data, &data);
+    }
+    Parse(parser, 0, data.token_data, &data);
+
+    /* move list data to array */
+    if (data.ssx_count > 0) {
+	data.ssx = (struct ssx *)bu_malloc(
+		sizeof(struct ssx) * data.ssx_count, "ssx array");
+
+	i = data.ssx_count - 1;
+	while (BU_LIST_WHILE(curr, ssx, &data.ssx_list)) {
+	    data.ssx[i--] = *curr;
+	    BU_LIST_DEQUEUE(&curr->l);
+	    BU_PUT(curr, struct ssx);
+	}
+    }
+
+    /* print data */
+    bu_log("brep1: %d surfaces\n", data.brep1_surface_count);
+    bu_log("brep2: %d surfaces\n", data.brep2_surface_count);
+
+    for (i = 0; i < data.ssx_count; ++i) {
+	bu_log("%d curve events between surface 1-%d and 2-%d\n",
+		data.ssx[i].final_curve_events, data.ssx[i].brep1_surface,
+		data.ssx[i].brep2_surface);
+	bu_log("\t%d brep1 isocurves intersected brep2\n",
+		data.ssx[i].intersecting_brep1_isocurves);
+	bu_log("\t%d brep2 isocurves intersected brep1\n",
+		data.ssx[i].intersecting_isocurves -
+		data.ssx[i].intersecting_brep1_isocurves);
+    }
+
+    /* clean up */
+    if (data.ssx) {
+	bu_free(data.ssx, "ssx array");
+    }
+    ParseFree(parser, dplot_free);
+    perplexFree(scanner);
+    fclose(input_file);
+
+    return 0;
+}
@@ -0,0 +1,56 @@
+#ifndef DPLOT_PARSER
+#define DPLOT_PARSER
+
+#define PERPLEX_ON_ENTER struct dplot_data *data = (struct dplot_data *)yyextra;
+
+#include "bu.h"
+#include "dplot_parser.h"
+#include "dplot_scanner.h"
+
+typedef struct {
+    int n;
+} token_t;
+
+struct isocsx {
+    struct bu_list l;
+    int events;
+};
+
+struct ssx {
+    struct bu_list l;
+    int brep1_surface;
+    int brep2_surface;
+    int final_curve_events;
+    int face1_clipped_curves;
+    int face2_clipped_curves;
+    int intersecting_brep1_isocurves;
+    int intersecting_isocurves;
+    struct bu_list isocsx_list; /* struct isocsx list */
+    int *isocsx_events;
+};
+
+struct split_face {
+    int outerloop_curves;
+    int innerloop_curves;
+};
+
+struct dplot_data {
+    token_t token_data;
+    int brep1_surface_count;
+    int brep2_surface_count;
+    int linked_curve_count;
+    int ssx_count;
+    int split_face_count;
+    struct bu_list ssx_list; /* struct ssx list */
+    struct ssx *ssx;
+    struct bu_list isocsx_list; /* struct ssx list */
+    struct split_face *face;
+};
+
+/* lemon prototypes */
+void *ParseAlloc(void *(*mallocProc)(size_t));
+void ParseFree(void *parser, void (*freeProc)(void *));
+void Parse(void *yyp, int yymajor, token_t tokenData, struct dplot_data *data);
+void ParseTrace(FILE *fp, char *s);
+
+#endif
@@ -0,0 +1,42 @@
+#include "bu.h"
+#include "dplot_reader.h"
+enum {INITIAL, normal};
+%%
+<> => normal { continue; }
+
+<normal> {
+"surfaces" {
+    return TOKEN_SURFACES;
+}
+
+"ssx" {
+    return TOKEN_SSX;
+}
+
+"linkedcurves" {
+    return TOKEN_LINKED_CURVES;
+}
+
+"splitfaces" {
+    return TOKEN_SPLIT_FACES;
+}
+
+"splitface" {
+    return TOKEN_SPLIT_FACE;
+}
+
+[0-9]|[1-9][0-9]+ {
+    data->token_data.n = 0;
+    bu_sscanf(yytext, "%d", &data->token_data.n);
+    return TOKEN_NONNEGATIVE_INT;
+}
+
+[\n]+ {
+    return TOKEN_EOL;
+}
+
+[^] {
+    /* ignore */
+    continue;
+}
+}
@@ -111,6 +111,8 @@ GED_EXPORT extern int brep_intersect_curve_curve(struct rt_db_internal *intern1,
 GED_EXPORT extern int brep_intersect_curve_surface(struct rt_db_internal *intern1, struct rt_db_internal *intern2, int i, int j);
 GED_EXPORT extern int brep_intersect_surface_surface(struct rt_db_internal *intern1, struct rt_db_internal *intern2, int i, int j, struct bn_vlblock *vbp);

+GED_EXPORT extern int ged_dplot_core(struct ged *gedp, int argc, const char *argv[]);
+
 using namespace brlcad;
 void
 plotface(const ON_BrepFace &face, struct bn_vlblock *vbp, int plotres, bool dim3d, const int red, const int green, const int blue);
@@ -127,6 +127,7 @@ GED_CMD(debugnmg)
 GED_CMD(decompose)
 GED_CMD(delay)
 GED_CMD(dir2ae)
+GED_CMD(dplot)
 GED_CMD(draw)
 GED_CMD(dsp)
 GED_CMD(dump)
@@ -13,6 +13,7 @@ CMAKEFILES(
  lint_tests.g
  obb_tests.g
  ged_tests.g
+  pnts_in_out.g
  )

 CMAKEFILES(CMakeLists.txt)
@@ -180,6 +180,12 @@ HIDDEN int to_data_vZ(struct ged *gedp,
 		      ged_func_ptr func,
 		      const char *usage,
 		      int maxargs);
+HIDDEN int to_dplot(struct ged *gedp,
+		       int argc,
+		       const char *argv[],
+		       ged_func_ptr func,
+		       const char *usage,
+		       int maxargs);
 HIDDEN int to_dlist_on(struct ged *gedp,
 		       int argc,
 		       const char *argv[],
@@ -604,6 +610,7 @@ static struct to_cmdtab ged_cmds[] = {
    {"debugnmg",	(char *)0, TO_UNLIMITED, to_pass_through_func, ged_debugnmg},
    {"decompose",	(char *)0, TO_UNLIMITED, to_pass_through_func, ged_decompose},
    {"delay",	(char *)0, TO_UNLIMITED, to_pass_through_func, ged_delay},
+    {"dplot",	"dplot_logfile", 1, to_dplot, ged_dplot},
    {"metaball_delete_pnt",	(char *)0, TO_UNLIMITED, to_pass_through_func, ged_metaball_delete_pnt},
    {"pipe_delete_pnt",	(char *)0, TO_UNLIMITED, to_pass_through_func, ged_pipe_delete_pnt},
    {"dir2ae",	(char *)0, TO_UNLIMITED, to_pass_through_func, ged_dir2ae},
@@ -3217,6 +3224,144 @@ to_dlist_on(struct ged *gedp,
    return GED_OK;
 }

+HIDDEN int
+to_dplot(struct ged *gedp,
+	    int argc,
+	    const char *argv[],
+	    ged_func_ptr func,
+	    const char *UNUSED(usage),
+	    int UNUSED(maxargs))
+{
+    int ac;
+    int ret;
+    char *av[256];
+    struct bu_vls callback_cmd = BU_VLS_INIT_ZERO;
+    struct bu_vls temp = BU_VLS_INIT_ZERO;
+    struct bu_vls result_copy = BU_VLS_INIT_ZERO;
+    struct bview *gdvp;
+    struct tclcad_ged_data *tgd = (struct tclcad_ged_data *)current_top->to_gedp->u_data;
+    char *who_av[3] = {"who", "b", NULL};
+    int first = 1;
+    int aflag = 0;
+
+    /* copy all args */
+    ac = argc;
+    for (int i = 0; i < ac; ++i)
+	av[i] = bu_strdup((char *)argv[i]);
+    av[ac] = (char *)0;
+
+    /* check for displayed objects */
+    ret = ged_who(gedp, 2, (const char **)who_av);
+    if (ret == GED_OK && strlen(bu_vls_addr(gedp->ged_result_str)) == 0)
+	aflag = 1;
+    bu_vls_trunc(gedp->ged_result_str, 0);
+
+    while ((ret = (*func)(gedp, ac, (const char **)av)) & GED_MORE) {
+	int ac_more;
+	const char **avmp;
+	const char **av_more = NULL;
+
+	/* save result string */
+	bu_vls_substr(&result_copy, gedp->ged_result_str, 0, bu_vls_strlen(gedp->ged_result_str));
+	bu_vls_trunc(gedp->ged_result_str, 0);
+
+	ret = ged_who(gedp, 1, (const char **)who_av);
+	if (ret == GED_OK && strlen(bu_vls_addr(gedp->ged_result_str)) == 0)
+	    aflag = 1;
+
+	bu_vls_trunc(gedp->ged_result_str, 0);
+
+	for (size_t i = 0; i < BU_PTBL_LEN(&current_top->to_gedp->ged_views); i++) {
+	    gdvp = (struct bview *)BU_PTBL_GET(&current_top->to_gedp->ged_views, i);
+	    if (to_is_viewable(gdvp)) {
+		gedp->ged_gvp->gv_x_samples = dm_get_width((struct dm *)gdvp->dmp);
+		gedp->ged_gvp->gv_y_samples = dm_get_height((struct dm *)gdvp->dmp);
+	    }
+	}
+
+	if (first && aflag) {
+	    first = 0;
+	    to_autoview_all_views(current_top);
+	} else {
+	    to_refresh_all_views(current_top);
+	}
+	/* restore result string */
+	bu_vls_substr(gedp->ged_result_str, &result_copy, 0, bu_vls_strlen(&result_copy));
+	bu_vls_free(&result_copy);
+
+	if (0 < bu_vls_strlen(&tgd->go_more_args_callback)) {
+	    bu_vls_trunc(&callback_cmd, 0);
+	    bu_vls_printf(&callback_cmd, "%s [string range {%s} 0 end]",
+			  bu_vls_addr(&tgd->go_more_args_callback),
+			  bu_vls_addr(gedp->ged_result_str));
+
+	    if (Tcl_Eval(current_top->to_interp, bu_vls_addr(&callback_cmd)) != TCL_OK) {
+		bu_vls_trunc(gedp->ged_result_str, 0);
+		bu_vls_printf(gedp->ged_result_str, "%s", Tcl_GetStringResult(current_top->to_interp));
+		Tcl_ResetResult(current_top->to_interp);
+		return GED_ERROR;
+	    }
+
+	    bu_vls_trunc(&temp, 0);
+	    bu_vls_printf(&temp, "%s", Tcl_GetStringResult(current_top->to_interp));
+	    Tcl_ResetResult(current_top->to_interp);
+	} else {
+	    bu_log("\r%s", bu_vls_addr(gedp->ged_result_str));
+	    bu_vls_trunc(&temp, 0);
+	    if (bu_vls_gets(&temp, stdin) < 0) {
+		break;
+	    }
+	}
+
+	if (Tcl_SplitList(current_top->to_interp, bu_vls_addr(&temp), &ac_more, &av_more) != TCL_OK) {
+	    continue;
+	}
+
+	if (ac_more < 1) {
+	    /* space has still been allocated */
+	    Tcl_Free((char *)av_more);
+
+	    continue;
+	}
+
+	/* skip first element if empty */
+	avmp = av_more;
+	if (*avmp[0] == '\0') {
+	    --ac_more;
+	    ++avmp;
+	}
+
+	/* ignore last element if empty */
+	if (*avmp[ac_more-1] == '\0')
+	    --ac_more;
+
+	/* copy additional args */
+	for (int i = 0; i < ac_more; ++i)
+	    av[ac++] = bu_strdup(avmp[i]);
+	av[ac+1] = (char *)0;
+
+	Tcl_Free((char *)av_more);
+    }
+
+    for (size_t i = 0; i < BU_PTBL_LEN(&current_top->to_gedp->ged_views); i++) {
+	gdvp = (struct bview *)BU_PTBL_GET(&current_top->to_gedp->ged_views, i);
+	if (to_is_viewable(gdvp)) {
+	    gedp->ged_gvp->gv_x_samples = dm_get_width((struct dm *)gdvp->dmp);
+	    gedp->ged_gvp->gv_y_samples = dm_get_height((struct dm *)gdvp->dmp);
+	}
+    }
+    to_refresh_all_views(current_top);
+
+    bu_vls_free(&callback_cmd);
+    bu_vls_free(&temp);
+
+    bu_vls_printf(gedp->ged_result_str, "BUILT_BY_MORE_ARGS");
+    for (int i = 0; i < ac; ++i) {
+	bu_vls_printf(gedp->ged_result_str, "%s ", av[i]);
+	bu_free((void *)av[i], "to_more_args_func");
+    }
+    return GED_OK;
+}

 HIDDEN int
 to_fontsize(struct ged *gedp,
@@ -212,6 +212,7 @@ package provide cadwidgets::Ged 1.0
 	method delay {args}
 	method dir2ae {args}
 	method dlist_on {args}
+	method dplot {args}
 	method draw {args}
 	method draw_ray {_start _partitions}
 	method dsp {args}
@@ -1636,6 +1637,10 @@ package provide cadwidgets::Ged 1.0
    eval $mGed dlist_on $args
 }

+::itcl::body cadwidgets::Ged::dplot {args} {
+    eval $mGed dplot $args
+}
+
 ::itcl::body cadwidgets::Ged::draw {args} {
    set mRayNeedGettrees 1

@@ -6259,6 +6264,7 @@ package provide cadwidgets::Ged 1.0
    $help add decompose		{{nmg_solid [prefix]}	{decompose nmg_solid into maximally connected shells}}
    $help add delay		{{sec usec} {delay processing for the specified amount of time}}
    $help add dir2ae		{{az el} {returns a direction vector given the azimuth and elevation}}
+    $help add dplot		{{dplot_log} {draws plotfiles from the given logfile}}
    $help add draw		{{"-C#/#/# <objects>"} {draw objects}}
    $help add dsp		{{obj [command]} {work with DSP primitives}}
    $help add dump		{{file} {write current state of database object to file}}
Author	SHA1	Message	Date
Clifford Yapp	e5403a286d	Changes merged, delete brep-debug (svn branch delete) svn:revision:77232 svn:branch:brep-debug svn:account:starseeker	2020-09-25 16:52:15 +00:00
Clifford Yapp	4198b5d007	Add libged symbol for libtclcad's use of ged_dplot svn:revision:77229 svn:branch:brep-debug svn:account:starseeker	2020-09-25 16:29:19 +00:00
Clifford Yapp	827c71210e	Port the dplot logic from tclcad_obj.c to the new setup. svn:revision:77228 svn:branch:brep-debug svn:account:starseeker	2020-09-25 16:28:34 +00:00
Clifford Yapp	24a06afaf6	dplot command needs to make the transition to the new ged command structure. svn:revision:77227 svn:branch:brep-debug svn:account:starseeker	2020-09-25 15:57:21 +00:00
Clifford Yapp	3536c92db7	Move bool development guide svn:revision:77226 svn:branch:brep-debug svn:account:starseeker	2020-09-25 15:35:15 +00:00
Clifford Yapp	25bcf35115	Move bool development guide (preliminary file move commit)	2020-09-25 15:35:15 +00:00
Clifford Yapp	4cbc21b73a	Merge trunk r76151 through r77224 svn:revision:77225 svn:branch:brep-debug svn:account:starseeker	2020-09-25 14:49:43 +00:00
Clifford Yapp	abe3f08bb6	Merge analyze command rework from branch, add doc update and NEWS items. Immediate functionality is to use geometric inside/outside testing to peform boolean ops on point clouds and generate new point sets, but should be implemented in such a way as to allow for expansion of additional capabilities. svn:revision:77222 svn:branch:trunk svn:account:starseeker	2020-09-25 13:41:24 +00:00
Clifford Yapp	a20b076073	Merge changes from trunk r71095 through r76151 svn:revision:76153 svn:branch:brep-debug svn:account:starseeker	2020-06-18 15:37:23 +00:00
Christopher Sean Morrison	a3064251ee	clarify svn:revision:72643 svn:branch:brep-debug svn:account:brlcad	2019-03-15 20:34:00 +00:00