# Camera/motion tracking
option(WITH_LIBMV "Enable libmv structure from motion library" ON)
+# Mesh boolean lib
+option(WITH_CARVE "Enable Carve library to handle mesh boolean operations" ON)
+
# Misc
option(WITH_INPUT_NDOF "Enable NDOF input devices (SpaceNavigator and friends)" ON)
option(WITH_RAYOPTIMIZATION "Enable use of SIMD (SSE) optimizations for the raytracer" ON)
set(WITH_HEADLESS ON)
endif()
-# auto enable openimageio and boost for cycles
+# auto enable openimageio for cycles
if(WITH_CYCLES)
set(WITH_OPENIMAGEIO ON)
+endif()
+
+# auto enable boost for cycles and carve
+if(WITH_CYCLES OR WITH_CARVE)
set(WITH_BOOST ON)
endif()
'WITH_BF_CYCLES', 'WITH_BF_CYCLES_CUDA_BINARIES' 'BF_CYCLES_CUDA_NVCC', 'BF_CYCLES_CUDA_NVCC', 'WITH_BF_CYCLES_CUDA_THREADED_COMPILE',
'WITH_BF_OIIO', 'WITH_BF_STATICOIIO', 'BF_OIIO', 'BF_OIIO_INC', 'BF_OIIO_LIB', 'BF_OIIO_LIB_STATIC', 'BF_OIIO_LIBPATH',
'WITH_BF_BOOST', 'WITH_BF_STATICBOOST', 'BF_BOOST', 'BF_BOOST_INC', 'BF_BOOST_LIB', 'BF_BOOST_LIB_STATIC', 'BF_BOOST_LIBPATH',
- 'WITH_BF_LIBMV'
+ 'WITH_BF_LIBMV', 'WITH_BF_CARVE'
]
# Have options here that scons expects to be lists
(BoolVariable('WITH_BF_LZO', 'Enable fast LZO pointcache compression', True)),
(BoolVariable('WITH_BF_LZMA', 'Enable best LZMA pointcache compression', True)),
+ (BoolVariable('WITH_BF_CARVE', 'Enable carve library for mesh boolean operations', True)),
(BoolVariable('WITH_BF_LIBMV', 'Enable libmv structure from motion library', True)),
if(WITH_LIBMV)
add_subdirectory(libmv)
endif()
+
+if(WITH_CARVE)
+ add_subdirectory(carve)
+endif()
if env['WITH_BF_LIBMV']:
SConscript(['libmv/SConscript'])
+
+if env['WITH_BF_CARVE']:
+ SConscript(['carve/SConscript'])
--- /dev/null
+# ***** BEGIN GPL LICENSE BLOCK *****
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+#
+# The Original Code is Copyright (C) 2006, Blender Foundation
+# All rights reserved.
+#
+# The Original Code is: all of this file.
+#
+# Contributor(s): Jacques Beaurai, Erwin Coumans
+#
+# ***** END GPL LICENSE BLOCK *****
+
+# NOTE: This file is automatically generated by bundle.sh script
+# If you're doing changes in this file, please update template
+# in that script too
+
+set(INC
+ include
+)
+
+set(INC_SYS
+)
+
+set(SRC
+ lib/carve.cpp
+ lib/mesh.cpp
+ lib/intersect_group.cpp
+ lib/intersect_classify_edge.cpp
+ lib/intersect_classify_group.cpp
+ lib/polyhedron.cpp
+ lib/geom3d.cpp
+ lib/polyline.cpp
+ lib/csg_collector.cpp
+ lib/triangulator.cpp
+ lib/intersect_face_division.cpp
+ lib/intersect_half_classify_group.cpp
+ lib/edge.cpp
+ lib/math.cpp
+ lib/geom2d.cpp
+ lib/tag.cpp
+ lib/intersection.cpp
+ lib/convex_hull.cpp
+ lib/csg.cpp
+ lib/intersect.cpp
+ lib/face.cpp
+ lib/pointset.cpp
+ lib/timing.cpp
+ lib/octree.cpp
+ lib/aabb.cpp
+ lib/intersect_debug.cpp
+
+ lib/intersect_classify_common.hpp
+ lib/csg_data.hpp
+ lib/csg_collector.hpp
+ lib/intersect_common.hpp
+ lib/intersect_classify_common_impl.hpp
+ lib/csg_detail.hpp
+ lib/intersect_debug.hpp
+
+ include/carve/polyhedron_decl.hpp
+ include/carve/geom2d.hpp
+ include/carve/exact.hpp
+ include/carve/triangulator_impl.hpp
+ include/carve/collection.hpp
+ include/carve/pointset.hpp
+ include/carve/djset.hpp
+ include/carve/kd_node.hpp
+ include/carve/polyline.hpp
+ include/carve/polyline_iter.hpp
+ include/carve/geom3d.hpp
+ include/carve/edge_decl.hpp
+ include/carve/face_decl.hpp
+ include/carve/aabb_impl.hpp
+ include/carve/colour.hpp
+ include/carve/pointset_iter.hpp
+ include/carve/polyline_decl.hpp
+ include/carve/rescale.hpp
+ include/carve/mesh_impl.hpp
+ include/carve/classification.hpp
+ include/carve/util.hpp
+ include/carve/triangulator.hpp
+ include/carve/polyhedron_base.hpp
+ include/carve/rtree.hpp
+ include/carve/math.hpp
+ include/carve/math_constants.hpp
+ include/carve/octree_decl.hpp
+ include/carve/input.hpp
+ include/carve/mesh_ops.hpp
+ include/carve/debug_hooks.hpp
+ include/carve/mesh_simplify.hpp
+ include/carve/interpolator.hpp
+ include/carve/poly_decl.hpp
+ include/carve/csg.hpp
+ include/carve/mesh.hpp
+ include/carve/carve.hpp
+ include/carve/gnu_cxx.h
+ include/carve/polyhedron_impl.hpp
+ include/carve/poly_impl.hpp
+ include/carve/aabb.hpp
+ include/carve/convex_hull.hpp
+ include/carve/vertex_decl.hpp
+ include/carve/win32.h
+ include/carve/edge_impl.hpp
+ include/carve/tag.hpp
+ include/carve/tree.hpp
+ include/carve/heap.hpp
+ include/carve/matrix.hpp
+ include/carve/poly.hpp
+ include/carve/vector.hpp
+ include/carve/intersection.hpp
+ include/carve/faceloop.hpp
+ include/carve/geom_impl.hpp
+ include/carve/octree_impl.hpp
+ include/carve/spacetree.hpp
+ include/carve/collection/unordered/std_impl.hpp
+ include/carve/collection/unordered/tr1_impl.hpp
+ include/carve/collection/unordered/libstdcpp_impl.hpp
+ include/carve/collection/unordered/boost_impl.hpp
+ include/carve/collection/unordered/vcpp_impl.hpp
+ include/carve/collection/unordered/fallback_impl.hpp
+ include/carve/collection/unordered.hpp
+ include/carve/face_impl.hpp
+ include/carve/pointset_impl.hpp
+ include/carve/cbrt.h
+ include/carve/vcpp_config.h
+ include/carve/geom.hpp
+ include/carve/vertex_impl.hpp
+ include/carve/polyline_impl.hpp
+ include/carve/pointset_decl.hpp
+ include/carve/timing.hpp
+ include/carve/csg_triangulator.hpp
+ include/carve/iobj.hpp
+ include/carve/collection_types.hpp
+)
+
+if(WITH_BOOST)
+ if(NOT MSVC)
+ # Boost is setting as preferred collections library in the Carve code when using MSVC compiler
+ add_definitions(
+ -DHAVE_BOOST_UNORDERED_COLLECTIONS
+ )
+ endif()
+
+ add_definitions(
+ -DCARVE_SYSTEM_BOOST
+ )
+
+ list(APPEND INC
+ ${BOOST_INCLUDE_DIR}
+ )
+endif()
+
+blender_add_lib(extern_carve "${SRC}" "${INC}" "${INC_SYS}")
--- /dev/null
+ GNU GENERAL PUBLIC LICENSE
+
+ The Qt GUI Toolkit is Copyright (C) 1994-2008 Trolltech ASA.
+
+ You may use, distribute and copy the Qt GUI Toolkit under the terms of
+ GNU General Public License version 2, which is displayed below.
+
+-------------------------------------------------------------------------
+
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+ 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.) You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. You must make sure that they, too, receive or can get the
+source code. And you must show them these terms so they know their
+rights.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+\f
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+\f
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+\f
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+\f
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+\f
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs. If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library. If this is what you want to do, use the GNU Library General
+Public License instead of this License.
+
+-------------------------------------------------------------------------
+
+In addition, as a special exception, Trolltech gives permission to link the
+code of its release of Qt with the OpenSSL project's "OpenSSL" library (or
+modified versions of it that use the same license as the "OpenSSL"
+library), and distribute the linked executables. You must comply with the GNU
+General Public License version 2 or the GNU General Public License version 3
+in all respects for all of the code used other than the "OpenSSL" code. If
+you modify this file, you may extend this exception to your version of the
+file, but you are not obligated to do so. If you do not wish to do so,
+delete this exception statement from your version of this file.
--- /dev/null
+#!/usr/bin/python
+
+# NOTE: This file is automatically generated by bundle.sh script
+# If you're doing changes in this file, please update template
+# in that script too
+
+Import ('env')
+
+sources = env.Glob('lib/*.cpp')
+
+defs = []
+incs = ['include']
+
+if env['WITH_BF_BOOST']:
+ if env['OURPLATFORM'] not in ('win32-vc', 'win64-vc'):
+ # Boost is setting as preferred collections library in the Carve code when using MSVC compiler
+ defs.append('HAVE_BOOST_UNORDERED_COLLECTIONS')
+
+ defs.append('CARVE_SYSTEM_BOOST')
+ incs.append(env['BF_BOOST_INC'])
+
+env.BlenderLib ('extern_carve', Split(sources), incs, defs, libtype=['extern'], priority=[40] )
--- /dev/null
+#!/bin/sh
+
+if [ -d ./.svn ]; then
+ echo "This script is supposed to work only when using git-svn"
+ exit 1
+fi
+
+tmp=`mktemp -d`
+
+hg clone https://code.google.com/p/carve/ $tmp/carve
+
+for p in `cat ./patches/series`; do
+ echo "Applying patch $p..."
+ cat ./patches/$p | patch -d $tmp/carve -p1
+done
+
+rm -rf include
+rm -rf lib
+
+cat "files.txt" | while f=`line`; do
+ mkdir -p `dirname $f`
+ cp $tmp/carve/$f $f
+done
+
+rm -rf $tmp
+
+sources=`find ./lib -type f -iname '*.cc' -or -iname '*.cpp' -or -iname '*.c' | sed -r 's/^\.\//\t/'`
+headers=`find ./lib -type f -iname '*.h' -or -iname '*.hpp' | sed -r 's/^\.\//\t/'`
+includes=`find ./include -type f -iname '*.h' -or -iname '*.hpp' | sed -r 's/^\.\//\t/'`
+
+mkdir -p include/carve/external/boost
+cp patches/files/random.hpp include/carve/external/boost/random.hpp
+cp patches/files/config.h include/carve/config.h
+
+cat > CMakeLists.txt << EOF
+# ***** BEGIN GPL LICENSE BLOCK *****
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+#
+# The Original Code is Copyright (C) 2006, Blender Foundation
+# All rights reserved.
+#
+# The Original Code is: all of this file.
+#
+# Contributor(s): Jacques Beaurai, Erwin Coumans
+#
+# ***** END GPL LICENSE BLOCK *****
+
+# NOTE: This file is automatically generated by bundle.sh script
+# If you're doing changes in this file, please update template
+# in that script too
+
+set(INC
+ include
+)
+
+set(INC_SYS
+)
+
+set(SRC
+${sources}
+
+${headers}
+
+${includes}
+)
+
+if(WITH_BOOST)
+ if(NOT MSVC)
+ # Boost is setting as preferred collections library in the Carve code when using MSVC compiler
+ add_definitions(
+ -DHAVE_BOOST_UNORDERED_COLLECTIONS
+ )
+ endif()
+
+ add_definitions(
+ -DCARVE_SYSTEM_BOOST
+ )
+
+ list(APPEND INC
+ \${BOOST_INCLUDE_DIR}
+ )
+endif()
+
+blender_add_lib(extern_carve "\${SRC}" "\${INC}" "\${INC_SYS}")
+EOF
+
+cat > SConscript << EOF
+#!/usr/bin/python
+
+# NOTE: This file is automatically generated by bundle.sh script
+# If you're doing changes in this file, please update template
+# in that script too
+
+Import ('env')
+
+sources = env.Glob('lib/*.cpp')
+
+defs = []
+incs = ['include']
+
+if env['WITH_BF_BOOST']:
+ if env['OURPLATFORM'] not in ('win32-vc', 'win64-vc'):
+ # Boost is setting as preferred collections library in the Carve code when using MSVC compiler
+ defs.append('HAVE_BOOST_UNORDERED_COLLECTIONS')
+
+ defs.append('CARVE_SYSTEM_BOOST')
+ incs.append(env['BF_BOOST_INC'])
+
+env.BlenderLib ('extern_carve', Split(sources), incs, defs, libtype=['extern'], priority=[40] )
+EOF
--- /dev/null
+include/carve/polyhedron_decl.hpp
+include/carve/geom2d.hpp
+include/carve/exact.hpp
+include/carve/triangulator_impl.hpp
+include/carve/collection.hpp
+include/carve/pointset.hpp
+include/carve/djset.hpp
+include/carve/kd_node.hpp
+include/carve/polyline.hpp
+include/carve/polyline_iter.hpp
+include/carve/geom3d.hpp
+include/carve/edge_decl.hpp
+include/carve/face_decl.hpp
+include/carve/aabb_impl.hpp
+include/carve/colour.hpp
+include/carve/pointset_iter.hpp
+include/carve/polyline_decl.hpp
+include/carve/rescale.hpp
+include/carve/mesh_impl.hpp
+include/carve/classification.hpp
+include/carve/util.hpp
+include/carve/triangulator.hpp
+include/carve/polyhedron_base.hpp
+include/carve/rtree.hpp
+include/carve/math.hpp
+include/carve/math_constants.hpp
+include/carve/octree_decl.hpp
+include/carve/input.hpp
+include/carve/mesh_ops.hpp
+include/carve/debug_hooks.hpp
+include/carve/mesh_simplify.hpp
+include/carve/interpolator.hpp
+include/carve/poly_decl.hpp
+include/carve/csg.hpp
+include/carve/mesh.hpp
+include/carve/carve.hpp
+include/carve/gnu_cxx.h
+include/carve/polyhedron_impl.hpp
+include/carve/poly_impl.hpp
+include/carve/aabb.hpp
+include/carve/convex_hull.hpp
+include/carve/vertex_decl.hpp
+include/carve/win32.h
+include/carve/edge_impl.hpp
+include/carve/tag.hpp
+include/carve/tree.hpp
+include/carve/heap.hpp
+include/carve/matrix.hpp
+include/carve/poly.hpp
+include/carve/vector.hpp
+include/carve/intersection.hpp
+include/carve/faceloop.hpp
+include/carve/geom_impl.hpp
+include/carve/octree_impl.hpp
+include/carve/spacetree.hpp
+include/carve/collection/unordered/std_impl.hpp
+include/carve/collection/unordered/tr1_impl.hpp
+include/carve/collection/unordered/libstdcpp_impl.hpp
+include/carve/collection/unordered/boost_impl.hpp
+include/carve/collection/unordered/vcpp_impl.hpp
+include/carve/collection/unordered/fallback_impl.hpp
+include/carve/collection/unordered.hpp
+include/carve/face_impl.hpp
+include/carve/pointset_impl.hpp
+include/carve/cbrt.h
+include/carve/vcpp_config.h
+include/carve/geom.hpp
+include/carve/vertex_impl.hpp
+include/carve/polyline_impl.hpp
+include/carve/pointset_decl.hpp
+include/carve/timing.hpp
+include/carve/csg_triangulator.hpp
+include/carve/iobj.hpp
+include/carve/collection_types.hpp
+lib/carve.cpp
+lib/mesh.cpp
+lib/intersect_group.cpp
+lib/intersect_classify_common.hpp
+lib/intersect_classify_edge.cpp
+lib/intersect_classify_group.cpp
+lib/csg_data.hpp
+lib/polyhedron.cpp
+lib/csg_collector.hpp
+lib/geom3d.cpp
+lib/polyline.cpp
+lib/csg_collector.cpp
+lib/triangulator.cpp
+lib/intersect_face_division.cpp
+lib/intersect_half_classify_group.cpp
+lib/edge.cpp
+lib/math.cpp
+lib/geom2d.cpp
+lib/tag.cpp
+lib/intersection.cpp
+lib/convex_hull.cpp
+lib/intersect_common.hpp
+lib/intersect_classify_common_impl.hpp
+lib/csg.cpp
+lib/intersect.cpp
+lib/csg_detail.hpp
+lib/face.cpp
+lib/pointset.cpp
+lib/timing.cpp
+lib/octree.cpp
+lib/aabb.cpp
+lib/intersect_debug.hpp
+lib/intersect_debug.cpp
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/geom.hpp>
+
+#include <vector>
+
+namespace carve {
+ namespace geom {
+
+
+
+ // n-dimensional AABB
+ template<unsigned ndim>
+ struct aabb {
+ typedef vector<ndim> vector_t;
+ typedef aabb<ndim> aabb_t;
+
+ vector_t pos; // the centre of the AABB
+ vector_t extent; // the extent of the AABB - the vector from the centre to the maximal vertex.
+
+ void empty();
+
+ bool isEmpty() const;
+
+ void fit(const vector_t &v1);
+ void fit(const vector_t &v1, const vector_t &v2);
+ void fit(const vector_t &v1, const vector_t &v2, const vector_t &v3);
+
+ template<typename iter_t, typename value_type>
+ void _fit(iter_t begin, iter_t end, value_type);
+
+ template<typename iter_t>
+ void _fit(iter_t begin, iter_t end, vector_t);
+
+ template<typename iter_t>
+ void _fit(iter_t begin, iter_t end, aabb_t);
+
+ template<typename iter_t>
+ void fit(iter_t begin, iter_t end);
+
+ template<typename iter_t, typename adapt_t>
+ void fit(iter_t begin, iter_t end, adapt_t adapt);
+
+ void unionAABB(const aabb<ndim> &a);
+
+ void expand(double pad);
+
+ bool completelyContains(const aabb<ndim> &other) const;
+
+ bool containsPoint(const vector_t &v) const;
+
+ bool intersectsLineSegment(const vector_t &v1, const vector_t &v2) const;
+
+ double axisSeparation(const aabb<ndim> &other, unsigned axis) const;
+
+ double maxAxisSeparation(const aabb<ndim> &other) const;
+
+ bool intersects(const aabb<ndim> &other) const;
+ bool intersects(const sphere<ndim> &s) const;
+ bool intersects(const plane<ndim> &plane) const;
+ bool intersects(const ray<ndim> &ray) const;
+ bool intersects(tri<ndim> tri) const;
+ bool intersects(const linesegment<ndim> &ls) const;
+
+ std::pair<double, double> rangeInDirection(const carve::geom::vector<ndim> &v) const;
+
+ vector_t min() const;
+ vector_t mid() const;
+ vector_t max() const;
+
+ double min(unsigned dim) const;
+ double mid(unsigned dim) const;
+ double max(unsigned dim) const;
+
+ double volume() const;
+
+ int compareAxis(const axis_pos &ap) const;
+
+ void constrainMax(const axis_pos &ap);
+ void constrainMin(const axis_pos &ap);
+
+ aabb getAABB() const;
+
+ aabb(const vector_t &_pos = vector_t::ZERO(),
+ const vector_t &_extent = vector_t::ZERO());
+
+ template<typename iter_t, typename adapt_t>
+ aabb(iter_t begin, iter_t end, adapt_t adapt);
+
+ template<typename iter_t>
+ aabb(iter_t begin, iter_t end);
+
+ aabb(const aabb<ndim> &a, const aabb<ndim> &b);
+ };
+
+ template<unsigned ndim>
+ bool operator==(const aabb<ndim> &a, const aabb<ndim> &b);
+
+ template<unsigned ndim>
+ bool operator!=(const aabb<ndim> &a, const aabb<ndim> &b);
+
+ template<unsigned ndim>
+ std::ostream &operator<<(std::ostream &o, const aabb<ndim> &a);
+
+
+
+ template<unsigned ndim, typename obj_t>
+ struct get_aabb {
+ aabb<ndim> operator()(const obj_t &obj) const {
+ return obj.getAABB();
+ }
+ };
+
+ template<unsigned ndim, typename obj_t>
+ struct get_aabb<ndim, obj_t *> {
+ aabb<ndim> operator()(const obj_t *obj) const {
+ return obj->getAABB();
+ }
+ };
+
+
+
+ }
+}
+
+namespace carve {
+ namespace geom3d {
+ typedef carve::geom::aabb<3> AABB;
+ }
+}
+
+#include <carve/aabb_impl.hpp>
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/vector.hpp>
+#include <carve/geom3d.hpp>
+
+#include <carve/geom.hpp>
+
+#include <vector>
+
+namespace carve {
+ namespace geom {
+
+ template<unsigned ndim>
+ void aabb<ndim>::empty() {
+ pos.setZero();
+ extent.setZero();
+ }
+
+ template<unsigned ndim>
+ bool aabb<ndim>::isEmpty() const {
+ return extent.exactlyZero();
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t, typename value_type>
+ void aabb<ndim>::_fit(iter_t begin, iter_t end, value_type) {
+ if (begin == end) {
+ empty();
+ return;
+ }
+
+ vector_t min, max;
+ aabb<ndim> a = get_aabb<ndim, value_type>()(*begin); ++begin;
+ min = a.min();
+ max = a.max();
+ while (begin != end) {
+ aabb<ndim> a = get_aabb<ndim, value_type>()(*begin); ++begin;
+ assign_op(min, min, a.min(), carve::util::min_functor());
+ assign_op(max, max, a.max(), carve::util::max_functor());
+ }
+
+ pos = (min + max) / 2.0;
+ assign_op(extent, max - pos, pos - min, carve::util::max_functor());
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t>
+ void aabb<ndim>::_fit(iter_t begin, iter_t end, vector_t) {
+ if (begin == end) {
+ empty();
+ return;
+ }
+
+ vector_t min, max;
+ bounds(begin, end, min, max);
+ pos = (min + max) / 2.0;
+ assign_op(extent, max - pos, pos - min, carve::util::max_functor());
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t>
+ void aabb<ndim>::_fit(iter_t begin, iter_t end, aabb_t) {
+ if (begin == end) {
+ empty();
+ return;
+ }
+
+ vector_t min, max;
+ aabb<ndim> a = *begin++;
+ min = a.min();
+ max = a.max();
+ while (begin != end) {
+ aabb<ndim> a = *begin; ++begin;
+ assign_op(min, min, a.min(), carve::util::min_functor());
+ assign_op(max, max, a.max(), carve::util::max_functor());
+ }
+
+ pos = (min + max) / 2.0;
+ assign_op(extent, max - pos, pos - min, carve::util::max_functor());
+ }
+
+ template<unsigned ndim>
+ void aabb<ndim>::fit(const vector_t &v1) {
+ pos = v1;
+ extent.setZero();
+ }
+
+ template<unsigned ndim>
+ void aabb<ndim>::fit(const vector_t &v1, const vector_t &v2) {
+ vector_t min, max;
+ assign_op(min, v1, v2, carve::util::min_functor());
+ assign_op(max, v1, v2, carve::util::max_functor());
+
+ pos = (min + max) / 2.0;
+ assign_op(extent, max - pos, pos - min, carve::util::max_functor());
+ }
+
+ template<unsigned ndim>
+ void aabb<ndim>::fit(const vector_t &v1, const vector_t &v2, const vector_t &v3) {
+ vector_t min, max;
+ min = max = v1;
+
+ assign_op(min, min, v2, carve::util::min_functor());
+ assign_op(max, max, v2, carve::util::max_functor());
+ assign_op(min, min, v3, carve::util::min_functor());
+ assign_op(max, max, v3, carve::util::max_functor());
+
+ pos = (min + max) / 2.0;
+ assign_op(extent, max - pos, pos - min, carve::util::max_functor());
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t, typename adapt_t>
+ void aabb<ndim>::fit(iter_t begin, iter_t end, adapt_t adapt) {
+ vector_t min, max;
+
+ bounds(begin, end, adapt, min, max);
+ pos = (min + max) / 2.0;
+ assign_op(extent, max - pos, pos - min, carve::util::max_functor());
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t>
+ void aabb<ndim>::fit(iter_t begin, iter_t end) {
+ _fit(begin, end, typename std::iterator_traits<iter_t>::value_type());
+ }
+
+ template<unsigned ndim>
+ void aabb<ndim>::expand(double pad) {
+ extent += pad;
+ }
+
+ template<unsigned ndim>
+ void aabb<ndim>::unionAABB(const aabb<ndim> &a) {
+ vector_t vmin, vmax;
+
+ assign_op(vmin, min(), a.min(), carve::util::min_functor());
+ assign_op(vmax, max(), a.max(), carve::util::max_functor());
+ pos = (vmin + vmax) / 2.0;
+ assign_op(extent, vmax - pos, pos - vmin, carve::util::max_functor());
+ }
+
+ template<unsigned ndim>
+ bool aabb<ndim>::completelyContains(const aabb<ndim> &other) const {
+ for (unsigned i = 0; i < ndim; ++i) {
+ if (fabs(other.pos.v[i] - pos.v[i]) + other.extent.v[i] > extent.v[i]) return false;
+ }
+ return true;
+ }
+
+ template<unsigned ndim>
+ bool aabb<ndim>::containsPoint(const vector_t &v) const {
+ for (unsigned i = 0; i < ndim; ++i) {
+ if (fabs(v.v[i] - pos.v[i]) > extent.v[i]) return false;
+ }
+ return true;
+ }
+
+ template<unsigned ndim>
+ double aabb<ndim>::axisSeparation(const aabb<ndim> &other, unsigned axis) const {
+ return fabs(other.pos.v[axis] - pos.v[axis]) - extent.v[axis] - other.extent.v[axis];
+ }
+
+ template<unsigned ndim>
+ double aabb<ndim>::maxAxisSeparation(const aabb<ndim> &other) const {
+ double m = axisSeparation(other, 0);
+ for (unsigned i = 1; i < ndim; ++i) {
+ m = std::max(m, axisSeparation(other, i));
+ }
+ return m;
+ }
+
+ template<unsigned ndim>
+ bool aabb<ndim>::intersects(const aabb<ndim> &other) const {
+ return maxAxisSeparation(other) <= 0.0;
+ }
+
+ template<unsigned ndim>
+ bool aabb<ndim>::intersects(const sphere<ndim> &s) const {
+ double r = 0.0;
+ for (unsigned i = 0; i < ndim; ++i) {
+ double t = fabs(s.C[i] - pos[i]) - extent[i]; if (t > 0.0) r += t*t;
+ }
+ return r <= s.r*s.r;
+ }
+
+ template<unsigned ndim>
+ bool aabb<ndim>::intersects(const plane<ndim> &plane) const {
+ double d1 = fabs(distance(plane, pos));
+ double d2 = dot(abs(plane.N), extent);
+ return d1 <= d2;
+ }
+
+ template<unsigned ndim>
+ bool aabb<ndim>::intersects(const linesegment<ndim> &ls) const {
+ return intersectsLineSegment(ls.v1, ls.v2);
+ }
+
+ template<unsigned ndim>
+ std::pair<double, double> aabb<ndim>::rangeInDirection(const carve::geom::vector<ndim> &v) const {
+ double d1 = dot(v, pos);
+ double d2 = dot(abs(v), extent);
+
+ return std::make_pair(d1 - d2, d1 + d2);
+ }
+
+ template<unsigned ndim>
+ typename aabb<ndim>::vector_t aabb<ndim>::min() const { return pos - extent; }
+
+ template<unsigned ndim>
+ typename aabb<ndim>::vector_t aabb<ndim>::mid() const { return pos; }
+
+ template<unsigned ndim>
+ typename aabb<ndim>::vector_t aabb<ndim>::max() const { return pos + extent; }
+
+ template<unsigned ndim>
+ double aabb<ndim>::min(unsigned dim) const { return pos.v[dim] - extent.v[dim]; }
+
+ template<unsigned ndim>
+ double aabb<ndim>::mid(unsigned dim) const { return pos.v[dim]; }
+
+ template<unsigned ndim>
+ double aabb<ndim>::max(unsigned dim) const { return pos.v[dim] + extent.v[dim]; }
+
+ template<unsigned ndim>
+ double aabb<ndim>::volume() const {
+ double v = 1.0;
+ for (size_t dim = 0; dim < ndim; ++dim) { v *= 2.0 * extent.v[dim]; }
+ return v;
+ }
+
+ template<unsigned ndim>
+ int aabb<ndim>::compareAxis(const axis_pos &ap) const {
+ double p = ap.pos - pos[ap.axis];
+ if (p > extent[ap.axis]) return -1;
+ if (p < -extent[ap.axis]) return +1;
+ return 0;
+ }
+
+ template<unsigned ndim>
+ void aabb<ndim>::constrainMax(const axis_pos &ap) {
+ if (pos[ap.axis] + extent[ap.axis] > ap.pos) {
+ double min = std::min(ap.pos, pos[ap.axis] - extent[ap.axis]);
+ pos[ap.axis] = (min + ap.pos) / 2.0;
+ extent[ap.axis] = ap.pos - pos[ap.axis];
+ }
+ }
+
+ template<unsigned ndim>
+ void aabb<ndim>::constrainMin(const axis_pos &ap) {
+ if (pos[ap.axis] - extent[ap.axis] < ap.pos) {
+ double max = std::max(ap.pos, pos[ap.axis] + extent[ap.axis]);
+ pos[ap.axis] = (ap.pos + max) / 2.0;
+ extent[ap.axis] = pos[ap.axis] - ap.pos;
+ }
+ }
+
+ template<unsigned ndim>
+ aabb<ndim> aabb<ndim>::getAABB() const {
+ return *this;
+ }
+
+ template<unsigned ndim>
+ aabb<ndim>::aabb(const vector_t &_pos,
+ const vector_t &_extent) : pos(_pos), extent(_extent) {
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t, typename adapt_t>
+ aabb<ndim>::aabb(iter_t begin, iter_t end, adapt_t adapt) {
+ fit(begin, end, adapt);
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t>
+ aabb<ndim>::aabb(iter_t begin, iter_t end) {
+ fit(begin, end);
+ }
+
+ template<unsigned ndim>
+ aabb<ndim>::aabb(const aabb<ndim> &a, const aabb<ndim> &b) {
+ fit(a, b);
+ }
+
+ template<unsigned ndim>
+ bool operator==(const aabb<ndim> &a, const aabb<ndim> &b) {
+ return a.pos == b.pos && a.extent == b.extent;
+ }
+
+ template<unsigned ndim>
+ bool operator!=(const aabb<ndim> &a, const aabb<ndim> &b) {
+ return a.pos != b.pos || a.extent != b.extent;
+ }
+
+ template<unsigned ndim>
+ std::ostream &operator<<(std::ostream &o, const aabb<ndim> &a) {
+ o << (a.pos - a.extent) << "--" << (a.pos + a.extent);
+ return o;
+ }
+
+ template<>
+ inline bool aabb<3>::intersects(const ray<3> &ray) const {
+ vector<3> t = pos - ray.v;
+ double r;
+
+ //l.cross(x-axis)?
+ r = extent.y * fabs(ray.D.z) + extent.z * fabs(ray.D.y);
+ if (fabs(t.y * ray.D.z - t.z * ray.D.y) > r) return false;
+
+ //ray.D.cross(y-axis)?
+ r = extent.x * fabs(ray.D.z) + extent.z * fabs(ray.D.x);
+ if (fabs(t.z * ray.D.x - t.x * ray.D.z) > r) return false;
+
+ //ray.D.cross(z-axis)?
+ r = extent.x*fabs(ray.D.y) + extent.y*fabs(ray.D.x);
+ if (fabs(t.x * ray.D.y - t.y * ray.D.x) > r) return false;
+
+ return true;
+ }
+
+ template<>
+ inline bool aabb<3>::intersectsLineSegment(const vector<3> &v1, const vector<3> &v2) const {
+ vector<3> half_length = 0.5 * (v2 - v1);
+ vector<3> t = pos - half_length - v1;
+ double r;
+
+ //do any of the principal axes form a separating axis?
+ if(fabs(t.x) > extent.x + fabs(half_length.x)) return false;
+ if(fabs(t.y) > extent.y + fabs(half_length.y)) return false;
+ if(fabs(t.z) > extent.z + fabs(half_length.z)) return false;
+
+ // NOTE: Since the separating axis is perpendicular to the line in
+ // these last four cases, the line does not contribute to the
+ // projection.
+
+ //line.cross(x-axis)?
+ r = extent.y * fabs(half_length.z) + extent.z * fabs(half_length.y);
+ if (fabs(t.y * half_length.z - t.z * half_length.y) > r) return false;
+
+ //half_length.cross(y-axis)?
+ r = extent.x * fabs(half_length.z) + extent.z * fabs(half_length.x);
+ if (fabs(t.z * half_length.x - t.x * half_length.z) > r) return false;
+
+ //half_length.cross(z-axis)?
+ r = extent.x*fabs(half_length.y) + extent.y*fabs(half_length.x);
+ if (fabs(t.x * half_length.y - t.y * half_length.x) > r) return false;
+
+ return true;
+ }
+
+ template<int Ax, int Ay, int Az, int c>
+ static inline bool intersectsTriangle_axisTest_3(const aabb<3> &aabb, const tri<3> &tri) {
+ const int d = (c+1) % 3, e = (c+2) % 3;
+ const vector<3> a = cross(VECTOR(Ax, Ay, Az), tri.v[d] - tri.v[c]);
+ double p1 = dot(a, tri.v[c]), p2 = dot(a, tri.v[e]);
+ if (p1 > p2) std::swap(p1, p2);
+ const double r = dot(abs(a), aabb.extent);
+ return !(p1 > r || p2 < -r);
+ }
+
+ template<int c>
+ static inline bool intersectsTriangle_axisTest_2(const aabb<3> &aabb, const tri<3> &tri) {
+ double vmin = std::min(std::min(tri.v[0][c], tri.v[1][c]), tri.v[2][c]),
+ vmax = std::max(std::max(tri.v[0][c], tri.v[1][c]), tri.v[2][c]);
+ return !(vmin > aabb.extent[c] || vmax < -aabb.extent[c]);
+ }
+
+ static inline bool intersectsTriangle_axisTest_1(const aabb<3> &aabb, const tri<3> &tri) {
+ vector<3> n = cross(tri.v[1] - tri.v[0], tri.v[2] - tri.v[0]);
+ double d1 = fabs(dot(n, tri.v[0]));
+ double d2 = dot(abs(n), aabb.extent);
+ return d1 <= d2;
+ }
+
+ template<>
+ inline bool aabb<3>::intersects(tri<3> tri) const {
+ tri.v[0] -= pos;
+ tri.v[1] -= pos;
+ tri.v[2] -= pos;
+
+ if (!intersectsTriangle_axisTest_2<0>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_2<1>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_2<2>(*this, tri)) return false;
+
+ if (!intersectsTriangle_axisTest_3<1,0,0,0>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_3<1,0,0,1>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_3<1,0,0,2>(*this, tri)) return false;
+
+ if (!intersectsTriangle_axisTest_3<0,1,0,0>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_3<0,1,0,1>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_3<0,1,0,2>(*this, tri)) return false;
+
+ if (!intersectsTriangle_axisTest_3<0,0,1,0>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_3<0,0,1,1>(*this, tri)) return false;
+ if (!intersectsTriangle_axisTest_3<0,0,1,2>(*this, tri)) return false;
+
+ if (!intersectsTriangle_axisTest_1(*this, tri)) return false;
+
+ return true;
+ }
+
+
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+#pragma once
+
+#if defined(CMAKE_BUILD)
+# include <carve/config.h>
+#elif defined(XCODE_BUILD)
+# include <carve/xcode_config.h>
+#elif defined(_MSC_VER)
+# include <carve/vcpp_config.h>
+#else
+# include <carve/config.h>
+#endif
+
+#if defined(WIN32)
+# include <carve/win32.h>
+#elif defined(__GNUC__)
+# include <carve/gnu_cxx.h>
+#endif
+
+#if defined(CARVE_SYSTEM_BOOST)
+# define BOOST_INCLUDE(x) <boost/x>
+#else
+# define BOOST_INCLUDE(x) <carve/external/boost/x>
+#endif
+
+#include <math.h>
+
+#include <string>
+#include <set>
+#include <map>
+#include <vector>
+#include <list>
+#include <sstream>
+#include <iomanip>
+
+#include <carve/collection.hpp>
+
+#include <carve/util.hpp>
+
+#include <stdarg.h>
+
+#define STR(x) #x
+#define XSTR(x) STR(x)
+
+/**
+ * \brief Top level Carve namespace.
+ */
+namespace carve {
+ static struct noinit_t {} NOINIT;
+
+ inline std::string fmtstring(const char *fmt, ...);
+
+ /**
+ * \brief Base class for all Carve exceptions.
+ */
+ struct exception {
+ private:
+ mutable std::string err;
+ mutable std::ostringstream accum;
+
+ public:
+ exception(const std::string &e) : err(e), accum() { }
+ exception() : err(), accum() { }
+ exception(const exception &e) : err(e.str()), accum() { }
+ exception &operator=(const exception &e) {
+ if (this != &e) {
+ err = e.str();
+ accum.str("");
+ }
+ return *this;
+ }
+
+ const std::string &str() const {
+ if (accum.str().size() > 0) {
+ err += accum.str();
+ accum.str("");
+ }
+ return err;
+ }
+
+ template<typename T>
+ exception &operator<<(const T &t) {
+ accum << t;
+ return *this;
+ }
+ };
+
+ template<typename iter_t, typename order_t = std::less<typename std::iterator_traits<iter_t>::value_type > >
+ struct index_sort {
+ iter_t base;
+ order_t order;
+ index_sort(const iter_t &_base) : base(_base), order() { }
+ index_sort(const iter_t &_base, const order_t &_order) : base(_base), order(_order) { }
+ template<typename U>
+ bool operator()(const U &a, const U &b) {
+ return order(*(base + a), *(base + b));
+ }
+ };
+
+ template<typename iter_t, typename order_t>
+ index_sort<iter_t, order_t> make_index_sort(const iter_t &base, const order_t &order) {
+ return index_sort<iter_t, order_t>(base, order);
+ }
+
+ template<typename iter_t>
+ index_sort<iter_t> make_index_sort(const iter_t &base) {
+ return index_sort<iter_t>(base);
+ }
+
+
+ enum RayIntersectionClass {
+ RR_DEGENERATE = -2,
+ RR_PARALLEL = -1,
+ RR_NO_INTERSECTION = 0,
+ RR_INTERSECTION = 1
+ };
+
+ enum LineIntersectionClass {
+ COLINEAR = -1,
+ NO_INTERSECTION = 0,
+ INTERSECTION_LL = 1,
+ INTERSECTION_PL = 2,
+ INTERSECTION_LP = 3,
+ INTERSECTION_PP = 4
+ };
+
+ enum PointClass {
+ POINT_UNK = -2,
+ POINT_OUT = -1,
+ POINT_ON = 0,
+ POINT_IN = 1,
+ POINT_VERTEX = 2,
+ POINT_EDGE = 3
+ };
+
+ enum IntersectionClass {
+ INTERSECT_BAD = -1,
+ INTERSECT_NONE = 0,
+ INTERSECT_FACE = 1,
+ INTERSECT_VERTEX = 2,
+ INTERSECT_EDGE = 3,
+ INTERSECT_PLANE = 4,
+ };
+
+
+
+ extern double EPSILON;
+ extern double EPSILON2;
+
+ static inline void setEpsilon(double ep) { EPSILON = ep; EPSILON2 = ep * ep; }
+
+
+
+ template<typename T>
+ struct identity_t {
+ typedef T argument_type;
+ typedef T result_type;
+ const T &operator()(const T &t) const { return t; }
+ };
+
+
+
+ template<typename iter_t>
+ inline bool is_sorted(iter_t first, iter_t last) {
+ if (first == last) return true;
+
+ iter_t iter = first;
+ iter_t next = first; ++next;
+ for (; next != last; iter = next, ++next) {
+ if (*next < *iter) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+
+
+ template<typename iter_t,
+ typename pred_t>
+ inline bool is_sorted(iter_t first, iter_t last, pred_t pred) {
+ if (first == last) return true;
+
+ iter_t iter = first;
+ iter_t next = first; ++next;
+ for (; next != last; iter = next, ++next) {
+ if (pred(*next, *iter)) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+
+
+ inline double rangeSeparation(const std::pair<double, double> &a,
+ const std::pair<double, double> &b) {
+ if (a.second < b.first) {
+ return b.first - a.second;
+ } else if (b.second < a.first) {
+ return a.first - b.second;
+ } else {
+ return 0.0;
+ }
+ }
+}
+
+
+#if defined(_MSC_VER)
+# define MACRO_BEGIN do {
+# define MACRO_END __pragma(warning(push)) __pragma(warning(disable:4127)) } while(0) __pragma(warning(pop))
+#else
+# define MACRO_BEGIN do {
+# define MACRO_END } while(0)
+#endif
+
+#if !defined(CARVE_NODEBUG)
+# define CARVE_ASSERT(x) MACRO_BEGIN if (!(x)) throw carve::exception() << __FILE__ << ":" << __LINE__ << " " << #x; MACRO_END
+#else
+# define CARVE_ASSERT(X)
+#endif
+
+#define CARVE_FAIL(x) MACRO_BEGIN throw carve::exception() << __FILE__ << ":" << __LINE__ << " " << #x; MACRO_END
--- /dev/null
+// N.B. only appropriate for IEEE doubles.
+// Cube root implementation obtained from code with the following notice:
+
+/* @(#)s_cbrt.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+/* Sometimes it's necessary to define __LITTLE_ENDIAN explicitly
+ but these catch some common cases. */
+
+#if defined(i386) || defined(i486) || \
+ defined(intel) || defined(x86) || defined(i86pc) || \
+ defined(__alpha) || defined(__osf__)
+#define __LITTLE_ENDIAN
+#endif
+
+#ifdef __LITTLE_ENDIAN
+#define __HI(x) *(1+(int*)&x)
+#define __LO(x) *(int*)&x
+#define __HIp(x) *(1+(int*)x)
+#define __LOp(x) *(int*)x
+#else
+#define __HI(x) *(int*)&x
+#define __LO(x) *(1+(int*)&x)
+#define __HIp(x) *(int*)x
+#define __LOp(x) *(1+(int*)x)
+#endif
+
+/* cbrt(x)
+ * Return cube root of x
+ */
+
+inline double cbrt(double x) {
+
+ static const unsigned
+ B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
+ B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
+ static const double
+ C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */
+ D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */
+ E = 1.41428571428571436819e+00, /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */
+ F = 1.60714285714285720630e+00, /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */
+ G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
+
+ int hx;
+ double r,s,t=0.0,w;
+ unsigned sign;
+
+ hx = __HI(x); /* high word of x */
+ sign=hx&0x80000000; /* sign= sign(x) */
+ hx ^=sign;
+ if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
+ if((hx|__LO(x))==0)
+ return(x); /* cbrt(0) is itself */
+
+ __HI(x) = hx; /* x <- |x| */
+ /* rough cbrt to 5 bits */
+ if(hx<0x00100000) /* subnormal number */
+ {__HI(t)=0x43500000; /* set t= 2**54 */
+ t*=x; __HI(t)=__HI(t)/3+B2;
+ }
+ else
+ __HI(t)=hx/3+B1;
+
+
+ /* new cbrt to 23 bits, may be implemented in single precision */
+ r=t*t/x;
+ s=C+r*t;
+ t*=G+F/(s+E+D/s);
+
+ /* chopped to 20 bits and make it larger than cbrt(x) */
+ __LO(t)=0; __HI(t)+=0x00000001;
+
+ /* one step newton iteration to 53 bits with error less than 0.667 ulps */
+ s=t*t; /* t*t is exact */
+ r=x/s;
+ w=t+t;
+ r=(r-t)/(w+r); /* r-s is exact */
+ t=t+t*r;
+
+ /* retore the sign bit */
+ __HI(t) |= sign;
+ return(t);
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/collection_types.hpp>
+
+namespace carve {
+ namespace csg {
+
+ enum FaceClass {
+ FACE_UNCLASSIFIED = -3,
+ FACE_ON_ORIENT_OUT = -2,
+ FACE_OUT = -1,
+ FACE_ON = 0,
+ FACE_IN = +1,
+ FACE_ON_ORIENT_IN = +2
+ };
+
+ enum FaceClassBit {
+ FACE_ON_ORIENT_OUT_BIT = 0x01,
+ FACE_OUT_BIT = 0x02,
+ FACE_IN_BIT = 0x04,
+ FACE_ON_ORIENT_IN_BIT = 0x08,
+
+ FACE_ANY_BIT = 0x0f,
+ FACE_ON_BIT = 0x09,
+ FACE_NOT_ON_BIT = 0x06
+ };
+
+ static inline FaceClass class_bit_to_class(unsigned i) {
+ if (i & FACE_ON_ORIENT_OUT_BIT) return FACE_ON_ORIENT_OUT;
+ if (i & FACE_OUT_BIT) return FACE_OUT;
+ if (i & FACE_IN_BIT) return FACE_IN;
+ if (i & FACE_ON_ORIENT_IN_BIT) return FACE_ON_ORIENT_IN;
+ return FACE_UNCLASSIFIED;
+ }
+
+ static inline unsigned class_to_class_bit(FaceClass f) {
+ switch (f) {
+ case FACE_ON_ORIENT_OUT: return FACE_ON_ORIENT_OUT_BIT;
+ case FACE_OUT: return FACE_OUT_BIT;
+ case FACE_ON: return FACE_ON_BIT;
+ case FACE_IN: return FACE_IN_BIT;
+ case FACE_ON_ORIENT_IN: return FACE_ON_ORIENT_IN_BIT;
+ case FACE_UNCLASSIFIED: return FACE_ANY_BIT;
+ default: return 0;
+ }
+ }
+
+ enum EdgeClass {
+ EDGE_UNK = -2,
+ EDGE_OUT = -1,
+ EDGE_ON = 0,
+ EDGE_IN = 1
+ };
+
+
+
+ const char *ENUM(FaceClass f);
+ const char *ENUM(PointClass p);
+
+
+
+ struct ClassificationInfo {
+ const carve::mesh::Mesh<3> *intersected_mesh;
+ FaceClass classification;
+
+ ClassificationInfo() : intersected_mesh(NULL), classification(FACE_UNCLASSIFIED) { }
+ ClassificationInfo(const carve::mesh::Mesh<3> *_intersected_mesh,
+ FaceClass _classification) :
+ intersected_mesh(_intersected_mesh),
+ classification(_classification) {
+ }
+ bool intersectedMeshIsClosed() const {
+ return intersected_mesh->isClosed();
+ }
+ };
+
+
+
+ struct EC2 {
+ EdgeClass cls[2];
+ EC2() { cls[0] = cls[1] = EDGE_UNK; }
+ EC2(EdgeClass a, EdgeClass b) { cls[0] = a; cls[1] = b; }
+ };
+
+ struct PC2 {
+ PointClass cls[2];
+ PC2() { cls[0] = cls[1] = POINT_UNK; }
+ PC2(PointClass a, PointClass b) { cls[0] = a; cls[1] = b; }
+ };
+
+ typedef std::unordered_map<std::pair<const carve::mesh::MeshSet<3>::vertex_t *, const carve::mesh::MeshSet<3>::vertex_t *>,
+ EC2> EdgeClassification;
+
+ typedef std::unordered_map<const carve::mesh::Vertex<3> *, PC2> VertexClassification;
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/collection/unordered.hpp>
+
+namespace carve {
+
+ template<typename set_t>
+ class set_insert_iterator : public std::iterator<std::output_iterator_tag, void, void, void, void> {
+
+ protected:
+ set_t *set;
+ public:
+
+ set_insert_iterator(set_t &s) : set(&s) {
+ }
+
+ set_insert_iterator &
+ operator=(typename set_t::const_reference value) {
+ set->insert(value);
+ return *this;
+ }
+
+ set_insert_iterator &operator*() { return *this; }
+ set_insert_iterator &operator++() { return *this; }
+ set_insert_iterator &operator++(int) { return *this; }
+ };
+
+ template<typename set_t>
+ inline set_insert_iterator<set_t>
+ set_inserter(set_t &s) {
+ return set_insert_iterator<set_t>(s);
+ }
+
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+#pragma once
+
+#if defined(HAVE_STD_UNORDERED_COLLECTIONS)
+
+# include <carve/collection/unordered/std_impl.hpp>
+
+#elif defined(HAVE_TR1_UNORDERED_COLLECTIONS)
+
+# include <carve/collection/unordered/tr1_impl.hpp>
+
+#elif defined(HAVE_BOOST_UNORDERED_COLLECTIONS)
+
+# include <carve/collection/unordered/boost_impl.hpp>
+
+#elif defined(HAVE_LIBSTDCPP_UNORDERED_COLLECTIONS)
+
+# include <carve/collection/unordered/libstdcpp_impl.hpp>
+
+#elif defined(_MSC_VER) && _MSC_VER >= 1300
+
+# include <carve/collection/unordered/vcpp_impl.hpp>
+
+#else
+
+# include <carve/collection/unordered/fallback_impl.hpp>
+
+#endif
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include BOOST_INCLUDE(unordered_set.hpp)
+#include BOOST_INCLUDE(unordered_map.hpp)
+
+#include <functional>
+
+namespace std {
+ template <typename Key, typename T, typename Hash = boost::hash<Key>,
+ typename Pred = std::equal_to<Key> >
+ class unordered_map : public boost::unordered_map<Key, T, Hash, Pred> {
+
+ public:
+ typedef T data_type;
+ };
+
+ template <typename Key, typename T, typename Hash = boost::hash<Key>,
+ typename Pred = std::equal_to<Key> >
+ class unordered_multimap : public boost::unordered_multimap<Key, T, Hash, Pred> {
+ };
+
+ template <typename Value, typename Hash = boost::hash<Value>,
+ typename Pred = std::equal_to<Value> >
+ class unordered_set : public boost::unordered_set<Value, Hash, Pred> {
+ };
+}
+
+#undef UNORDERED_COLLECTIONS_SUPPORT_RESIZE
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <set>
+#include <map>
+
+namespace std {
+
+ template<typename K, typename T, typename H = int>
+ class unordered_map : public std::map<K, T> {
+ typedef std::map<K, T> super;
+ public:
+ typedef T data_type;
+ };
+
+ template<typename K, typename H = int>
+ class unordered_set : public std::set<K> {
+ typedef std::set<K> super;
+ public:
+ };
+
+}
+
+#undef UNORDERED_COLLECTIONS_SUPPORT_RESIZE
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+#pragma once
+
+#include <ext/hash_map>
+#include <ext/hash_set>
+
+namespace __gnu_cxx {
+ template <typename T>
+ struct hash<T *> : public std::unary_function<T *, size_t> {
+ size_t operator()(T *v) const {
+ size_t x = (size_t)(v);
+ return x + (x>>3);
+ }
+ };
+
+ template <typename A, typename B>
+ struct hash<std::pair<A, B> > : public std::unary_function<std::pair<A, B>, size_t> {
+ size_t operator()(const std::pair<A, B> &v) const {
+ std::size_t seed = 0;
+
+ seed ^= hash<A>()(v.first);
+ seed ^= hash<B>()(v.second) + (seed<<6) + (seed>>2);
+
+ return seed;
+ }
+ };
+}
+
+namespace std {
+
+ template<typename K, typename V, typename H = __gnu_cxx::hash<K> >
+ class unordered_map : public __gnu_cxx::hash_map<K, V, H> {
+ typedef __gnu_cxx::hash_map<K, V, H> super;
+ public:
+ typedef typename super::mapped_type data_type;
+ };
+
+ template<typename K, typename H = __gnu_cxx::hash<K> >
+ class unordered_set : public __gnu_cxx::hash_set<K, H> {
+ typedef __gnu_cxx::hash_set<K, H> super;
+ public:
+ };
+
+}
+
+#define UNORDERED_COLLECTIONS_SUPPORT_RESIZE 1
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <unordered_map>
+#include <unordered_set>
+
+#undef UNORDERED_COLLECTIONS_SUPPORT_RESIZE
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <tr1/unordered_map>
+#include <tr1/unordered_set>
+#include <tr1/functional>
+
+namespace std {
+ namespace tr1 {
+ template <typename A, typename B>
+ struct hash<std::pair<A, B> > : public std::unary_function<std::pair<A, B>, size_t> {
+ size_t operator()(const std::pair<A, B> &v) const {
+ std::size_t seed = 0;
+
+ seed ^= hash<A>()(v.first);
+ seed ^= hash<B>()(v.second) + (seed<<6) + (seed>>2);
+
+ return seed;
+ }
+ };
+ }
+
+
+
+ template <typename Key, typename T,
+ typename Hash = tr1::hash<Key>,
+ typename Pred = std::equal_to<Key> >
+ class unordered_map : public std::tr1::unordered_map<Key, T, Hash, Pred> {
+ public:
+ typedef T data_type;
+ };
+
+ template <typename Value,
+ typename Hash = tr1::hash<Value>,
+ typename Pred = std::equal_to<Value> >
+ class unordered_set : public std::tr1::unordered_set<Value, Hash, Pred> {
+ public:
+ };
+
+}
+
+#undef UNORDERED_COLLECTIONS_SUPPORT_RESIZE
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <hash_map>
+#include <hash_set>
+
+namespace std {
+
+ namespace {
+
+ template<class Value, class Hash> class hash_traits {
+ Hash hash_value;
+ std::less<Value> comp;
+ public:
+ enum {
+ bucket_size = 4,
+ min_buckets = 8
+ };
+ // hash _Keyval to size_t value
+ size_t operator()(const Value& v) const {
+ return ((size_t)hash_value(v));
+ }
+ // test if _Keyval1 ordered before _Keyval2
+ bool operator()(const Value& v1, const Value& v2) const {
+ return (comp(v1, v2));
+ }
+ };
+
+ }
+
+ template <typename Key, typename T, typename Hash = stdext::hash_compare<Key, less<Key> >, typename Pred = std::equal_to<Key> >
+ class unordered_map
+ : public stdext::hash_map<Key, T, hash_traits<Key, Hash> > {
+ typedef stdext::hash_map<Key, T, hash_traits<Key, Hash> > super;
+ public:
+ unordered_map() : super() {}
+ };
+
+ template <typename Value, typename Hash = stdext::hash_compare<Key, less<Key> >, typename Pred = std::equal_to<Value> >
+ class unordered_set
+ : public stdext::hash_set<Value, hash_traits<Value, Hash> > {
+ typedef stdext::hash_set<Value, hash_traits<Value, Hash> > super;
+ public:
+ unordered_set() : super() {}
+ };
+
+}
+
+#undef UNORDERED_COLLECTIONS_SUPPORT_RESIZE
--- /dev/null
+
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/mesh.hpp>
+
+namespace carve {
+ namespace csg {
+
+ typedef std::pair<
+ carve::mesh::MeshSet<3>::vertex_t *,
+ carve::mesh::MeshSet<3>::vertex_t *> V2;
+
+ typedef std::pair<
+ carve::mesh::MeshSet<3>::face_t *,
+ carve::mesh::MeshSet<3>::face_t *> F2;
+
+ static inline V2 ordered_edge(
+ carve::mesh::MeshSet<3>::vertex_t *a,
+ carve::mesh::MeshSet<3>::vertex_t *b) {
+ return V2(std::min(a, b), std::max(a, b));
+ }
+
+ static inline V2 flip(const V2 &v) {
+ return V2(v.second, v.first);
+ }
+
+ // include/carve/csg.hpp include/carve/faceloop.hpp
+ // lib/intersect.cpp lib/intersect_classify_common_impl.hpp
+ // lib/intersect_classify_edge.cpp
+ // lib/intersect_classify_group.cpp
+ // lib/intersect_classify_simple.cpp
+ // lib/intersect_face_division.cpp lib/intersect_group.cpp
+ // lib/intersect_half_classify_group.cpp
+ typedef std::unordered_set<V2> V2Set;
+
+ // include/carve/csg.hpp include/carve/polyhedron_decl.hpp
+ // lib/csg_collector.cpp lib/intersect.cpp
+ // lib/intersect_common.hpp lib/intersect_face_division.cpp
+ // lib/polyhedron.cpp
+ typedef std::unordered_map<
+ carve::mesh::MeshSet<3>::vertex_t *,
+ carve::mesh::MeshSet<3>::vertex_t *> VVMap;
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/geom.hpp>
+
+namespace carve {
+ namespace colour {
+ static inline void HSV2RGB(float H, float S, float V, float &r, float &g, float &b) {
+ H = 6.0f * H;
+ if (S < 5.0e-6) {
+ r = g = b = V; return;
+ } else {
+ int i = (int)H;
+ float f = H - i;
+ float p1 = V * (1.0f - S);
+ float p2 = V * (1.0f - S * f);
+ float p3 = V * (1.0f - S * (1.0f - f));
+ switch (i) {
+ case 0: r = V; g = p3; b = p1; return;
+ case 1: r = p2; g = V; b = p1; return;
+ case 2: r = p1; g = V; b = p3; return;
+ case 3: r = p1; g = p2; b = V; return;
+ case 4: r = p3; g = p1; b = V; return;
+ case 5: r = V; g = p1; b = p2; return;
+ }
+ }
+ r = g = b = 0.0;
+ }
+ }
+}
--- /dev/null
+#define CARVE_VERSION "2.0.0a"
+
+#undef CARVE_DEBUG
+#undef CARVE_DEBUG_WRITE_PLY_DATA
+
+#if defined(__GNUC__)
+# if !defined(HAVE_BOOST_UNORDERED_COLLECTIONS)
+# define HAVE_TR1_UNORDERED_COLLECTIONS
+# endif
+
+# define HAVE_STDINT_H
+#endif
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <list>
+#include <vector>
+#include <algorithm>
+
+#include <carve/carve.hpp>
+
+#include <carve/geom2d.hpp>
+
+namespace carve {
+ namespace geom {
+ std::vector<int> convexHull(const std::vector<carve::geom2d::P2> &points);
+
+ template<typename project_t, typename polygon_container_t>
+ std::vector<int> convexHull(const project_t &project, const polygon_container_t &points) {
+ std::vector<carve::geom2d::P2> proj;
+ proj.reserve(points.size());
+ for (typename polygon_container_t::const_iterator i = points.begin(); i != points.end(); ++i) {
+ proj.push_back(project(*i));
+ }
+ return convexHull(proj);
+ }
+
+ template<typename project_t, typename iter_t>
+ std::vector<int> convexHull(const project_t &project, iter_t beg, iter_t end, size_t size_hint = 0) {
+ std::vector<carve::geom2d::P2> proj;
+ if (size_hint) proj.reserve(size_hint);
+ for (; beg != end; ++beg) {
+ proj.push_back(project(*beg));
+ }
+ return convexHull(proj);
+ }
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <list>
+#include <vector>
+#include <algorithm>
+
+#include <carve/carve.hpp>
+
+#include <carve/geom3d.hpp>
+
+#include <carve/mesh.hpp>
+
+#include <carve/collection_types.hpp>
+#include <carve/classification.hpp>
+#include <carve/iobj.hpp>
+#include <carve/faceloop.hpp>
+#include <carve/intersection.hpp>
+#include <carve/rtree.hpp>
+
+namespace carve {
+ namespace csg {
+
+ class VertexPool {
+ typedef carve::mesh::MeshSet<3>::vertex_t vertex_t;
+
+ const static unsigned blocksize = 1024;
+ typedef std::list<std::vector<vertex_t> > pool_t;
+ pool_t pool;
+ public:
+ void reset();
+ vertex_t *get(const vertex_t::vector_t &v = vertex_t::vector_t::ZERO());
+ bool inPool(vertex_t *v) const;
+
+ VertexPool();
+ ~VertexPool();
+ };
+
+
+
+ namespace detail {
+ struct Data;
+ class LoopEdges;
+ }
+
+ /**
+ * \class CSG
+ * \brief The class responsible for the computation of CSG operations.
+ *
+ */
+ class CSG {
+ private:
+
+ public:
+ typedef carve::mesh::MeshSet<3> meshset_t;
+
+ struct Hook {
+ /**
+ * \class Hook
+ * \brief Provides API access to intermediate steps in CSG calculation.
+ *
+ */
+ virtual void intersectionVertex(const meshset_t::vertex_t * /* vertex */,
+ const IObjPairSet & /* intersections */) {
+ }
+ virtual void processOutputFace(std::vector<meshset_t::face_t *> & /* faces */,
+ const meshset_t::face_t * /* orig_face */,
+ bool /* flipped */) {
+ }
+ virtual void resultFace(const meshset_t::face_t * /* new_face */,
+ const meshset_t::face_t * /* orig_face */,
+ bool /* flipped */) {
+ }
+
+ virtual ~Hook() {
+ }
+ };
+
+ /**
+ * \class Hooks
+ * \brief Management of API hooks.
+ *
+ */
+ class Hooks {
+ public:
+ enum {
+ RESULT_FACE_HOOK = 0,
+ PROCESS_OUTPUT_FACE_HOOK = 1,
+ INTERSECTION_VERTEX_HOOK = 2,
+ HOOK_MAX = 3,
+
+ RESULT_FACE_BIT = 0x0001,
+ PROCESS_OUTPUT_FACE_BIT = 0x0002,
+ INTERSECTION_VERTEX_BIT = 0x0004
+ };
+
+ std::vector<std::list<Hook *> > hooks;
+
+ bool hasHook(unsigned hook_num);
+
+ void intersectionVertex(const meshset_t::vertex_t *vertex,
+ const IObjPairSet &intersections);
+
+ void processOutputFace(std::vector<meshset_t::face_t *> &faces,
+ const meshset_t::face_t *orig_face,
+ bool flipped);
+
+ void resultFace(const meshset_t::face_t *new_face,
+ const meshset_t::face_t *orig_face,
+ bool flipped);
+
+ void registerHook(Hook *hook, unsigned hook_bits);
+ void unregisterHook(Hook *hook);
+
+ void reset();
+
+ Hooks();
+ ~Hooks();
+ };
+
+ /**
+ * \class Collector
+ * \brief Base class for objects responsible for selecting result from which form the result polyhedron.
+ *
+ */
+ class Collector {
+ Collector(const Collector &);
+ Collector &operator=(const Collector &);
+
+ protected:
+
+ public:
+ virtual void collect(FaceLoopGroup *group, CSG::Hooks &) =0;
+ virtual meshset_t *done(CSG::Hooks &) =0;
+
+ Collector() {}
+ virtual ~Collector() {}
+ };
+
+ private:
+ typedef carve::geom::RTreeNode<3, carve::mesh::Face<3> *> face_rtree_t;
+ typedef std::unordered_map<carve::mesh::Face<3> *, std::vector<carve::mesh::Face<3> *> > face_pairs_t;
+
+ /// The computed intersection data.
+ Intersections intersections;
+
+ /// A map from intersection point to a set of intersections
+ /// represented by pairs of intersection objects.
+ VertexIntersections vertex_intersections;
+
+ /// A pool from which temporary vertices are allocated. Also
+ /// provides testing for pool membership.
+ VertexPool vertex_pool;
+
+ void init();
+
+ void makeVertexIntersections();
+
+ void groupIntersections();
+
+ void _generateVertexVertexIntersections(carve::mesh::MeshSet<3>::vertex_t *va,
+ carve::mesh::MeshSet<3>::edge_t *eb);
+ void generateVertexVertexIntersections(carve::mesh::MeshSet<3>::face_t *a,
+ const std::vector<carve::mesh::MeshSet<3>::face_t *> &b);
+
+ void _generateVertexEdgeIntersections(carve::mesh::MeshSet<3>::vertex_t *va,
+ carve::mesh::MeshSet<3>::edge_t *eb);
+ void generateVertexEdgeIntersections(carve::mesh::MeshSet<3>::face_t *a,
+ const std::vector<carve::mesh::MeshSet<3>::face_t *> &b);
+
+ void _generateEdgeEdgeIntersections(carve::mesh::MeshSet<3>::edge_t *ea,
+ carve::mesh::MeshSet<3>::edge_t *eb);
+ void generateEdgeEdgeIntersections(carve::mesh::MeshSet<3>::face_t *a,
+ const std::vector<carve::mesh::MeshSet<3>::face_t *> &b);
+
+ void _generateVertexFaceIntersections(carve::mesh::MeshSet<3>::face_t *fa,
+ carve::mesh::MeshSet<3>::edge_t *eb);
+ void generateVertexFaceIntersections(carve::mesh::MeshSet<3>::face_t *a,
+ const std::vector<carve::mesh::MeshSet<3>::face_t *> &b);
+
+ void _generateEdgeFaceIntersections(carve::mesh::MeshSet<3>::face_t *fa,
+ carve::mesh::MeshSet<3>::edge_t *eb);
+ void generateEdgeFaceIntersections(carve::mesh::MeshSet<3>::face_t *a,
+ const std::vector<carve::mesh::MeshSet<3>::face_t *> &b);
+
+ void generateIntersectionCandidates(carve::mesh::MeshSet<3> *a,
+ const face_rtree_t *a_node,
+ carve::mesh::MeshSet<3> *b,
+ const face_rtree_t *b_node,
+ face_pairs_t &face_pairs,
+ bool descend_a = true);
+ /**
+ * \brief Compute all points of intersection between poly \a a and poly \a b
+ *
+ * @param a Polyhedron a.
+ * @param b Polyhedron b.
+ */
+ void generateIntersections(meshset_t *a,
+ const face_rtree_t *a_node,
+ meshset_t *b,
+ const face_rtree_t *b_node,
+ detail::Data &data);
+
+ /**
+ * \brief Generate tables of intersecting pairs of faces.
+ *
+ * @param[out] data Internal data-structure holding intersection info.
+ */
+ void intersectingFacePairs(detail::Data &data);
+
+ /**
+ * \brief Divide edges in \a edges that are intersected by polyhedron \a poly
+ *
+ * @param edges The edges to divide.
+ * @param[in] poly The polyhedron to divide against.
+ * @param[in,out] data Intersection information.
+ */
+ void divideEdges(
+ const std::vector<meshset_t::edge_t> &edges,
+ meshset_t *poly,
+ detail::Data &data);
+
+ void divideIntersectedEdges(detail::Data &data);
+
+ /**
+ * \brief From the intersection points of pairs of intersecting faces, compute intersection edges.
+ *
+ * @param[out] eclass Classification information about created edges.
+ * @param[in,out] data Intersection information.
+ */
+ void makeFaceEdges(
+ EdgeClassification &eclass,
+ detail::Data &data);
+
+ friend void classifyEasyFaces(
+ FaceLoopList &face_loops,
+ VertexClassification &vclass,
+ meshset_t *other_poly,
+ int other_poly_num,
+ CSG &csg,
+ CSG::Collector &collector);
+
+ size_t generateFaceLoops(
+ meshset_t *poly,
+ const detail::Data &data,
+ FaceLoopList &face_loops_out);
+
+
+
+ // intersect_group.cpp
+
+ /**
+ * \brief Build a loop edge mapping from a list of face loops.
+ *
+ * @param[in] loops A list of face loops.
+ * @param[in] edge_count A hint as to the number of edges in \a loops.
+ * @param[out] edge_map The calculated map of edges to loops.
+ */
+ void makeEdgeMap(
+ const FaceLoopList &loops,
+ size_t edge_count,
+ detail::LoopEdges &edge_map);
+
+ /**
+ * \brief Divide a list of face loops into groups that are connected by at least one edge not present in \a no_cross.
+ *
+ * @param[in] src The source mesh from which these loops derive.
+ * @param[in,out] face_loops The list of loops (will be emptied as a side effect)
+ * @param[in] loop_edges A loop edge map used for traversing connected loops.
+ * @param[in] no_cross A set of edges not to cross.
+ * @param[out] out_loops A list of grouped face loops.
+ */
+ void groupFaceLoops(
+ carve::mesh::MeshSet<3> *src,
+ FaceLoopList &face_loops,
+ const detail::LoopEdges &loop_edges,
+ const V2Set &no_cross,
+ FLGroupList &out_loops);
+
+ /**
+ * \brief Find the set of edges shared between two edge maps.
+ *
+ * @param[in] edge_map_a The first edge map.
+ * @param[in] edge_map_b The second edge map.
+ * @param[out] shared_edges The resulting set of common edges.
+ */
+ void findSharedEdges(
+ const detail::LoopEdges &edge_map_a,
+ const detail::LoopEdges &edge_map_b,
+ V2Set &shared_edges);
+
+
+ // intersect_classify_edge.cpp
+
+ /**
+ *
+ *
+ * @param shared_edges
+ * @param vclass
+ * @param poly_a
+ * @param a_loops_grouped
+ * @param a_edge_map
+ * @param poly_b
+ * @param b_loops_grouped
+ * @param b_edge_map
+ * @param collector
+ */
+ void classifyFaceGroupsEdge(
+ const V2Set &shared_edges,
+ VertexClassification &vclass,
+ meshset_t *poly_a,
+ const face_rtree_t *poly_a_rtree,
+ FLGroupList &a_loops_grouped,
+ const detail::LoopEdges &a_edge_map,
+ meshset_t *poly_b,
+ const face_rtree_t *poly_b_rtree,
+ FLGroupList &b_loops_grouped,
+ const detail::LoopEdges &b_edge_map,
+ CSG::Collector &collector);
+
+ // intersect_classify_group.cpp
+
+ /**
+ *
+ *
+ * @param shared_edges
+ * @param vclass
+ * @param poly_a
+ * @param a_loops_grouped
+ * @param a_edge_map
+ * @param poly_b
+ * @param b_loops_grouped
+ * @param b_edge_map
+ * @param collector
+ */
+ void classifyFaceGroups(
+ const V2Set &shared_edges,
+ VertexClassification &vclass,
+ meshset_t *poly_a,
+ const face_rtree_t *poly_a_rtree,
+ FLGroupList &a_loops_grouped,
+ const detail::LoopEdges &a_edge_map,
+ meshset_t *poly_b,
+ const face_rtree_t *poly_b_rtree,
+ FLGroupList &b_loops_grouped,
+ const detail::LoopEdges &b_edge_map,
+ CSG::Collector &collector);
+
+ // intersect_half_classify_group.cpp
+
+ /**
+ *
+ *
+ * @param shared_edges
+ * @param vclass
+ * @param poly_a
+ * @param a_loops_grouped
+ * @param a_edge_map
+ * @param poly_b
+ * @param b_loops_grouped
+ * @param b_edge_map
+ * @param FaceClass
+ * @param b_out
+ */
+ void halfClassifyFaceGroups(
+ const V2Set &shared_edges,
+ VertexClassification &vclass,
+ meshset_t *poly_a,
+ const face_rtree_t *poly_a_rtree,
+ FLGroupList &a_loops_grouped,
+ const detail::LoopEdges &a_edge_map,
+ meshset_t *poly_b,
+ const face_rtree_t *poly_b_rtree,
+ FLGroupList &b_loops_grouped,
+ const detail::LoopEdges &b_edge_map,
+ std::list<std::pair<FaceClass, meshset_t *> > &b_out);
+
+ // intersect.cpp
+
+ /**
+ * \brief The main calculation method for CSG.
+ *
+ * @param[in] a Polyhedron a
+ * @param[in] b Polyhedron b
+ * @param[out] vclass
+ * @param[out] eclass
+ * @param[out] a_face_loops
+ * @param[out] b_face_loops
+ * @param[out] a_edge_count
+ * @param[out] b_edge_count
+ */
+ void calc(
+ meshset_t *a,
+ const face_rtree_t *a_rtree,
+ meshset_t *b,
+ const face_rtree_t *b_rtree,
+ VertexClassification &vclass,
+ EdgeClassification &eclass,
+ FaceLoopList &a_face_loops,
+ FaceLoopList &b_face_loops,
+ size_t &a_edge_count,
+ size_t &b_edge_count);
+
+ public:
+ /**
+ * \enum OP
+ * \brief Enumeration of the supported CSG operations.
+ */
+ enum OP {
+ UNION, /**< in a or b. */
+ INTERSECTION, /**< in a and b. */
+ A_MINUS_B, /**< in a, but not b. */
+ B_MINUS_A, /**< in b, but not a. */
+ SYMMETRIC_DIFFERENCE, /**< in a or b, but not both. */
+ ALL /**< all split faces from a and b */
+ };
+
+ /**
+ * \enum CLASSIFY_TYPE
+ * \brief The type of classification algorithm to use.
+ */
+ enum CLASSIFY_TYPE {
+ CLASSIFY_NORMAL, /**< Normal (group) classifier. */
+ CLASSIFY_EDGE /**< Edge classifier. */
+ };
+
+ CSG::Hooks hooks; /**< The manager for calculation hooks. */
+
+ CSG();
+ ~CSG();
+
+ /**
+ * \brief Compute a CSG operation between two polyhedra, \a a and \a b.
+ *
+ * @param a Polyhedron a
+ * @param b Polyhedron b
+ * @param collector The collector (determines the CSG operation performed)
+ * @param shared_edges A pointer to a set that will be populated with shared edges (if not NULL).
+ * @param classify_type The type of classifier to use.
+ *
+ * @return
+ */
+ meshset_t *compute(
+ meshset_t *a,
+ meshset_t *b,
+ CSG::Collector &collector,
+ V2Set *shared_edges = NULL,
+ CLASSIFY_TYPE classify_type = CLASSIFY_NORMAL);
+
+ /**
+ * \brief Compute a CSG operation between two closed polyhedra, \a a and \a b.
+ *
+ * @param a Polyhedron a
+ * @param b Polyhedron b
+ * @param op The CSG operation (A collector is created automatically).
+ * @param shared_edges A pointer to a set that will be populated with shared edges (if not NULL).
+ * @param classify_type The type of classifier to use.
+ *
+ * @return
+ */
+ meshset_t *compute(
+ meshset_t *a,
+ meshset_t *b,
+ OP op,
+ V2Set *shared_edges = NULL,
+ CLASSIFY_TYPE classify_type = CLASSIFY_NORMAL);
+
+ void slice(
+ meshset_t *a,
+ meshset_t *b,
+ std::list<meshset_t *> &a_sliced,
+ std::list<meshset_t *> &b_sliced,
+ V2Set *shared_edges = NULL);
+
+ bool sliceAndClassify(
+ meshset_t *closed,
+ meshset_t *open,
+ std::list<std::pair<FaceClass, meshset_t *> > &result,
+ V2Set *shared_edges = NULL);
+ };
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+#include <carve/csg.hpp>
+#include <carve/tag.hpp>
+#include <carve/poly.hpp>
+#include <carve/triangulator.hpp>
+#include <deque>
+
+namespace carve {
+ namespace csg {
+
+ namespace detail {
+ template<bool with_improvement>
+ class CarveTriangulator : public csg::CSG::Hook {
+
+ public:
+ CarveTriangulator() {
+ }
+
+ virtual ~CarveTriangulator() {
+ }
+
+ virtual void processOutputFace(std::vector<carve::mesh::MeshSet<3>::face_t *> &faces,
+ const carve::mesh::MeshSet<3>::face_t *orig,
+ bool flipped) {
+ std::vector<carve::mesh::MeshSet<3>::face_t *> out_faces;
+
+ size_t n_tris = 0;
+ for (size_t f = 0; f < faces.size(); ++f) {
+ CARVE_ASSERT(faces[f]->nVertices() >= 3);
+ n_tris += faces[f]->nVertices() - 2;
+ }
+
+ out_faces.reserve(n_tris);
+
+ for (size_t f = 0; f < faces.size(); ++f) {
+ carve::mesh::MeshSet<3>::face_t *face = faces[f];
+
+ if (face->nVertices() == 3) {
+ out_faces.push_back(face);
+ continue;
+ }
+
+ std::vector<triangulate::tri_idx> result;
+
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> vloop;
+ face->getVertices(vloop);
+
+ triangulate::triangulate(
+ carve::mesh::MeshSet<3>::face_t::projection_mapping(face->project),
+ vloop,
+ result);
+
+ if (with_improvement) {
+ triangulate::improve(
+ carve::mesh::MeshSet<3>::face_t::projection_mapping(face->project),
+ vloop,
+ carve::mesh::vertex_distance(),
+ result);
+ }
+
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> fv;
+ fv.resize(3);
+ for (size_t i = 0; i < result.size(); ++i) {
+ fv[0] = vloop[result[i].a];
+ fv[1] = vloop[result[i].b];
+ fv[2] = vloop[result[i].c];
+ out_faces.push_back(face->create(fv.begin(), fv.end(), false));
+ }
+ delete face;
+ }
+ std::swap(faces, out_faces);
+ }
+ };
+ }
+
+ typedef detail::CarveTriangulator<false> CarveTriangulator;
+ typedef detail::CarveTriangulator<true> CarveTriangulatorWithImprovement;
+
+ class CarveTriangulationImprover : public csg::CSG::Hook {
+ public:
+ CarveTriangulationImprover() {
+ }
+
+ virtual ~CarveTriangulationImprover() {
+ }
+
+ virtual void processOutputFace(std::vector<carve::mesh::MeshSet<3>::face_t *> &faces,
+ const carve::mesh::MeshSet<3>::face_t *orig,
+ bool flipped) {
+ if (faces.size() == 1) return;
+
+ // doing improvement as a separate hook is much messier than
+ // just incorporating it into the triangulation hook.
+
+ typedef std::map<carve::mesh::MeshSet<3>::vertex_t *, size_t> vert_map_t;
+ std::vector<carve::mesh::MeshSet<3>::face_t *> out_faces;
+ vert_map_t vert_map;
+
+ out_faces.reserve(faces.size());
+
+
+ carve::mesh::MeshSet<3>::face_t::projection_mapping projector(faces[0]->project);
+
+ std::vector<triangulate::tri_idx> result;
+
+ for (size_t f = 0; f < faces.size(); ++f) {
+ carve::mesh::MeshSet<3>::face_t *face = faces[f];
+ if (face->nVertices() != 3) {
+ out_faces.push_back(face);
+ } else {
+ triangulate::tri_idx tri;
+ for (carve::mesh::MeshSet<3>::face_t::edge_iter_t i = face->begin(); i != face->end(); ++i) {
+ size_t v = 0;
+ vert_map_t::iterator j = vert_map.find(i->vert);
+ if (j == vert_map.end()) {
+ v = vert_map.size();
+ vert_map[i->vert] = v;
+ } else {
+ v = (*j).second;
+ }
+ tri.v[i.idx()] = v;
+ }
+ result.push_back(tri);
+ delete face;
+ }
+ }
+
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> verts;
+ verts.resize(vert_map.size());
+ for (vert_map_t::iterator i = vert_map.begin(); i != vert_map.end(); ++i) {
+ verts[(*i).second] = (*i).first;
+ }
+
+ triangulate::improve(projector, verts, carve::mesh::vertex_distance(), result);
+
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> fv;
+ fv.resize(3);
+ for (size_t i = 0; i < result.size(); ++i) {
+ fv[0] = verts[result[i].a];
+ fv[1] = verts[result[i].b];
+ fv[2] = verts[result[i].c];
+ out_faces.push_back(orig->create(fv.begin(), fv.end(), false));
+ }
+
+ std::swap(faces, out_faces);
+ }
+ };
+
+ class CarveTriangulationQuadMerger : public csg::CSG::Hook {
+ // this code is incomplete.
+ typedef std::map<V2, F2> edge_map_t;
+
+ public:
+ CarveTriangulationQuadMerger() {
+ }
+
+ virtual ~CarveTriangulationQuadMerger() {
+ }
+
+ double scoreQuad(edge_map_t::iterator i, edge_map_t &edge_map) {
+ if (!(*i).second.first || !(*i).second.second) return -1;
+ return 0;
+ }
+
+ carve::mesh::MeshSet<3>::face_t *mergeQuad(edge_map_t::iterator i, edge_map_t &edge_map) {
+ return NULL;
+ }
+
+ void recordEdge(carve::mesh::MeshSet<3>::vertex_t *v1,
+ carve::mesh::MeshSet<3>::vertex_t *v2,
+ carve::mesh::MeshSet<3>::face_t *f,
+ edge_map_t &edge_map) {
+ if (v1 < v2) {
+ edge_map[V2(v1, v2)].first = f;
+ } else {
+ edge_map[V2(v2, v1)].second = f;
+ }
+ }
+
+ virtual void processOutputFace(std::vector<carve::mesh::MeshSet<3>::face_t *> &faces,
+ const carve::mesh::MeshSet<3>::face_t *orig,
+ bool flipped) {
+ if (faces.size() == 1) return;
+
+ std::vector<carve::mesh::MeshSet<3>::face_t *> out_faces;
+ edge_map_t edge_map;
+
+ out_faces.reserve(faces.size());
+
+ poly::p2_adapt_project<3> projector(faces[0]->project);
+
+ for (size_t f = 0; f < faces.size(); ++f) {
+ carve::mesh::MeshSet<3>::face_t *face = faces[f];
+ if (face->nVertices() != 3) {
+ out_faces.push_back(face);
+ } else {
+ carve::mesh::MeshSet<3>::face_t::vertex_t *v1, *v2, *v3;
+ v1 = face->edge->vert;
+ v2 = face->edge->next->vert;
+ v3 = face->edge->next->next->vert;
+ recordEdge(v1, v2, face, edge_map);
+ recordEdge(v2, v3, face, edge_map);
+ recordEdge(v3, v1, face, edge_map);
+ }
+ }
+
+ for (edge_map_t::iterator i = edge_map.begin(); i != edge_map.end();) {
+ if ((*i).second.first && (*i).second.second) {
+ ++i;
+ } else {
+ edge_map.erase(i++);
+ }
+ }
+
+ while (edge_map.size()) {
+ edge_map_t::iterator i = edge_map.begin();
+ edge_map_t::iterator best = i;
+ double best_score = scoreQuad(i, edge_map);
+ for (++i; i != edge_map.end(); ++i) {
+ double score = scoreQuad(i, edge_map);
+ if (score > best_score) best = i;
+ }
+ if (best_score < 0) break;
+ out_faces.push_back(mergeQuad(best, edge_map));
+ }
+
+ if (edge_map.size()) {
+ tagable::tag_begin();
+ for (edge_map_t::iterator i = edge_map.begin(); i != edge_map.end(); ++i) {
+ carve::mesh::MeshSet<3>::face_t *a = const_cast<carve::mesh::MeshSet<3>::face_t *>((*i).second.first);
+ carve::mesh::MeshSet<3>::face_t *b = const_cast<carve::mesh::MeshSet<3>::face_t *>((*i).second.first);
+ if (a && a->tag_once()) out_faces.push_back(a);
+ if (b && b->tag_once()) out_faces.push_back(b);
+ }
+ }
+
+ std::swap(faces, out_faces);
+ }
+ };
+
+ class CarveHoleResolver : public csg::CSG::Hook {
+
+ public:
+ CarveHoleResolver() {
+ }
+
+ virtual ~CarveHoleResolver() {
+ }
+
+ bool findRepeatedEdges(const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &vertices,
+ std::list<std::pair<size_t, size_t> > &edge_pos) {
+ std::map<V2, size_t> edges;
+ for (size_t i = 0; i < vertices.size() - 1; ++i) {
+ edges[std::make_pair(vertices[i], vertices[i+1])] = i;
+ }
+ edges[std::make_pair(vertices[vertices.size()-1], vertices[0])] = vertices.size() - 1;
+
+ for (std::map<V2, size_t>::iterator i = edges.begin(); i != edges.end(); ++i) {
+ V2 rev = V2((*i).first.second, (*i).first.first);
+ std::map<V2, size_t>::iterator j = edges.find(rev);
+ if (j != edges.end()) {
+ edge_pos.push_back(std::make_pair((*i).second, (*j).second));
+ }
+ }
+ return edge_pos.size() > 0;
+ }
+
+ void flood(size_t t1,
+ size_t t2,
+ size_t old_grp,
+ size_t new_grp_1,
+ size_t new_grp_2,
+ std::vector<size_t> &grp,
+ const std::vector<triangulate::tri_idx> &tris,
+ const std::map<std::pair<size_t, size_t>, size_t> &tri_edge) {
+ grp[t1] = new_grp_1;
+ grp[t2] = new_grp_2;
+
+ std::deque<size_t> to_visit;
+ to_visit.push_back(t1);
+ to_visit.push_back(t2);
+ std::vector<std::pair<size_t, size_t> > rev;
+ rev.resize(3);
+ while (to_visit.size()) {
+ size_t curr = to_visit.front();
+ to_visit.pop_front();
+ triangulate::tri_idx ct = tris[curr];
+ rev[0] = std::make_pair(ct.b, ct.a);
+ rev[1] = std::make_pair(ct.c, ct.b);
+ rev[2] = std::make_pair(ct.a, ct.c);
+
+ for (size_t i = 0; i < 3; ++i) {
+ std::map<std::pair<size_t, size_t>, size_t>::const_iterator adj = tri_edge.find(rev[i]);
+ if (adj == tri_edge.end()) continue;
+ size_t next = (*adj).second;
+ if (grp[next] != old_grp) continue;
+ grp[next] = grp[curr];
+ to_visit.push_back(next);
+ }
+ }
+ }
+
+ void findPerimeter(const std::vector<triangulate::tri_idx> &tris,
+ const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &verts,
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> &out) {
+ std::map<std::pair<size_t, size_t>, size_t> edges;
+ for (size_t i = 0; i < tris.size(); ++i) {
+ edges[std::make_pair(tris[i].a, tris[i].b)] = i;
+ edges[std::make_pair(tris[i].b, tris[i].c)] = i;
+ edges[std::make_pair(tris[i].c, tris[i].a)] = i;
+ }
+ std::map<size_t, size_t> unpaired;
+ for (std::map<std::pair<size_t, size_t>, size_t>::iterator i = edges.begin(); i != edges.end(); ++i) {
+ if (edges.find(std::make_pair((*i).first.second, (*i).first.first)) == edges.end()) {
+ CARVE_ASSERT(unpaired.find((*i).first.first) == unpaired.end());
+ unpaired[(*i).first.first] = (*i).first.second;
+ }
+ }
+ out.clear();
+ out.reserve(unpaired.size());
+ size_t start = (*unpaired.begin()).first;
+ size_t vert = start;
+ do {
+ out.push_back(verts[vert]);
+ CARVE_ASSERT(unpaired.find(vert) != unpaired.end());
+ vert = unpaired[vert];
+ } while (vert != start);
+ }
+
+ virtual void processOutputFace(std::vector<carve::mesh::MeshSet<3>::face_t *> &faces,
+ const carve::mesh::MeshSet<3>::face_t *orig,
+ bool flipped) {
+ std::vector<carve::mesh::MeshSet<3>::face_t *> out_faces;
+
+ for (size_t f = 0; f < faces.size(); ++f) {
+ carve::mesh::MeshSet<3>::face_t *face = faces[f];
+
+ if (face->nVertices() == 3) {
+ out_faces.push_back(face);
+ continue;
+ }
+
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> vloop;
+ face->getVertices(vloop);
+
+ std::list<std::pair<size_t, size_t> > rep_edges;
+ if (!findRepeatedEdges(vloop, rep_edges)) {
+ out_faces.push_back(face);
+ continue;
+ }
+
+ std::vector<triangulate::tri_idx> result;
+ triangulate::triangulate(
+ carve::mesh::MeshSet<3>::face_t::projection_mapping(face->project),
+ vloop,
+ result);
+
+ std::map<std::pair<size_t, size_t>, size_t> tri_edge;
+ for (size_t i = 0; i < result.size(); ++i) {
+ tri_edge[std::make_pair(result[i].a, result[i].b)] = i;
+ tri_edge[std::make_pair(result[i].b, result[i].c)] = i;
+ tri_edge[std::make_pair(result[i].c, result[i].a)] = i;
+ }
+
+ std::vector<size_t> grp;
+ grp.resize(result.size(), 0);
+
+ size_t grp_max = 0;
+
+ while (rep_edges.size()) {
+ std::pair<size_t, size_t> e1, e2;
+
+ e1.first = rep_edges.front().first;
+ e1.second = (e1.first + 1) % vloop.size();
+
+ e2.first = rep_edges.front().second;
+ e2.second = (e2.first + 1) % vloop.size();
+
+ rep_edges.pop_front();
+
+ CARVE_ASSERT(tri_edge.find(e1) != tri_edge.end());
+ size_t t1 = tri_edge[e1];
+ CARVE_ASSERT(tri_edge.find(e2) != tri_edge.end());
+ size_t t2 = tri_edge[e2];
+
+ if (grp[t1] != grp[t2]) {
+ continue;
+ }
+
+ size_t t1g = ++grp_max;
+ size_t t2g = ++grp_max;
+
+ flood(t1, t2, grp[t1], t1g, t2g, grp, result, tri_edge);
+ }
+
+ std::set<size_t> groups;
+ std::copy(grp.begin(), grp.end(), std::inserter(groups, groups.begin()));
+
+ // now construct perimeters for each group.
+ std::vector<triangulate::tri_idx> grp_tris;
+ grp_tris.reserve(result.size());
+ for (std::set<size_t>::iterator i = groups.begin(); i != groups.end(); ++i) {
+ size_t grp_id = *i;
+ grp_tris.clear();
+ for (size_t j = 0; j < grp.size(); ++j) {
+ if (grp[j] == grp_id) {
+ grp_tris.push_back(result[j]);
+ }
+ }
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> grp_perim;
+ findPerimeter(grp_tris, vloop, grp_perim);
+ out_faces.push_back(face->create(grp_perim.begin(), grp_perim.end(), false));
+ }
+ }
+ std::swap(faces, out_faces);
+ }
+ };
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/vector.hpp>
+#include <carve/geom3d.hpp>
+#include <carve/csg.hpp>
+
+#include <iomanip>
+
+template<typename MAP>
+void map_histogram(std::ostream &out, const MAP &map) {
+ std::vector<int> hist;
+ for (typename MAP::const_iterator i = map.begin(); i != map.end(); ++i) {
+ size_t n = (*i).second.size();
+ if (hist.size() <= n) {
+ hist.resize(n + 1);
+ }
+ hist[n]++;
+ }
+ int total = map.size();
+ std::string bar(50, '*');
+ for (size_t i = 0; i < hist.size(); i++) {
+ if (hist[i] > 0) {
+ out << std::setw(5) << i << " : " << std::setw(5) << hist[i] << " " << bar.substr(50 - hist[i] * 50 / total) << std::endl;
+ }
+ }
+}
+
+namespace carve {
+ namespace csg {
+ class IntersectDebugHooks {
+ public:
+ virtual void drawIntersections(const VertexIntersections & /* vint */) {
+ }
+
+ virtual void drawPoint(const carve::mesh::MeshSet<3>::vertex_t * /* v */,
+ float /* r */,
+ float /* g */,
+ float /* b */,
+ float /* a */,
+ float /* rad */) {
+ }
+ virtual void drawEdge(const carve::mesh::MeshSet<3>::vertex_t * /* v1 */,
+ const carve::mesh::MeshSet<3>::vertex_t * /* v2 */,
+ float /* rA */, float /* gA */, float /* bA */, float /* aA */,
+ float /* rB */, float /* gB */, float /* bB */, float /* aB */,
+ float /* thickness */ = 1.0) {
+ }
+
+ virtual void drawFaceLoopWireframe(const std::vector<carve::mesh::MeshSet<3>::vertex_t *> & /* face_loop */,
+ const carve::mesh::MeshSet<3>::vertex_t & /* normal */,
+ float /* r */, float /* g */, float /* b */, float /* a */,
+ bool /* inset */ = true) {
+ }
+
+ virtual void drawFaceLoop(const std::vector<carve::mesh::MeshSet<3>::vertex_t *> & /* face_loop */,
+ const carve::mesh::MeshSet<3>::vertex_t & /* normal */,
+ float /* r */, float /* g */, float /* b */, float /* a */,
+ bool /* offset */ = true,
+ bool /* lit */ = true) {
+ }
+
+ virtual void drawFaceLoop2(const std::vector<carve::mesh::MeshSet<3>::vertex_t *> & /* face_loop */,
+ const carve::mesh::MeshSet<3>::vertex_t & /* normal */,
+ float /* rF */, float /* gF */, float /* bF */, float /* aF */,
+ float /* rB */, float /* gB */, float /* bB */, float /* aB */,
+ bool /* offset */ = true,
+ bool /* lit */ = true) {
+ }
+
+ virtual ~IntersectDebugHooks() {
+ }
+ };
+
+ IntersectDebugHooks *intersect_installDebugHooks(IntersectDebugHooks *hooks);
+ bool intersect_debugEnabled();
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <vector>
+
+
+
+namespace carve {
+namespace djset {
+
+
+
+ class djset {
+
+ protected:
+ struct elem {
+ size_t parent, rank;
+ elem(size_t p, size_t r) : parent(p), rank(r) {}
+ elem() {}
+ };
+
+ std::vector<elem> set;
+ size_t n_sets;
+
+ public:
+ djset() : set(), n_sets(0) {
+ }
+
+ djset(size_t N) {
+ n_sets = N;
+ set.reserve(N);
+ for (size_t i = 0; i < N; ++i) {
+ set.push_back(elem(i,0));
+ }
+ }
+
+ void init(size_t N) {
+ if (N == set.size()) {
+ for (size_t i = 0; i < N; ++i) {
+ set[i] = elem(i,0);
+ }
+ n_sets = N;
+ } else {
+ djset temp(N);
+ std::swap(set, temp.set);
+ std::swap(n_sets, temp.n_sets);
+ }
+ }
+
+ size_t count() const {
+ return n_sets;
+ }
+
+ size_t find_set_head(size_t a) {
+ if (a == set[a].parent) return a;
+
+ size_t a_head = a;
+ while (set[a_head].parent != a_head) a_head = set[a_head].parent;
+ set[a].parent = a_head;
+ return a_head;
+ }
+
+ bool same_set(size_t a, size_t b) {
+ return find_set_head(a) == find_set_head(b);
+ }
+
+ void merge_sets(size_t a, size_t b) {
+ a = find_set_head(a);
+ b = find_set_head(b);
+ if (a != b) {
+ n_sets--;
+ if (set[a].rank < set[b].rank) {
+ set[a].parent = b;
+ } else if (set[b].rank < set[a].rank) {
+ set[b].parent = a;
+ } else {
+ set[a].rank++;
+ set[b].parent = a;
+ }
+ }
+ }
+
+ void get_index_to_set(std::vector<size_t> &index_set, std::vector<size_t> &set_size) {
+ index_set.clear();
+ index_set.resize(set.size(), n_sets);
+ set_size.clear();
+ set_size.resize(n_sets, 0);
+
+ size_t c = 0;
+ for (size_t i = 0; i < set.size(); ++i) {
+ size_t s = find_set_head(i);
+ if (index_set[s] == n_sets) index_set[s] = c++;
+ index_set[i] = index_set[s];
+ set_size[index_set[s]]++;
+ }
+ }
+
+ template<typename in_iter_t, typename out_collection_t>
+ void collate(in_iter_t in, out_collection_t &out) {
+ std::vector<size_t> set_id(set.size(), n_sets);
+ out.clear();
+ out.resize(n_sets);
+ size_t c = 0;
+ for (size_t i = 0; i < set.size(); ++i) {
+ size_t s = find_set_head(i);
+ if (set_id[s] == n_sets) set_id[s] = c++;
+ s = set_id[s];
+ std::insert_iterator<typename out_collection_t::value_type> j(out[s], out[s].end());
+ *j = *in++;
+ }
+ }
+ };
+
+
+
+}
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/vector.hpp>
+#include <carve/tag.hpp>
+
+#include <vector>
+#include <list>
+
+namespace carve {
+ namespace poly {
+
+
+
+ struct Object;
+
+
+
+ template<unsigned ndim>
+ class Edge : public tagable {
+ public:
+ typedef Vertex<ndim> vertex_t;
+ typedef typename Vertex<ndim>::vector_t vector_t;
+ typedef Object obj_t;
+
+ const vertex_t *v1, *v2;
+ const obj_t *owner;
+
+ Edge(const vertex_t *_v1, const vertex_t *_v2, const obj_t *_owner) :
+ tagable(), v1(_v1), v2(_v2), owner(_owner) {
+ }
+
+ ~Edge() {
+ }
+ };
+
+
+
+ struct hash_edge_ptr {
+ template<unsigned ndim>
+ size_t operator()(const Edge<ndim> * const &e) const {
+ return (size_t)e;
+ }
+ };
+
+
+
+ }
+}
+
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+namespace carve {
+ namespace poly {
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <vector>
+#include <numeric>
+#include <algorithm>
+
+
+namespace carve {
+ namespace exact {
+
+ class exact_t : public std::vector<double> {
+ typedef std::vector<double> super;
+
+ public:
+ exact_t() : super() {
+ }
+
+ exact_t(double v, size_t sz = 1) : super(sz, v) {
+ }
+
+ template<typename iter_t>
+ exact_t(iter_t a, iter_t b) : super(a, b) {
+ }
+
+ exact_t(double a, double b) : super() {
+ reserve(2);
+ push_back(a);
+ push_back(b);
+ }
+
+ exact_t(double a, double b, double c) : super() {
+ reserve(3);
+ push_back(a);
+ push_back(b);
+ push_back(c);
+ }
+
+ exact_t(double a, double b, double c, double d) : super() {
+ reserve(4);
+ push_back(a);
+ push_back(b);
+ push_back(c);
+ push_back(d);
+ }
+
+ exact_t(double a, double b, double c, double d, double e) : super() {
+ reserve(5);
+ push_back(a);
+ push_back(b);
+ push_back(c);
+ push_back(d);
+ push_back(e);
+ }
+
+ exact_t(double a, double b, double c, double d, double e, double f) : super() {
+ reserve(6);
+ push_back(a);
+ push_back(b);
+ push_back(c);
+ push_back(d);
+ push_back(e);
+ push_back(f);
+ }
+
+ exact_t(double a, double b, double c, double d, double e, double f, double g) : super() {
+ reserve(7);
+ push_back(a);
+ push_back(b);
+ push_back(c);
+ push_back(d);
+ push_back(e);
+ push_back(f);
+ push_back(g);
+ }
+
+ exact_t(double a, double b, double c, double d, double e, double f, double g, double h) : super() {
+ reserve(8);
+ push_back(a);
+ push_back(b);
+ push_back(c);
+ push_back(d);
+ push_back(e);
+ push_back(f);
+ push_back(g);
+ push_back(h);
+ }
+
+ void compress();
+
+ exact_t compressed() const {
+ exact_t result(*this);
+ result.compress();
+ return result;
+ }
+
+ operator double() const {
+ return std::accumulate(begin(), end(), 0.0);
+ }
+
+ void removeZeroes() {
+ erase(std::remove(begin(), end(), 0.0), end());
+ }
+ };
+
+ inline std::ostream &operator<<(std::ostream &out, const exact_t &p) {
+ out << '{';
+ out << p[0];
+ for (size_t i = 1; i < p.size(); ++i) out << ';' << p[i];
+ out << '}';
+ return out;
+ }
+
+
+
+ namespace detail {
+ const struct constants_t {
+ double splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */
+ double epsilon; /* = 2^(-p). Used to estimate roundoff errors. */
+ /* A set of coefficients used to calculate maximum roundoff errors. */
+ double resulterrbound;
+ double ccwerrboundA, ccwerrboundB, ccwerrboundC;
+ double o3derrboundA, o3derrboundB, o3derrboundC;
+ double iccerrboundA, iccerrboundB, iccerrboundC;
+ double isperrboundA, isperrboundB, isperrboundC;
+
+ constants_t() {
+ double half;
+ double check, lastcheck;
+ int every_other;
+
+ every_other = 1;
+ half = 0.5;
+ epsilon = 1.0;
+ splitter = 1.0;
+ check = 1.0;
+ /* Repeatedly divide `epsilon' by two until it is too small to add to */
+ /* one without causing roundoff. (Also check if the sum is equal to */
+ /* the previous sum, for machines that round up instead of using exact */
+ /* rounding. Not that this library will work on such machines anyway. */
+ do {
+ lastcheck = check;
+ epsilon *= half;
+ if (every_other) {
+ splitter *= 2.0;
+ }
+ every_other = !every_other;
+ check = 1.0 + epsilon;
+ } while ((check != 1.0) && (check != lastcheck));
+ splitter += 1.0;
+
+ /* Error bounds for orientation and incircle tests. */
+ resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
+ ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
+ ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
+ ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
+ o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
+ o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
+ o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
+ iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
+ iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
+ iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
+ isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
+ isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
+ isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
+ }
+ } constants;
+
+ template<unsigned U, unsigned V>
+ struct op {
+ enum {
+ Vlo = V / 2,
+ Vhi = V - Vlo
+ };
+
+ static inline void add(const double *a, const double *b, double *r) {
+ double t[U + Vlo];
+ op<U, Vlo>::add(a, b, t);
+ for (size_t i = 0; i < Vlo; ++i) r[i] = t[i];
+ op<U, Vhi>::add(t + Vlo, b + Vlo, r + Vlo);
+ }
+
+ static inline void sub(const double *a, const double *b, double *r) {
+ double t[U + Vlo];
+ op<U, Vlo>::sub(a, b, t);
+ for (size_t i = 0; i < Vlo; ++i) r[i] = t[i];
+ op<U, Vhi>::sub(t + Vlo, b + Vlo, r + Vlo);
+ }
+ };
+
+ template<unsigned U>
+ struct op<U, 1> {
+ enum {
+ Ulo = U / 2,
+ Uhi = U - Ulo
+ };
+ static void add(const double *a, const double *b, double *r) {
+ double t[Ulo + 1];
+ op<Ulo, 1>::add(a, b, t);
+ for (size_t i = 0; i < Ulo; ++i) r[i] = t[i];
+ op<Uhi, 1>::add(a + Ulo, t + Ulo, r + Ulo);
+ }
+
+ static void sub(const double *a, const double *b, double *r) {
+ double t[Ulo + 1];
+ op<Ulo, 1>::sub(a, b, t);
+ for (size_t i = 0; i < Ulo; ++i) r[i] = t[i];
+ op<Uhi, 1>::add(a + Ulo, t + Ulo, r + Ulo);
+ }
+ };
+
+ template<>
+ struct op<1, 1> {
+ static void add_fast(const double *a, const double *b, double *r) {
+ assert(fabs(a[0]) >= fabs(b[0]));
+ volatile double sum = a[0] + b[0];
+ volatile double bvirt = sum - a[0];
+ r[0] = b[0] - bvirt;
+ r[1] = sum;
+ }
+
+ static void sub_fast(const double *a, const double *b, double *r) {
+ assert(fabs(a[0]) >= fabs(b[0]));
+ volatile double diff = a[0] - b[0];
+ volatile double bvirt = a[0] - diff;
+ r[0] = bvirt - b[0];
+ r[1] = diff;
+ }
+
+ static void add(const double *a, const double *b, double *r) {
+ volatile double sum = a[0] + b[0];
+ volatile double bvirt = sum - a[0];
+ double avirt = sum - bvirt;
+ double bround = b[0] - bvirt;
+ double around = a[0] - avirt;
+ r[0] = around + bround;
+ r[1] = sum;
+ }
+
+ static void sub(const double *a, const double *b, double *r) {
+ volatile double diff = a[0] - b[0];
+ volatile double bvirt = a[0] - diff;
+ double avirt = diff + bvirt;
+ double bround = bvirt - b[0];
+ double around = a[0] - avirt;
+ r[0] = around + bround;
+ r[1] = diff;
+ }
+ };
+
+
+ template<unsigned U, unsigned V>
+ static exact_t add(const double *a, const double *b) {
+ exact_t result;
+ result.resize(U + V);
+ op<U,V>::add(a, b, &result[0]);
+ return result;
+ }
+
+
+ template<unsigned U, unsigned V>
+ static exact_t sub(const double *a, const double *b) {
+ exact_t result;
+ result.resize(U + V);
+ op<U,V>::sub(a, b, &result[0]);
+ return result;
+ }
+
+
+ template<unsigned U, unsigned V>
+ static exact_t add(const exact_t &a, const exact_t &b) {
+ assert(a.size() == U);
+ assert(b.size() == V);
+ exact_t result;
+ result.resize(U + V);
+ std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
+ op<U,V>::add(&a[0], &b[0], &result[0]);
+ return result;
+ }
+
+
+ template<unsigned U, unsigned V>
+ static exact_t add(const exact_t &a, const double *b) {
+ assert(a.size() == U);
+ exact_t result;
+ result.resize(U + V);
+ std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
+ op<U,V>::add(&a[0], b, &result[0]);
+ return result;
+ }
+
+
+ template<unsigned U, unsigned V>
+ static exact_t sub(const exact_t &a, const exact_t &b) {
+ assert(a.size() == U);
+ assert(b.size() == V);
+ exact_t result;
+ result.resize(U + V);
+ std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
+ op<U,V>::sub(&a[0], &b[0], &result[0]);
+ return result;
+ }
+
+
+ template<unsigned U, unsigned V>
+ static exact_t sub(const exact_t &a, const double *b) {
+ assert(a.size() == U);
+ exact_t result;
+ result.resize(U + V);
+ std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
+ op<U,V>::sub(&a[0], &b[0], &result[0]);
+ return result;
+ }
+
+
+ static inline void split(const double a, double *r) {
+ volatile double c = constants.splitter * a;
+ volatile double abig = c - a;
+ r[1] = c - abig;
+ r[0] = a - r[1];
+ }
+
+ static inline void prod_1_1(const double *a, const double *b, double *r) {
+ r[1] = a[0] * b[0];
+ double a_sp[2]; split(a[0], a_sp);
+ double b_sp[2]; split(b[0], b_sp);
+ double err1 = r[1] - a_sp[1] * b_sp[1];
+ double err2 = err1 - a_sp[0] * b_sp[1];
+ double err3 = err2 - a_sp[1] * b_sp[0];
+ r[0] = a_sp[0] * b_sp[0] - err3;
+ }
+
+ static inline void prod_1_1s(const double *a, const double *b, const double *b_sp, double *r) {
+ r[1] = a[0] * b[0];
+ double a_sp[2]; split(a[0], a_sp);
+ double err1 = r[1] - a_sp[1] * b_sp[1];
+ double err2 = err1 - a_sp[0] * b_sp[1];
+ double err3 = err2 - a_sp[1] * b_sp[0];
+ r[0] = a_sp[0] * b_sp[0] - err3;
+ }
+
+ static inline void prod_1s_1s(const double *a, const double *a_sp, const double *b, const double *b_sp, double *r) {
+ r[1] = a[0] * b[0];
+ double err1 = r[1] - a_sp[1] * b_sp[1];
+ double err2 = err1 - a_sp[0] * b_sp[1];
+ double err3 = err2 - a_sp[1] * b_sp[0];
+ r[0] = a_sp[0] * b_sp[0] - err3;
+ }
+
+ static inline void prod_2_1(const double *a, const double *b, double *r) {
+ double b_sp[2]; split(b[0], b_sp);
+ double t1[2]; prod_1_1s(a+0, b, b_sp, t1);
+ r[0] = t1[0];
+ double t2[2]; prod_1_1s(a+1, b, b_sp, t2);
+ double t3[2]; op<1,1>::add(t1+1, t2, t3);
+ r[1] = t3[0];
+ double t4[2]; op<1,1>::add_fast(t2+1, t3+1, r + 2);
+ }
+
+ static inline void prod_1_2(const double *a, const double *b, double *r) {
+ prod_2_1(b, a, r);
+ }
+
+ static inline void prod_4_1(const double *a, const double *b, double *r) {
+ double b_sp[2]; split(b[0], b_sp);
+ double t1[2]; prod_1_1s(a+0, b, b_sp, t1);
+ r[0] = t1[0];
+ double t2[2]; prod_1_1s(a+1, b, b_sp, t2);
+ double t3[2]; op<1,1>::add(t1+1, t2, t3);
+ r[1] = t3[0];
+ double t4[2]; op<1,1>::add_fast(t2+1, t3+1, t4);
+ r[2] = t4[0];
+ double t5[2]; prod_1_1s(a+2, b, b_sp, t5);
+ double t6[2]; op<1,1>::add(t4+1, t5, t6);
+ r[3] = t6[0];
+ double t7[2]; op<1,1>::add_fast(t5+1, t6+1, t7);
+ r[4] = t7[0];
+ double t8[2]; prod_1_1s(a+3, b, b_sp, t8);
+ double t9[2]; op<1,1>::add(t7+1, t8, t9);
+ r[5] = t9[0];
+ op<1,1>::add_fast(t8+1, t9+1, r + 6);
+ }
+
+ static inline void prod_1_4(const double *a, const double *b, double *r) {
+ prod_4_1(b, a, r);
+ }
+
+ static inline void prod_2_2(const double *a, const double *b, double *r) {
+ double a1_sp[2]; split(a[1], a1_sp);
+ double a0_sp[2]; split(a[0], a0_sp);
+ double b1_sp[2]; split(b[1], b1_sp);
+ double b0_sp[2]; split(b[0], b0_sp);
+
+ double t1[2]; prod_1s_1s(a+0, a0_sp, b+0, b0_sp, t1);
+ r[0] = t1[0];
+ double t2[2]; prod_1s_1s(a+1, a1_sp, b+0, b0_sp, t2);
+
+ double t3[2]; op<1,1>::add(t1+1, t2, t3);
+ double t4[2]; op<1,1>::add_fast(t2+1, t3+1, t4);
+
+ double t5[2]; prod_1s_1s(a+0, a0_sp, b+1, b1_sp, t5);
+
+ double t6[2]; op<1,1>::add(t3, t5, t6);
+ r[1] = t6[0];
+ double t7[2]; op<1,1>::add(t4, t6+1, t7);
+ double t8[2]; op<1,1>::add(t4+1, t7+1, t8);
+
+ double t9[2]; prod_1s_1s(a+1, a1_sp, b+1, b1_sp, t9);
+
+ double t10[2]; op<1,1>::add(t5+1, t9, t10);
+ double t11[2]; op<1,1>::add(t7, t10, t11);
+ r[2] = t11[0];
+ double t12[2]; op<1,1>::add(t8, t11+1, t12);
+ double t13[2]; op<1,1>::add(t8+1, t12+1, t13);
+ double t14[2]; op<1,1>::add(t9+1, t10+1, t14);
+ double t15[2]; op<1,1>::add(t12, t14, t15);
+ r[3] = t15[0];
+ double t16[2]; op<1,1>::add(t13, t15+1, t16);
+ double t17[2]; op<1,1>::add(t13+1, t16+1, t17);
+ double t18[2]; op<1,1>::add(t16, t14+1, t18);
+ r[4] = t18[0];
+ double t19[2]; op<1,1>::add(t17, t18+1, t19);
+ r[5] = t19[0];
+ double t20[2]; op<1,1>::add(t17+1, t19+1, t20);
+ r[6] = t20[0];
+ r[7] = t20[1];
+ }
+
+
+
+ static inline void square(const double a, double *r) {
+ r[1] = a * a;
+ double a_sp[2]; split(a, a_sp);
+ double err1 = r[1] - (a_sp[1] * a_sp[1]);
+ double err3 = err1 - ((a_sp[1] + a_sp[1]) * a_sp[0]);
+ r[0] = a_sp[0] * a_sp[0] - err3;
+ }
+
+ static inline void square_2(const double *a, double *r) {
+ double t1[2]; square(a[0], t1);
+ r[0] = t1[0];
+ double t2 = a[0] + a[0];
+ double t3[2]; prod_1_1(a+1, &t2, t3);
+ double t4[3]; op<2,1>::add(t3, t1 + 1, t4);
+ r[1] = t4[0];
+ double t5[2]; square(a[1], t5);
+ double t6[4]; op<2,2>::add(t5, t4 + 1, r + 2);
+ }
+ }
+
+
+
+ void exact_t::compress() {
+ double sum[2];
+
+ int j = size() - 1;
+ double Q = (*this)[j];
+ for (int i = (int)size()-2; i >= 0; --i) {
+ detail::op<1,1>::add_fast(&Q, &(*this)[i], sum);
+ if (sum[0] != 0) {
+ (*this)[j--] = sum[1];
+ Q = sum[0];
+ } else {
+ Q = sum[1];
+ }
+ }
+ int j2 = 0;
+ for (int i = j + 1; i < (int)size(); ++i) {
+ detail::op<1,1>::add_fast(&(*this)[i], &Q, sum);
+ if (sum[0] != 0) {
+ (*this)[j2++] = sum[0];
+ }
+ Q = sum[1];
+ }
+ (*this)[j2++] = Q;
+
+ erase(begin() + j2, end());
+ }
+
+ template<typename iter_t>
+ void negate(iter_t begin, iter_t end) {
+ while (begin != end) { *begin = -*begin; ++begin; }
+ }
+
+ void negate(exact_t &e) {
+ negate(&e[0], &e[e.size()]);
+ }
+
+ template<typename iter_t>
+ void scale_zeroelim(iter_t ebegin,
+ iter_t eend,
+ double b,
+ exact_t &h) {
+ double Q;
+
+ h.clear();
+ double b_sp[2]; detail::split(b, b_sp);
+
+ double prod[2], sum[2];
+
+ detail::prod_1_1s((double *)ebegin++, &b, b_sp, prod);
+ Q = prod[1];
+ if (prod[0] != 0.0) {
+ h.push_back(prod[0]);
+ }
+ while (ebegin != eend) {
+ double enow = *ebegin++;
+ detail::prod_1_1s(&enow, &b, b_sp, prod);
+ detail::op<1,1>::add(&Q, prod, sum);
+ if (sum[0] != 0) {
+ h.push_back(sum[0]);
+ }
+ detail::op<1,1>::add_fast(prod+1, sum+1, sum);
+ Q = sum[1];
+ if (sum[0] != 0) {
+ h.push_back(sum[0]);
+ }
+ }
+ if ((Q != 0.0) || (h.size() == 0)) {
+ h.push_back(Q);
+ }
+ }
+
+ void scale_zeroelim(const exact_t &e,
+ double b,
+ exact_t &h) {
+ scale_zeroelim(&e[0], &e[e.size()], b, h);
+ }
+
+ template<typename iter_t>
+ void sum_zeroelim(iter_t ebegin,
+ iter_t eend,
+ iter_t fbegin,
+ iter_t fend,
+ exact_t &h) {
+ double Q;
+ double enow, fnow;
+
+ double sum[2];
+
+ enow = *ebegin;
+ fnow = *fbegin;
+
+ h.clear();
+
+ if ((fnow > enow) == (fnow > -enow)) {
+ Q = enow;
+ enow = *++ebegin;
+ } else {
+ Q = fnow;
+ fnow = *++fbegin;
+ }
+
+ if (ebegin != eend && fbegin != fend) {
+ if ((fnow > enow) == (fnow > -enow)) {
+ detail::op<1,1>::add_fast(&enow, &Q, sum);
+ enow = *++ebegin;
+ } else {
+ detail::op<1,1>::add_fast(&fnow, &Q, sum);
+ fnow = *++fbegin;
+ }
+ Q = sum[1];
+ if (sum[0] != 0.0) {
+ h.push_back(sum[0]);
+ }
+ while (ebegin != eend && fbegin != fend) {
+ if ((fnow > enow) == (fnow > -enow)) {
+ detail::op<1,1>::add(&Q, &enow, sum);
+ enow = *++ebegin;
+ } else {
+ detail::op<1,1>::add(&Q, &fnow, sum);
+ fnow = *++fbegin;
+ }
+ Q = sum[1];
+ if (sum[0] != 0.0) {
+ h.push_back(sum[0]);
+ }
+ }
+ }
+
+ while (ebegin != eend) {
+ detail::op<1,1>::add(&Q, &enow, sum);
+ enow = *++ebegin;
+ Q = sum[1];
+ if (sum[0] != 0.0) {
+ h.push_back(sum[0]);
+ }
+ }
+ while (fbegin != fend) {
+ detail::op<1,1>::add(&Q, &fnow, sum);
+ fnow = *++fbegin;
+ Q = sum[1];
+ if (sum[0] != 0.0) {
+ h.push_back(sum[0]);
+ }
+ }
+
+ if (Q != 0.0 || !h.size()) {
+ h.push_back(Q);
+ }
+ }
+
+ void sum_zeroelim(const exact_t &e,
+ const exact_t &f,
+ exact_t &h) {
+ sum_zeroelim(&e[0], &e[e.size()], &f[0], &f[f.size()], h);
+ }
+
+ void sum_zeroelim(const double *ebegin,
+ const double *eend,
+ const exact_t &f,
+ exact_t &h) {
+ sum_zeroelim(ebegin, eend, &f[0], &f[f.size()], h);
+ }
+
+ void sum_zeroelim(const exact_t &e,
+ const double *fbegin,
+ const double *fend,
+ exact_t &h) {
+ sum_zeroelim(&e[0], &e[e.size()], fbegin, fend, h);
+ }
+
+
+ // XXX: not implemented yet
+ //exact_t operator+(const exact_t &a, const exact_t &b) {
+ //}
+
+
+
+ void diffprod(const double a, const double b, const double c, const double d, double *r) {
+ // return ab - cd;
+ double ab[2], cd[2];
+ detail::prod_1_1(&a, &b, ab);
+ detail::prod_1_1(&c, &d, cd);
+ detail::op<2,2>::sub(ab, cd, r);
+ }
+
+ double orient3dexact(const double *pa,
+ const double *pb,
+ const double *pc,
+ const double *pd) {
+ using namespace detail;
+
+ double ab[4]; diffprod(pa[0], pb[1], pb[0], pa[1], ab);
+ double bc[4]; diffprod(pb[0], pc[1], pc[0], pb[1], bc);
+ double cd[4]; diffprod(pc[0], pd[1], pd[0], pc[1], cd);
+ double da[4]; diffprod(pd[0], pa[1], pa[0], pd[1], da);
+ double ac[4]; diffprod(pa[0], pc[1], pc[0], pa[1], ac);
+ double bd[4]; diffprod(pb[0], pd[1], pd[0], pb[1], bd);
+
+ exact_t temp;
+ exact_t cda, dab, abc, bcd;
+ exact_t adet, bdet, cdet, ddet, abdet, cddet, det;
+
+ sum_zeroelim(cd, cd + 4, da, da + 4, temp);
+ sum_zeroelim(temp, ac, ac + 4, cda);
+
+ sum_zeroelim(da, da + 4, ab, ab + 4, temp);
+ sum_zeroelim(temp, bd, bd + 4, dab);
+
+ negate(bd, bd + 4);
+ negate(ac, bd + 4);
+
+ sum_zeroelim(ab, ab + 4, bc, bc + 4, temp);
+ sum_zeroelim(temp, ac, ac + 4, abc);
+
+ sum_zeroelim(bc, bc + 4, cd, cd + 4, temp);
+ sum_zeroelim(temp, bd, bd + 4, bcd);
+
+ scale_zeroelim(bcd, +pa[2], adet);
+ scale_zeroelim(cda, -pb[2], bdet);
+ scale_zeroelim(dab, +pc[2], cdet);
+ scale_zeroelim(abc, -pd[2], ddet);
+
+ sum_zeroelim(adet, bdet, abdet);
+ sum_zeroelim(cdet, ddet, cddet);
+
+ sum_zeroelim(abdet, cddet, det);
+
+ return det[det.size() - 1];
+ }
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/geom2d.hpp>
+#include <carve/vector.hpp>
+#include <carve/matrix.hpp>
+#include <carve/geom3d.hpp>
+#include <carve/aabb.hpp>
+#include <carve/tag.hpp>
+
+#include <vector>
+#include <list>
+#include <map>
+
+namespace carve {
+ namespace poly {
+
+
+
+ struct Object;
+
+ template<unsigned ndim>
+ class Edge;
+
+
+
+ template<unsigned ndim>
+ struct p2_adapt_project {
+ typedef carve::geom2d::P2 (*proj_t)(const carve::geom::vector<ndim> &);
+ proj_t proj;
+ p2_adapt_project(proj_t _proj) : proj(_proj) { }
+ carve::geom2d::P2 operator()(const carve::geom::vector<ndim> &v) const { return proj(v); }
+ carve::geom2d::P2 operator()(const carve::geom::vector<ndim> *v) const { return proj(*v); }
+ carve::geom2d::P2 operator()(const Vertex<ndim> &v) const { return proj(v.v); }
+ carve::geom2d::P2 operator()(const Vertex<ndim> *v) const { return proj(v->v); }
+ };
+
+
+ template<unsigned ndim>
+ class Face : public tagable {
+ public:
+ typedef Vertex<ndim> vertex_t;
+ typedef typename Vertex<ndim>::vector_t vector_t;
+ typedef Edge<ndim> edge_t;
+ typedef Object obj_t;
+ typedef carve::geom::aabb<ndim> aabb_t;
+ typedef carve::geom::plane<ndim> plane_t;
+
+ typedef carve::geom2d::P2 (*project_t)(const vector_t &);
+ typedef vector_t (*unproject_t)(const carve::geom2d::P2 &, const plane_t &);
+
+ protected:
+ std::vector<const vertex_t *> vertices; // pointer into polyhedron.vertices
+ std::vector<const edge_t *> edges; // pointer into polyhedron.edges
+
+ project_t getProjector(bool positive_facing, int axis);
+ unproject_t getUnprojector(bool positive_facing, int axis);
+
+ public:
+ typedef typename std::vector<const vertex_t *>::iterator vertex_iter_t;
+ typedef typename std::vector<const vertex_t *>::const_iterator const_vertex_iter_t;
+
+ typedef typename std::vector<const edge_t *>::iterator edge_iter_t;
+ typedef typename std::vector<const edge_t *>::const_iterator const_edge_iter_t;
+
+ obj_t *owner;
+
+ aabb_t aabb;
+ plane_t plane_eqn;
+ int manifold_id;
+ int group_id;
+
+ project_t project;
+ unproject_t unproject;
+
+ Face(const std::vector<const vertex_t *> &_vertices, bool delay_recalc = false);
+ Face(const vertex_t *v1, const vertex_t *v2, const vertex_t *v3, bool delay_recalc = false);
+ Face(const vertex_t *v1, const vertex_t *v2, const vertex_t *v3, const vertex_t *v4, bool delay_recalc = false);
+
+ template <typename iter_t>
+ Face(const Face *base, iter_t vbegin, iter_t vend, bool flipped) {
+ init(base, vbegin, vend, flipped);
+ }
+
+ Face(const Face *base, const std::vector<const vertex_t *> &_vertices, bool flipped) {
+ init(base, _vertices, flipped);
+ }
+
+ Face() {}
+ ~Face() {}
+
+ bool recalc();
+
+ template<typename iter_t>
+ Face *init(const Face *base, iter_t vbegin, iter_t vend, bool flipped);
+ Face *init(const Face *base, const std::vector<const vertex_t *> &_vertices, bool flipped);
+
+ template<typename iter_t>
+ Face *create(iter_t vbegin, iter_t vend, bool flipped) const;
+ Face *create(const std::vector<const vertex_t *> &_vertices, bool flipped) const;
+
+ Face *clone(bool flipped = false) const;
+ void invert();
+
+ void getVertexLoop(std::vector<const vertex_t *> &loop) const;
+
+ const vertex_t *&vertex(size_t idx);
+ const vertex_t *vertex(size_t idx) const;
+ size_t nVertices() const;
+
+ vertex_iter_t vbegin() { return vertices.begin(); }
+ vertex_iter_t vend() { return vertices.end(); }
+ const_vertex_iter_t vbegin() const { return vertices.begin(); }
+ const_vertex_iter_t vend() const { return vertices.end(); }
+
+ std::vector<carve::geom::vector<2> > projectedVertices() const;
+
+ const edge_t *&edge(size_t idx);
+ const edge_t *edge(size_t idx) const;
+ size_t nEdges() const;
+
+ edge_iter_t ebegin() { return edges.begin(); }
+ edge_iter_t eend() { return edges.end(); }
+ const_edge_iter_t ebegin() const { return edges.begin(); }
+ const_edge_iter_t eend() const { return edges.end(); }
+
+ bool containsPoint(const vector_t &p) const;
+ bool containsPointInProjection(const vector_t &p) const;
+ bool simpleLineSegmentIntersection(const carve::geom::linesegment<ndim> &line,
+ vector_t &intersection) const;
+ IntersectionClass lineSegmentIntersection(const carve::geom::linesegment<ndim> &line,
+ vector_t &intersection) const;
+ vector_t centroid() const;
+
+ p2_adapt_project<ndim> projector() const {
+ return p2_adapt_project<ndim>(project);
+ }
+
+ void swap(Face<ndim> &other);
+ };
+
+
+
+ struct hash_face_ptr {
+ template<unsigned ndim>
+ size_t operator()(const Face<ndim> * const &f) const {
+ return (size_t)f;
+ }
+ };
+
+
+
+ namespace face {
+
+
+
+ template<unsigned ndim>
+ static inline carve::geom2d::P2 project(const Face<ndim> *f, const typename Face<ndim>::vector_t &v) {
+ return f->project(v);
+ }
+
+
+
+ template<unsigned ndim>
+ static inline carve::geom2d::P2 project(const Face<ndim> &f, const typename Face<ndim>::vector_t &v) {
+ return f.project(v);
+ }
+
+
+
+ template<unsigned ndim>
+ static inline typename Face<ndim>::vector_t unproject(const Face<ndim> *f, const carve::geom2d::P2 &p) {
+ return f->unproject(p, f->plane_eqn);
+ }
+
+
+
+ template<unsigned ndim>
+ static inline typename Face<ndim>::vector_t unproject(const Face<ndim> &f, const carve::geom2d::P2 &p) {
+ return f.unproject(p, f.plane_eqn);
+ }
+
+
+
+ }
+
+
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+namespace std {
+ template<unsigned ndim>
+ inline void swap(carve::poly::Face<ndim> &a, carve::poly::Face<ndim> &b) {
+ a.swap(b);
+ }
+}
+
+namespace carve {
+ namespace poly {
+ template<unsigned ndim>
+ void Face<ndim>::swap(Face<ndim> &other) {
+ std::swap(vertices, other.vertices);
+ std::swap(edges, other.edges);
+ std::swap(owner, other.owner);
+ std::swap(aabb, other.aabb);
+ std::swap(plane_eqn, other.plane_eqn);
+ std::swap(manifold_id, other.manifold_id);
+ std::swap(group_id, other.group_id);
+ std::swap(project, other.project);
+ std::swap(unproject, other.unproject);
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t>
+ Face<ndim> *Face<ndim>::init(const Face<ndim> *base, iter_t vbegin, iter_t vend, bool flipped) {
+ vertices.reserve(std::distance(vbegin, vend));
+
+ if (flipped) {
+ std::reverse_copy(vbegin, vend, std::back_inserter(vertices));
+ plane_eqn = -base->plane_eqn;
+ } else {
+ std::copy(vbegin, vend, std::back_inserter(vertices));
+ plane_eqn = base->plane_eqn;
+ }
+
+ edges.clear();
+ edges.resize(nVertices(), NULL);
+
+ aabb.fit(vertices.begin(), vertices.end(), vec_adapt_vertex_ptr());
+ untag();
+
+ int da = carve::geom::largestAxis(plane_eqn.N);
+
+ project = getProjector(plane_eqn.N.v[da] > 0, da);
+ unproject = getUnprojector(plane_eqn.N.v[da] > 0, da);
+
+ return this;
+ }
+
+ template<unsigned ndim>
+ template<typename iter_t>
+ Face<ndim> *Face<ndim>::create(iter_t vbegin, iter_t vend, bool flipped) const {
+ return (new Face)->init(this, vbegin, vend, flipped);
+ }
+
+ template<unsigned ndim>
+ Face<ndim> *Face<ndim>::create(const std::vector<const vertex_t *> &_vertices, bool flipped) const {
+ return (new Face)->init(this, _vertices.begin(), _vertices.end(), flipped);
+ }
+
+ template<unsigned ndim>
+ Face<ndim> *Face<ndim>::clone(bool flipped) const {
+ return (new Face)->init(this, vertices, flipped);
+ }
+
+ template<unsigned ndim>
+ void Face<ndim>::getVertexLoop(std::vector<const vertex_t *> &loop) const {
+ loop.resize(nVertices(), NULL);
+ std::copy(vbegin(), vend(), loop.begin());
+ }
+
+ template<unsigned ndim>
+ const typename Face<ndim>::edge_t *&Face<ndim>::edge(size_t idx) {
+ return edges[idx];
+ }
+
+ template<unsigned ndim>
+ const typename Face<ndim>::edge_t *Face<ndim>::edge(size_t idx) const {
+ return edges[idx];
+ }
+
+ template<unsigned ndim>
+ size_t Face<ndim>::nEdges() const {
+ return edges.size();
+ }
+
+ template<unsigned ndim>
+ const typename Face<ndim>::vertex_t *&Face<ndim>::vertex(size_t idx) {
+ return vertices[idx];
+ }
+
+ template<unsigned ndim>
+ const typename Face<ndim>::vertex_t *Face<ndim>::vertex(size_t idx) const {
+ return vertices[idx];
+ }
+
+ template<unsigned ndim>
+ size_t Face<ndim>::nVertices() const {
+ return vertices.size();
+ }
+
+ template<unsigned ndim>
+ typename Face<ndim>::vector_t Face<ndim>::centroid() const {
+ vector_t c;
+ carve::geom::centroid(vertices.begin(), vertices.end(), vec_adapt_vertex_ptr(), c);
+ return c;
+ }
+
+ template<unsigned ndim>
+ std::vector<carve::geom::vector<2> > Face<ndim>::projectedVertices() const {
+ p2_adapt_project<ndim> proj = projector();
+ std::vector<carve::geom::vector<2> > result;
+ result.reserve(nVertices());
+ for (size_t i = 0; i < nVertices(); ++i) {
+ result.push_back(proj(vertex(i)->v));
+ }
+ return result;
+ }
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/classification.hpp>
+#include <carve/collection_types.hpp>
+
+namespace carve {
+ namespace csg {
+
+ struct FaceLoopGroup;
+
+ struct FaceLoop {
+ FaceLoop *next, *prev;
+ const carve::mesh::MeshSet<3>::face_t *orig_face;
+ std::vector<carve::mesh::MeshSet<3>::vertex_t *> vertices;
+ FaceLoopGroup *group;
+
+ FaceLoop(const carve::mesh::MeshSet<3>::face_t *f, const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &v) : next(NULL), prev(NULL), orig_face(f), vertices(v), group(NULL) {}
+ };
+
+
+ struct FaceLoopList {
+ FaceLoop *head, *tail;
+ unsigned count;
+
+ FaceLoopList() : head(NULL), tail(NULL), count(0) { }
+
+ void append(FaceLoop *f) {
+ f->prev = tail;
+ f->next = NULL;
+ if (tail) tail->next = f;
+ tail = f;
+ if (!head) head = f;
+ count++;
+ }
+
+ void prepend(FaceLoop *f) {
+ f->next = head;
+ f->prev = NULL;
+ if (head) head->prev = f;
+ head = f;
+ if (!tail) tail = f;
+ count++;
+ }
+
+ unsigned size() const {
+ return count;
+ }
+
+ FaceLoop *remove(FaceLoop *f) {
+ FaceLoop *r = f->next;
+ if (f->prev) { f->prev->next = f->next; } else { head = f->next; }
+ if (f->next) { f->next->prev = f->prev; } else { tail = f->prev; }
+ f->next = f->prev = NULL;
+ count--;
+ return r;
+ }
+
+ ~FaceLoopList() {
+ FaceLoop *a = head, *b;
+ while (a) {
+ b = a;
+ a = a->next;
+ delete b;
+ }
+ }
+ };
+
+ struct FaceLoopGroup {
+ const carve::mesh::MeshSet<3> *src;
+ FaceLoopList face_loops;
+ V2Set perimeter;
+ std::list<ClassificationInfo> classification;
+
+ FaceLoopGroup(const carve::mesh::MeshSet<3> *_src) : src(_src) {
+ }
+
+ FaceClass classificationAgainst(const carve::mesh::MeshSet<3>::mesh_t *mesh) const;
+ };
+
+
+
+ typedef std::list<FaceLoopGroup> FLGroupList;
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <vector>
+
+namespace carve {
+ namespace geom {
+
+ template<unsigned ndim> struct aabb;
+
+ // ========================================================================
+ struct _uninitialized { };
+
+ template<unsigned ndim>
+ struct base {
+ double v[ndim];
+ };
+
+ template<> struct base<2> {union { double v[2]; struct { double x, y; }; }; };
+ template<> struct base<3> {union { double v[3]; struct { double x, y, z; }; }; };
+ template<> struct base<4> {union { double v[4]; struct { double x, y, z, w; }; }; };
+
+ template<unsigned ndim>
+ struct vector : public base<ndim> {
+ enum { __ndim = ndim };
+
+ static vector ZERO();
+ double length2() const;
+ double length() const;
+ vector<ndim> &normalize();
+ vector<ndim> normalized() const;
+ bool exactlyZero() const;
+ bool isZero(double epsilon = EPSILON) const;
+ void setZero();
+ void fill(double val);
+ vector<ndim> &scaleBy(double d);
+ vector<ndim> &invscaleBy(double d);
+ vector<ndim> scaled(double d) const;
+ vector<ndim> invscaled(double d) const;
+ vector<ndim> &negate();
+ vector<ndim> negated() const;
+ double &operator[](unsigned i);
+ const double &operator[](unsigned i) const;
+ template<typename assign_t>
+ vector<ndim> &operator=(const assign_t &t);
+ std::string asStr() const;
+
+ aabb<ndim> getAABB() const;
+
+ vector() { setZero(); }
+ vector(noinit_t) { }
+ };
+
+ template<unsigned ndim>
+ vector<ndim> vector<ndim>::ZERO() { vector<ndim> r; r.setZero(); return r; }
+
+ static inline vector<2> VECTOR(double x, double y) { vector<2> r; r.x = x; r.y = y; return r; }
+ static inline vector<3> VECTOR(double x, double y, double z) { vector<3> r; r.x = x; r.y = y; r.z = z; return r; }
+ static inline vector<4> VECTOR(double x, double y, double z, double w) { vector<4> r; r.x = x; r.y = y; r.z = z; r.w = w; return r; }
+
+ template<unsigned ndim> vector<ndim> operator+(const vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> vector<ndim> operator+(const vector<ndim> &a, double b);
+ template<unsigned ndim, typename val_t> vector<ndim> operator+(const vector<ndim> &a, const val_t &b);
+ template<unsigned ndim, typename val_t> vector<ndim> operator+(const val_t &a, const vector<ndim> &b);
+
+ template<unsigned ndim> vector<ndim> &operator+=(vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> vector<ndim> &operator+=(vector<ndim> &a, double b);
+ template<unsigned ndim, typename val_t> vector<ndim> &operator+=(vector<ndim> &a, const val_t &b);
+
+ template<unsigned ndim> vector<ndim> operator-(const vector<ndim> &a);
+
+ template<unsigned ndim> vector<ndim> operator-(const vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> vector<ndim> operator-(const vector<ndim> &a, double b);
+ template<unsigned ndim, typename val_t> vector<ndim> operator-(const vector<ndim> &a, const val_t &b);
+ template<unsigned ndim, typename val_t> vector<ndim> operator-(const val_t &a, const vector<ndim> &b);
+
+ template<unsigned ndim> vector<ndim> &operator-=(vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> vector<ndim> &operator-=(vector<ndim> &a, double b);
+ template<unsigned ndim, typename val_t> vector<ndim> &operator-=(vector<ndim> &a, const val_t &b);
+
+ template<unsigned ndim> vector<ndim> operator*(const vector<ndim> &a, double s);
+ template<unsigned ndim> vector<ndim> operator*(double s, const vector<ndim> &a);
+ template<unsigned ndim> vector<ndim> &operator*=(vector<ndim> &a, double s);
+
+ template<unsigned ndim> vector<ndim> operator/(const vector<ndim> &a, double s);
+ template<unsigned ndim> vector<ndim> &operator/=(vector<ndim> &a, double s);
+
+ template<unsigned ndim> bool operator==(const vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator!=(const vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator<(const vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator<=(const vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator>(const vector<ndim> &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator>=(const vector<ndim> &a, const vector<ndim> &b);
+
+ template<unsigned ndim> vector<ndim> abs(const vector<ndim> &a);
+
+ template<unsigned ndim> double distance2(const vector<ndim> &a, const vector<ndim> &b);
+
+ template<unsigned ndim> double distance(const vector<ndim> &a, const vector<ndim> &b);
+
+ template<unsigned ndim> bool equal(const vector<ndim> &a, const vector<ndim> &b);
+
+ template<unsigned ndim> int smallestAxis(const vector<ndim> &a);
+
+ template<unsigned ndim> int largestAxis(const vector<ndim> &a);
+
+ template<unsigned ndim> vector<2> select(const vector<ndim> &a, int a1, int a2);
+
+ template<unsigned ndim> vector<3> select(const vector<ndim> &a, int a1, int a2, int a3);
+
+ template<unsigned ndim, typename assign_t, typename oper_t>
+ vector<ndim> &assign_op(vector<ndim> &a, const assign_t &t, oper_t op);
+
+ template<unsigned ndim, typename assign1_t, typename assign2_t, typename oper_t>
+ vector<ndim> &assign_op(vector<ndim> &a, const assign1_t &t1, const assign2_t &t2, oper_t op);
+
+ template<unsigned ndim, typename iter_t>
+ void bounds(iter_t begin, iter_t end, vector<ndim> &min, vector<ndim> &max);
+
+ template<unsigned ndim, typename iter_t, typename adapt_t>
+ void bounds(iter_t begin, iter_t end, adapt_t adapt, vector<ndim> &min, vector<ndim> &max);
+
+ template<unsigned ndim, typename iter_t, typename adapt_t>
+ void centroid(iter_t begin, iter_t end, adapt_t adapt, vector<ndim> &c);
+
+ template<unsigned ndim, typename val_t> double dot(const vector<ndim> &a, const val_t &b);
+
+ static inline vector<3> cross(const vector<3> &a, const vector<3> &b);
+
+ static inline double cross(const vector<2> &a, const vector<2> &b);
+
+ static inline double dotcross(const vector<3> &a, const vector<3> &b, const vector<3> &c);
+
+
+
+ // ========================================================================
+ struct axis_pos {
+ int axis;
+ double pos;
+
+ axis_pos(int _axis, double _pos) : axis(_axis), pos(_pos) { }
+ };
+
+ template<unsigned ndim>
+ double distance(const axis_pos &a, const vector<ndim> &b);
+
+ template<unsigned ndim>
+ double distance2(const axis_pos &a, const vector<ndim> &b);
+
+ template<unsigned ndim> bool operator<(const axis_pos &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator<(const vector<ndim> &a, const axis_pos &b);
+
+ template<unsigned ndim> bool operator<=(const axis_pos &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator<=(const vector<ndim> &a, const axis_pos &b);
+
+ template<unsigned ndim> bool operator>(const axis_pos &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator>(const vector<ndim> &a, const axis_pos &b);
+
+ template<unsigned ndim> bool operator>=(const axis_pos &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator>=(const vector<ndim> &a, const axis_pos &b);
+
+ template<unsigned ndim> bool operator==(const axis_pos &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator==(const vector<ndim> &a, const axis_pos &b);
+
+ template<unsigned ndim> bool operator!=(const axis_pos &a, const vector<ndim> &b);
+ template<unsigned ndim> bool operator!=(const vector<ndim> &a, const axis_pos &b);
+
+
+
+ // ========================================================================
+ template<unsigned ndim>
+ struct ray {
+ typedef vector<ndim> vector_t;
+
+ vector_t D, v;
+
+ bool OK() const;
+
+ ray() { }
+ ray(vector_t _D, vector_t _v) : D(_D), v(_v) { }
+ };
+
+ template<unsigned ndim>
+ ray<ndim> rayThrough(const vector<ndim> &a, const vector<ndim> &b);
+
+ static inline double distance2(const ray<3> &r, const vector<3> &v);
+
+ static inline double distance(const ray<3> &r, const vector<3> &v);
+
+ static inline double distance2(const ray<2> &r, const vector<2> &v);
+
+ static inline double distance(const ray<2> &r, const vector<2> &v);
+
+
+
+ // ========================================================================
+ template<unsigned ndim>
+ struct linesegment {
+ typedef vector<ndim> vector_t;
+
+ vector_t v1;
+ vector_t v2;
+ vector_t midpoint;
+ vector_t half_length;
+
+ void update();
+ bool OK() const;
+ void flip();
+
+ aabb<ndim> getAABB() const;
+
+ linesegment(const vector_t &_v1, const vector_t &_v2);
+ };
+
+ template<unsigned ndim>
+ double distance2(const linesegment<ndim> &l, const vector<ndim> &v);
+
+ template<unsigned ndim>
+ double distance(const linesegment<ndim> &l, const vector<ndim> &v);
+
+
+
+ // ========================================================================
+ template<unsigned ndim>
+ struct plane {
+ typedef vector<ndim> vector_t;
+
+ vector_t N;
+ double d;
+
+ void negate();
+
+ plane();
+ plane(const vector_t &_N, vector_t _p);
+ plane(const vector_t &_N, double _d);
+ };
+
+ template<unsigned ndim>
+ inline plane<ndim> operator-(const plane<ndim> &p);
+
+ template<unsigned ndim, typename val_t>
+ double distance(const plane<ndim> &plane, const val_t &point);
+
+ template<unsigned ndim, typename val_t>
+ double distance2(const plane<ndim> &plane, const val_t &point);
+
+ template<unsigned ndim>
+ static inline vector<ndim> closestPoint(const plane<ndim> &p, const vector<ndim> &v);
+
+
+
+ // ========================================================================
+ template<unsigned ndim>
+ struct sphere {
+ typedef vector<ndim> vector_t;
+
+ vector_t C;
+ double r;
+
+ aabb<ndim> getAABB() const;
+
+ sphere();
+ sphere(const vector_t &_C, double _r);
+ };
+
+ template<unsigned ndim, typename val_t>
+ double distance(const sphere<ndim> &sphere, const val_t &point);
+
+ template<unsigned ndim, typename val_t>
+ double distance2(const sphere<ndim> &sphere, const val_t &point);
+
+ template<unsigned ndim>
+ static inline vector<ndim> closestPoint(const sphere<ndim> &sphere, const vector<ndim> &point);
+
+
+ // ========================================================================
+ template<unsigned ndim>
+ struct tri {
+ typedef vector<ndim> vector_t;
+
+ vector_t v[3];
+
+ aabb<ndim> getAABB() const;
+
+ tri(vector_t _v[3]);
+ tri(const vector_t &a, const vector_t &b, const vector_t &c);
+
+ vector_t normal() const {
+ return cross(v[1] - v[0], v[2] - v[1]).normalized();
+ }
+ };
+
+
+
+ template<unsigned ndim> std::ostream &operator<<(std::ostream &o, const vector<ndim> &v);
+ template<unsigned ndim> std::ostream &operator<<(std::ostream &o, const carve::geom::plane<ndim> &p);
+ template<unsigned ndim> std::ostream &operator<<(std::ostream &o, const carve::geom::sphere<ndim> &sphere);
+ template<unsigned ndim> std::ostream &operator<<(std::ostream &o, const carve::geom::tri<ndim> &tri);
+
+
+
+ template<unsigned ndim> vector<ndim> closestPoint(const tri<ndim> &tri, const vector<ndim> &pt);
+ template<unsigned ndim> double distance(const tri<ndim> &tri, const vector<ndim> &pt);
+ template<unsigned ndim> double distance2(const tri<ndim> &tri, const vector<ndim> &pt);
+
+
+
+ // ========================================================================
+ struct distance_functor {
+ template<typename obj1_t, typename obj2_t>
+ double operator()(const obj1_t &o1, const obj2_t &o2) {
+ return distance(o1, o2);
+ }
+ };
+
+
+
+ // ========================================================================
+ template<int base, int power> struct __pow__ { enum { val = __pow__<base, (power >> 1)>::val * __pow__<base, power - (power >> 1)>::val }; };
+ template<int base> struct __pow__<base, 1> { enum { val = base }; };
+ template<int base> struct __pow__<base, 0> { enum { val = 1 }; };
+
+ template<unsigned base, unsigned ndigits>
+ struct quantize {
+ typedef __pow__<base, ndigits> fac;
+
+ double operator()(double in) {
+ return round(in * fac::val) / fac::val;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> operator()(const vector<ndim> &in) {
+ vector<ndim> r(NOINIT);
+ assign_op(r, in, *this);
+ return r;
+ }
+ };
+
+
+
+ }
+}
+
+
+#include <carve/geom_impl.hpp>
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/math.hpp>
+#include <carve/math_constants.hpp>
+
+#include <carve/geom.hpp>
+
+#include <vector>
+#include <algorithm>
+
+#include <math.h>
+
+#if defined(CARVE_DEBUG)
+# include <iostream>
+#endif
+
+#if defined CARVE_USE_EXACT_PREDICATES
+# include <carve/shewchuk_predicates.hpp>
+#endif
+
+namespace carve {
+ namespace geom2d {
+
+ typedef carve::geom::vector<2> P2;
+ typedef carve::geom::ray<2> Ray2;
+ typedef carve::geom::linesegment<2> LineSegment2;
+
+
+
+ struct p2_adapt_ident {
+ P2 &operator()(P2 &p) const { return p; }
+ const P2 &operator()(const P2 &p) const { return p; }
+ };
+
+
+
+ typedef std::vector<P2> P2Vector;
+
+ /**
+ * \brief Return the orientation of c with respect to the ray defined by a->b.
+ *
+ * (Can be implemented exactly)
+ *
+ * @param[in] a
+ * @param[in] b
+ * @param[in] c
+ *
+ * @return positive, if c to the left of a->b.
+ * zero, if c is colinear with a->b.
+ * negative, if c to the right of a->b.
+ */
+#if defined CARVE_USE_EXACT_PREDICATES
+ inline double orient2d(const P2 &a, const P2 &b, const P2 &c) {
+ return shewchuk::orient2d(a.v, b.v, c.v);
+ }
+#else
+ inline double orient2d(const P2 &a, const P2 &b, const P2 &c) {
+ double acx = a.x - c.x;
+ double bcx = b.x - c.x;
+ double acy = a.y - c.y;
+ double bcy = b.y - c.y;
+ return acx * bcy - acy * bcx;
+ }
+#endif
+
+ /**
+ * \brief Determine whether p is internal to the anticlockwise
+ * angle abc, where b is the apex of the angle.
+ *
+ * @param[in] a
+ * @param[in] b
+ * @param[in] c
+ * @param[in] p
+ *
+ * @return true, if p is contained in the anticlockwise angle from
+ * b->a to b->c. Reflex angles contain p if p lies
+ * on b->a or on b->c. Acute angles do not contain p
+ * if p lies on b->a or on b->c. This is so that
+ * internalToAngle(a,b,c,p) = !internalToAngle(c,b,a,p)
+ */
+ inline bool internalToAngle(const P2 &a,
+ const P2 &b,
+ const P2 &c,
+ const P2 &p) {
+ bool reflex = (a < c) ? orient2d(b, a, c) <= 0.0 : orient2d(b, c, a) > 0.0;
+ double d1 = orient2d(b, a, p);
+ double d2 = orient2d(b, c, p);
+ if (reflex) {
+ return d1 >= 0.0 || d2 <= 0.0;
+ } else {
+ return d1 > 0.0 && d2 < 0.0;
+ }
+ }
+
+ /**
+ * \brief Determine whether p is internal to the anticlockwise
+ * angle ac, with apex at (0,0).
+ *
+ * @param[in] a
+ * @param[in] c
+ * @param[in] p
+ *
+ * @return true, if p is contained in a0c.
+ */
+ inline bool internalToAngle(const P2 &a,
+ const P2 &c,
+ const P2 &p) {
+ return internalToAngle(a, P2::ZERO(), c, p);
+ }
+
+ template<typename P2vec>
+ bool isAnticlockwise(const P2vec &tri) {
+ return orient2d(tri[0], tri[1], tri[2]) > 0.0;
+ }
+
+ template<typename P2vec>
+ bool pointIntersectsTriangle(const P2 &p, const P2vec &tri) {
+ int orient = isAnticlockwise(tri) ? +1 : -1;
+ if (orient2d(tri[0], tri[1], p) * orient < 0) return false;
+ if (orient2d(tri[1], tri[2], p) * orient < 0) return false;
+ if (orient2d(tri[2], tri[0], p) * orient < 0) return false;
+ return true;
+ }
+
+ template<typename P2vec>
+ bool lineIntersectsTriangle(const P2 &p1, const P2 &p2, const P2vec &tri) {
+ double s[3];
+ // does tri lie on one side or the other of p1-p2?
+ s[0] = orient2d(p1, p2, tri[0]);
+ s[1] = orient2d(p1, p2, tri[1]);
+ s[2] = orient2d(p1, p2, tri[2]);
+ if (*std::max_element(s, s+3) < 0) return false;
+ if (*std::min_element(s, s+3) > 0) return false;
+
+ // does line lie entirely to the right of a triangle edge?
+ int orient = isAnticlockwise(tri) ? +1 : -1;
+ if (orient2d(tri[0], tri[1], p1) * orient < 0 && orient2d(tri[0], tri[1], p2) * orient < 0) return false;
+ if (orient2d(tri[1], tri[2], p1) * orient < 0 && orient2d(tri[1], tri[2], p2) * orient < 0) return false;
+ if (orient2d(tri[2], tri[0], p1) * orient < 0 && orient2d(tri[2], tri[0], p2) * orient < 0) return false;
+ return true;
+ }
+
+ template<typename P2vec>
+ int triangleLineOrientation(const P2 &p1, const P2 &p2, const P2vec &tri) {
+ double lo, hi, tmp;
+ lo = hi = orient2d(p1, p2, tri[0]);
+ tmp = orient2d(p1, p2, tri[1]); lo = std::min(lo, tmp); hi = std::max(hi, tmp);
+ tmp = orient2d(p1, p2, tri[2]); lo = std::min(lo, tmp); hi = std::max(hi, tmp);
+ if (hi < 0.0) return -1;
+ if (lo > 0.0) return +1;
+ return 0;
+ }
+
+ template<typename P2vec>
+ bool triangleIntersectsTriangle(const P2vec &tri_b, const P2vec &tri_a) {
+ int orient_a = isAnticlockwise(tri_a) ? +1 : -1;
+ if (triangleLineOrientation(tri_a[0], tri_a[1], tri_b) * orient_a < 0) return false;
+ if (triangleLineOrientation(tri_a[1], tri_a[2], tri_b) * orient_a < 0) return false;
+ if (triangleLineOrientation(tri_a[2], tri_a[0], tri_b) * orient_a < 0) return false;
+
+ int orient_b = isAnticlockwise(tri_b) ? +1 : -1;
+ if (triangleLineOrientation(tri_b[0], tri_b[1], tri_a) * orient_b < 0) return false;
+ if (triangleLineOrientation(tri_b[1], tri_b[2], tri_a) * orient_b < 0) return false;
+ if (triangleLineOrientation(tri_b[2], tri_b[0], tri_a) * orient_b < 0) return false;
+
+ return true;
+ }
+
+
+
+ static inline double atan2(const P2 &p) {
+ return ::atan2(p.y, p.x);
+ }
+
+
+
+ struct LineIntersectionInfo {
+ LineIntersectionClass iclass;
+ P2 ipoint;
+ int p1, p2;
+
+ LineIntersectionInfo(LineIntersectionClass _iclass,
+ P2 _ipoint = P2::ZERO(),
+ int _p1 = -1,
+ int _p2 = -1) :
+ iclass(_iclass), ipoint(_ipoint), p1(_p1), p2(_p2) {
+ }
+ };
+
+ struct PolyInclusionInfo {
+ PointClass iclass;
+ int iobjnum;
+
+ PolyInclusionInfo(PointClass _iclass,
+ int _iobjnum = -1) :
+ iclass(_iclass), iobjnum(_iobjnum) {
+ }
+ };
+
+ struct PolyIntersectionInfo {
+ IntersectionClass iclass;
+ P2 ipoint;
+ size_t iobjnum;
+
+ PolyIntersectionInfo(IntersectionClass _iclass,
+ const P2 &_ipoint,
+ size_t _iobjnum) :
+ iclass(_iclass), ipoint(_ipoint), iobjnum(_iobjnum) {
+ }
+ };
+
+ bool lineSegmentIntersection_simple(const P2 &l1v1, const P2 &l1v2,
+ const P2 &l2v1, const P2 &l2v2);
+ bool lineSegmentIntersection_simple(const LineSegment2 &l1,
+ const LineSegment2 &l2);
+
+ LineIntersectionInfo lineSegmentIntersection(const P2 &l1v1, const P2 &l1v2,
+ const P2 &l2v1, const P2 &l2v2);
+ LineIntersectionInfo lineSegmentIntersection(const LineSegment2 &l1,
+ const LineSegment2 &l2);
+
+ int lineSegmentPolyIntersections(const std::vector<P2> &points,
+ LineSegment2 line,
+ std::vector<PolyInclusionInfo> &out);
+
+ int sortedLineSegmentPolyIntersections(const std::vector<P2> &points,
+ LineSegment2 line,
+ std::vector<PolyInclusionInfo> &out);
+
+
+
+ static inline bool quadIsConvex(const P2 &a, const P2 &b, const P2 &c, const P2 &d) {
+ double s_1, s_2;
+
+ s_1 = carve::geom2d::orient2d(a, c, b);
+ s_2 = carve::geom2d::orient2d(a, c, d);
+ if ((s_1 < 0.0 && s_2 < 0.0) || (s_1 > 0.0 && s_2 > 0.0)) return false;
+
+ s_1 = carve::geom2d::orient2d(b, d, a);
+ s_2 = carve::geom2d::orient2d(b, d, c);
+ if ((s_1 < 0.0 && s_2 < 0.0) || (s_1 > 0.0 && s_2 > 0.0)) return false;
+
+ return true;
+ }
+
+ template<typename T, typename adapt_t>
+ inline bool quadIsConvex(const T &a, const T &b, const T &c, const T &d, adapt_t adapt) {
+ return quadIsConvex(adapt(a), adapt(b), adapt(c), adapt(d));
+ }
+
+
+
+ double signedArea(const std::vector<P2> &points);
+
+ static inline double signedArea(const P2 &a, const P2 &b, const P2 &c) {
+ return ((b.y + a.y) * (b.x - a.x) + (c.y + b.y) * (c.x - b.x) + (a.y + c.y) * (a.x - c.x)) / 2.0;
+ }
+
+ template<typename T, typename adapt_t>
+ double signedArea(const std::vector<T> &points, adapt_t adapt) {
+ P2Vector::size_type l = points.size();
+ double A = 0.0;
+
+ for (P2Vector::size_type i = 0; i < l - 1; i++) {
+ A += (adapt(points[i + 1]).y + adapt(points[i]).y) * (adapt(points[i + 1]).x - adapt(points[i]).x);
+ }
+ A += (adapt(points[0]).y + adapt(points[l - 1]).y) * (adapt(points[0]).x - adapt(points[l - 1]).x);
+
+ return A / 2.0;
+ }
+
+
+
+ template<typename iter_t, typename adapt_t>
+ double signedArea(iter_t begin, iter_t end, adapt_t adapt) {
+ double A = 0.0;
+ P2 p, n;
+
+ if (begin == end) return 0.0;
+
+ p = adapt(*begin);
+ for (iter_t c = begin; ++c != end; ) {
+ P2 n = adapt(*c);
+ A += (n.y + p.y) * (n.x - p.x);
+ p = n;
+ }
+ n = adapt(*begin);
+ A += (n.y + p.y) * (n.x - p.x);
+
+ return A / 2.0;
+ }
+
+
+
+ bool pointInPolySimple(const std::vector<P2> &points, const P2 &p);
+
+ template<typename T, typename adapt_t>
+ bool pointInPolySimple(const std::vector<T> &points, adapt_t adapt, const P2 &p) {
+ CARVE_ASSERT(points.size() > 0);
+ P2Vector::size_type l = points.size();
+ double s = 0.0;
+ double rp, r0, d;
+
+ rp = r0 = atan2(adapt(points[0]) - p);
+
+ for (P2Vector::size_type i = 1; i < l; i++) {
+ double r = atan2(adapt(points[i]) - p);
+ d = r - rp;
+ if (d > M_PI) d -= M_TWOPI;
+ if (d < -M_PI) d += M_TWOPI;
+ s = s + d;
+ rp = r;
+ }
+
+ d = r0 - rp;
+ if (d > M_PI) d -= M_TWOPI;
+ if (d < -M_PI) d += M_TWOPI;
+ s = s + d;
+
+ return !carve::math::ZERO(s);
+ }
+
+
+
+ PolyInclusionInfo pointInPoly(const std::vector<P2> &points, const P2 &p);
+
+ template<typename T, typename adapt_t>
+ PolyInclusionInfo pointInPoly(const std::vector<T> &points, adapt_t adapt, const P2 &p) {
+ P2Vector::size_type l = points.size();
+ for (unsigned i = 0; i < l; i++) {
+ if (equal(adapt(points[i]), p)) return PolyInclusionInfo(POINT_VERTEX, i);
+ }
+
+ for (unsigned i = 0; i < l; i++) {
+ unsigned j = (i + 1) % l;
+
+ if (std::min(adapt(points[i]).x, adapt(points[j]).x) - EPSILON < p.x &&
+ std::max(adapt(points[i]).x, adapt(points[j]).x) + EPSILON > p.x &&
+ std::min(adapt(points[i]).y, adapt(points[j]).y) - EPSILON < p.y &&
+ std::max(adapt(points[i]).y, adapt(points[j]).y) + EPSILON > p.y &&
+ distance2(carve::geom::rayThrough(adapt(points[i]), adapt(points[j])), p) < EPSILON2) {
+ return PolyInclusionInfo(POINT_EDGE, i);
+ }
+ }
+
+ if (pointInPolySimple(points, adapt, p)) {
+ return PolyInclusionInfo(POINT_IN);
+ }
+
+ return PolyInclusionInfo(POINT_OUT);
+ }
+
+
+
+ bool pickContainedPoint(const std::vector<P2> &poly, P2 &result);
+
+ template<typename T, typename adapt_t>
+ bool pickContainedPoint(const std::vector<T> &poly, adapt_t adapt, P2 &result) {
+#if defined(CARVE_DEBUG)
+ std::cerr << "pickContainedPoint ";
+ for (unsigned i = 0; i < poly.size(); ++i) std::cerr << " " << adapt(poly[i]);
+ std::cerr << std::endl;
+#endif
+
+ const size_t S = poly.size();
+ P2 a, b, c;
+ for (unsigned i = 0; i < S; ++i) {
+ a = adapt(poly[i]);
+ b = adapt(poly[(i + 1) % S]);
+ c = adapt(poly[(i + 2) % S]);
+
+ if (cross(a - b, c - b) < 0) {
+ P2 p = (a + b + c) / 3;
+ if (pointInPolySimple(poly, adapt, p)) {
+ result = p;
+ return true;
+ }
+ }
+ }
+ return false;
+ }
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/geom.hpp>
+
+#include <math.h>
+#include <carve/math_constants.hpp>
+
+#include <vector>
+#include <list>
+#include <map>
+
+#if defined(CARVE_DEBUG)
+# include <iostream>
+#endif
+
+#if defined CARVE_USE_EXACT_PREDICATES
+# include <carve/shewchuk_predicates.hpp>
+#endif
+
+namespace carve {
+ namespace geom3d {
+
+ typedef carve::geom::plane<3> Plane;
+ typedef carve::geom::ray<3> Ray;
+ typedef carve::geom::linesegment<3> LineSegment;
+ typedef carve::geom::vector<3> Vector;
+
+ template<typename iter_t, typename adapt_t>
+ bool fitPlane(iter_t begin, iter_t end, adapt_t adapt, Plane &plane) {
+ Vector centroid;
+ carve::geom::centroid(begin, end, adapt, centroid);
+ iter_t i;
+
+ Vector n = Vector::ZERO();
+ Vector v;
+ Vector p1, p2, p3, c1, c2;
+ if (begin == end) return false;
+
+ i = begin;
+ p1 = c1 = adapt(*i); if (++i == end) return false;
+ p2 = c2 = adapt(*i); if (++i == end) return false;
+
+#if defined(CARVE_DEBUG)
+ size_t N = 2;
+#endif
+ do {
+ p3 = adapt(*i);
+ v = cross(p3 - p2, p1 - p2);
+ if (v.v[largestAxis(v)]) v.negate();
+ n += v;
+ p1 = p2; p2 = p3;
+#if defined(CARVE_DEBUG)
+ ++N;
+#endif
+ } while (++i != end);
+
+ p1 = p2; p2 = p3; p3 = c1;
+ v = cross(p3 - p2, p1 - p2);
+ if (v.v[largestAxis(v)]) v.negate();
+ n += v;
+
+ p1 = p2; p2 = p3; p3 = c2;
+ v = cross(p3 - p2, p1 - p2);
+ if (v.v[largestAxis(v)]) v.negate();
+ n += v;
+
+ n.normalize();
+ plane.N = n;
+ plane.d = -dot(n, centroid);
+#if defined(CARVE_DEBUG)
+ if (N > 3) {
+ std::cerr << "N = " << N << " fitted distance:";
+ for (i = begin; i != end; ++i) {
+ Vector p = adapt(*i);
+ std::cerr << " {" << p << "} " << distance(plane, p);
+ }
+ std::cerr << std::endl;
+ }
+#endif
+ return true;
+ }
+
+ bool planeIntersection(const Plane &a, const Plane &b, Ray &r);
+
+ IntersectionClass rayPlaneIntersection(const Plane &p,
+ const Vector &v1,
+ const Vector &v2,
+ Vector &v,
+ double &t);
+
+ IntersectionClass lineSegmentPlaneIntersection(const Plane &p,
+ const LineSegment &line,
+ Vector &v);
+
+ RayIntersectionClass rayRayIntersection(const Ray &r1,
+ const Ray &r2,
+ Vector &v1,
+ Vector &v2,
+ double &mu1,
+ double &mu2);
+
+
+
+ // test whether point d is above, below or on the plane formed by the triangle a,b,c.
+ // return: +ve = d is below a,b,c
+ // -ve = d is above a,b,c
+ // 0 = d is on a,b,c
+#if defined CARVE_USE_EXACT_PREDICATES
+ inline double orient3d(const Vector &a,
+ const Vector &b,
+ const Vector &c,
+ const Vector &d) {
+ return shewchuk::orient3d(a.v, b.v, c.v, d.v);
+ }
+#else
+ inline double orient3d(const Vector &a,
+ const Vector &b,
+ const Vector &c,
+ const Vector &d) {
+ return dotcross((a - d), (b - d), (c - d));
+ }
+#endif
+
+ // Volume of a tetrahedron described by 4 points. Will be
+ // positive if the anticlockwise normal of a,b,c is oriented out
+ // of the tetrahedron.
+ //
+ // see: http://mathworld.wolfram.com/Tetrahedron.html
+ inline double tetrahedronVolume(const Vector &a,
+ const Vector &b,
+ const Vector &c,
+ const Vector &d) {
+ return dotcross((a - d), (b - d), (c - d)) / 6.0;
+ }
+
+ /**
+ * \brief Determine whether p is internal to the wedge defined by
+ * the area between the planes defined by a,b,c and a,b,d
+ * angle abc, where ab is the apex of the angle.
+ *
+ * @param[in] a
+ * @param[in] b
+ * @param[in] c
+ * @param[in] d
+ * @param[in] p
+ *
+ * @return true, if p is contained in the wedge defined by the
+ * area between the planes defined by a,b,c and
+ * a,b,d. If the wedge is reflex, p is considered to
+ * be contained if it lies on either plane. Acute
+ * wdges do not contain p if p lies on either
+ * plane. This is so that internalToWedge(a,b,c,d,p) =
+ * !internalToWedge(a,b,d,c,p)
+ */
+ inline bool internalToWedge(const Vector &a,
+ const Vector &b,
+ const Vector &c,
+ const Vector &d,
+ const Vector &p) {
+ bool reflex = (c < d) ?
+ orient3d(a, b, c, d) >= 0.0 :
+ orient3d(a, b, d, c) < 0.0;
+
+ double d1 = orient3d(a, b, c, p);
+ double d2 = orient3d(a, b, d, p);
+
+ if (reflex) {
+ // above a,b,c or below a,b,d (or coplanar with either)
+ return d1 <= 0.0 || d2 >= 0.0;
+ } else {
+ // above a,b,c and below a,b,d
+ return d1 < 0.0 && d2 > 0.0;
+ }
+ }
+
+ /**
+ * \brief Determine the ordering relationship of a and b, when
+ * rotating around direction, starting from base.
+ *
+ * @param[in] adirection
+ * @param[in] base
+ * @param[in] a
+ * @param[in] b
+ *
+ * @return
+ * * -1, if a is ordered before b around, rotating about direction.
+ * * 0, if a and b are equal in angle.
+ * * +1, if a is ordered after b around, rotating about direction.
+ */
+ inline int compareAngles(const Vector &direction, const Vector &base, const Vector &a, const Vector &b) {
+ const double d1 = carve::geom3d::orient3d(carve::geom::VECTOR(0,0,0), direction, a, b);
+ const double d2 = carve::geom3d::orient3d(carve::geom::VECTOR(0,0,0), direction, base, a);
+ const double d3 = carve::geom3d::orient3d(carve::geom::VECTOR(0,0,0), direction, base, b);
+
+ // CASE: a and b are coplanar wrt. direction.
+ if (d1 == 0.0) {
+ // a and b point in the same direction.
+ if (dot(a, b) > 0.0) {
+ // Neither is less than the other.
+ return 0;
+ }
+
+ // a and b point in opposite directions.
+ // * if d2 < 0.0, a is above plane(direction, base) and is less
+ // than b.
+ // * if d2 == 0.0 a is coplanar with plane(direction, base) and is
+ // less than b if it points in the same direction as base.
+ // * if d2 > 0.0, a is below plane(direction, base) and is greater
+ // than b.
+
+ if (d2 == 0.0) { return dot(a, base) > 0.0 ? -1 : +1; }
+ if (d3 == 0.0) { return dot(b, base) > 0.0 ? +1 : -1; }
+ if (d2 < 0.0 && d3 > 0.0) return -1;
+ if (d2 > 0.0 && d3 < 0.0) return +1;
+
+ // both a and b are to one side of plane(direction, base) -
+ // rounding error (if a and b are truly coplanar with
+ // direction, one should be above, and one should be below any
+ // other plane that is not itself coplanar with
+ // plane(direction, a|b) - which would imply d2 and d3 == 0.0).
+
+ // If both are below plane(direction, base) then the one that
+ // points in the same direction as base is greater.
+ // If both are above plane(direction, base) then the one that
+ // points in the same direction as base is lesser.
+ if (d2 > 0.0) { return dot(a, base) > 0.0 ? +1 : -1; }
+ else { return dot(a, base) > 0.0 ? -1 : +1; }
+ }
+
+ // CASE: a and b are not coplanar wrt. direction
+
+ if (d2 < 0.0) {
+ // if a is above plane(direction,base), then a is less than b if
+ // b is below plane(direction,base) or b is above plane(direction,a)
+ return (d3 > 0.0 || d1 < 0.0) ? -1 : +1;
+ } else if (d2 == 0.0) {
+ // if a is on plane(direction,base) then a is less than b if a
+ // points in the same direction as base, or b is below
+ // plane(direction,base)
+ return (dot(a, base) > 0.0 || d3 > 0.0) ? -1 : +1;
+ } else {
+ // if a is below plane(direction,base), then a is less than b if b
+ // is below plane(direction,base) and b is above plane(direction,a)
+ return (d3 > 0.0 && d1 < 0.0) ? -1 : +1;
+ }
+ }
+
+ // The anticlockwise angle from vector "from" to vector "to", oriented around the vector "orient".
+ static inline double antiClockwiseAngle(const Vector &from, const Vector &to, const Vector &orient) {
+ double dp = dot(from, to);
+ Vector cp = cross(from, to);
+ if (cp.isZero()) {
+ if (dp < 0) {
+ return M_PI;
+ } else {
+ return 0.0;
+ }
+ } else {
+ if (dot(cp, orient) > 0.0) {
+ return acos(dp);
+ } else {
+ return M_TWOPI - acos(dp);
+ }
+ }
+ }
+
+
+
+ static inline double antiClockwiseOrdering(const Vector &from, const Vector &to, const Vector &orient) {
+ double dp = dot(from, to);
+ Vector cp = cross(from, to);
+ if (cp.isZero()) {
+ if (dp < 0) {
+ return 2.0;
+ } else {
+ return 0.0;
+ }
+ } else {
+ if (dot(cp, orient) > 0.0) {
+ // 1..-1 -> 0..2
+ return 1.0 - dp;
+ } else {
+ // -1..1 -> 2..4
+ return dp + 1.0;
+ }
+ }
+ }
+
+
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/math.hpp>
+
+namespace carve {
+ namespace geom {
+
+
+
+ template<unsigned ndim>
+ double vector<ndim>::length2() const { return dot(*this, *this); }
+ template<unsigned ndim>
+ double vector<ndim>::length() const { return sqrt(dot(*this, *this)); }
+
+ template<unsigned ndim>
+ vector<ndim> &vector<ndim>::normalize() { *this /= length(); return *this; }
+ template<unsigned ndim>
+ vector<ndim> vector<ndim>::normalized() const { return *this / length(); }
+
+ template<unsigned ndim>
+ bool vector<ndim>::exactlyZero() const {
+ for (unsigned i = 0; i < ndim; ++i) if (this->v[i]) return false;
+ return true;
+ }
+ template<unsigned ndim>
+ bool vector<ndim>::isZero(double epsilon) const {
+ return length2() < epsilon * epsilon;
+ }
+
+ template<unsigned ndim>
+ void vector<ndim>::setZero() { for (size_t i = 0; i < ndim; ++i) this->v[i] = 0.0; }
+
+ template<unsigned ndim>
+ void vector<ndim>::fill(double val) { for (size_t i = 0; i < ndim; ++i) this->v[i] = val; }
+
+ template<unsigned ndim>
+ vector<ndim> &vector<ndim>::scaleBy(double d) { for (unsigned i = 0; i < ndim; ++i) this->v[i] *= d; return *this; }
+ template<unsigned ndim>
+ vector<ndim> &vector<ndim>::invscaleBy(double d) { for (unsigned i = 0; i < ndim; ++i) this->v[i] /= d; return *this; }
+
+ template<unsigned ndim>
+ vector<ndim> vector<ndim>::scaled(double d) const { return *this * d; }
+ template<unsigned ndim>
+ vector<ndim> vector<ndim>::invscaled(double d) const { return *this / d; }
+
+ template<unsigned ndim>
+ vector<ndim> &vector<ndim>::negate() { for (unsigned i = 0; i < ndim; ++i) this->v[i] = -this->v[i]; return *this; }
+ template<unsigned ndim>
+ vector<ndim> vector<ndim>::negated() const { return -*this; }
+
+ template<unsigned ndim>
+ double &vector<ndim>::operator[](unsigned i) { return this->v[i]; }
+ template<unsigned ndim>
+ const double &vector<ndim>::operator[](unsigned i) const { return this->v[i]; }
+
+ template<unsigned ndim>
+ template<typename assign_t>
+ vector<ndim> &vector<ndim>::operator=(const assign_t &t) {
+ for (unsigned i = 0; i < ndim; ++i) this->v[i] = t[i];
+ return *this;
+ }
+
+ template<unsigned ndim>
+ std::string vector<ndim>::asStr() const {
+ std::ostringstream out;
+ out << '<';
+ out << std::setprecision(24);
+ for (unsigned i = 0; i < ndim; ++i) { if (i) out << ','; out << this->v[i]; }
+ out << '>';
+ return out.str();
+ }
+
+
+
+ template<unsigned ndim>
+ vector<ndim> operator+(const vector<ndim> &a, const vector<ndim> &b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] + b[i];
+ return c;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> operator+(const vector<ndim> &a, double b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] + b;
+ return c;
+ }
+
+ template<unsigned ndim, typename val_t>
+ vector<ndim> operator+(const vector<ndim> &a, const val_t &b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] + b[i];
+ return c;
+ }
+
+ template<unsigned ndim, typename val_t>
+ vector<ndim> operator+(const val_t &a, const vector<ndim> &b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] + b[i];
+ return c;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> &operator+=(vector<ndim> &a, const vector<ndim> &b) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] += b[i];
+ return a;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> &operator+=(vector<ndim> &a, double b) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] += b;
+ return a;
+ }
+
+ template<unsigned ndim, typename val_t>
+ vector<ndim> &operator+=(vector<ndim> &a, const val_t &b) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] += b[i];
+ return a;
+ }
+
+
+
+ template<unsigned ndim>
+ vector<ndim> operator-(const vector<ndim> &a) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = -a[i];
+ return c;
+ }
+
+
+
+ template<unsigned ndim>
+ vector<ndim> operator-(const vector<ndim> &a, double b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] - b;
+ return c;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> operator-(const vector<ndim> &a, const vector<ndim> &b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] - b[i];
+ return c;
+ }
+
+ template<unsigned ndim, typename val_t>
+ vector<ndim> operator-(const vector<ndim> &a, const val_t &b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] - b[i];
+ return c;
+ }
+
+ template<unsigned ndim, typename val_t>
+ vector<ndim> operator-(const val_t &a, const vector<ndim> &b) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] - b[i];
+ return c;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> &operator-=(vector<ndim> &a, const vector<ndim> &b) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] -= b[i];
+ return a;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> &operator-=(vector<ndim> &a, double b) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] -= b;
+ return a;
+ }
+
+ template<unsigned ndim, typename val_t>
+ vector<ndim> &operator-=(vector<ndim> &a, const val_t &b) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] -= b[i];
+ return a;
+ }
+
+
+
+ template<unsigned ndim>
+ vector<ndim> operator*(const vector<ndim> &a, double s) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] * s;
+ return c;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> operator*(double s, const vector<ndim> &a) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] * s;
+ return c;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> &operator*=(vector<ndim> &a, double s) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] *= s;
+ return a;
+ }
+
+
+
+ template<unsigned ndim>
+ vector<ndim> operator/(const vector<ndim> &a, double s) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = a[i] / s;
+ return c;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> &operator/=(vector<ndim> &a, double s) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] /= s;
+ return a;
+ }
+
+
+
+ template<unsigned ndim>
+ bool operator==(const vector<ndim> &a, const vector<ndim> &b) {
+ for (unsigned i = 0; i < ndim; ++i) { if (a[i] != b[i]) return false; }
+ return true;
+ }
+
+ template<unsigned ndim>
+ bool operator!=(const vector<ndim> &a, const vector<ndim> &b) {
+ return !(a == b);
+ }
+
+ template<unsigned ndim>
+ bool operator<(const vector<ndim> &a, const vector<ndim> &b) {
+ for (unsigned i = 0; i < ndim; ++i) { if (a[i] < b[i]) return true; if (a[i] > b[i]) return false; }
+ return false;
+ }
+
+ template<unsigned ndim>
+ bool operator<=(const vector<ndim> &a, const vector<ndim> &b) {
+ return !(b < a);
+ }
+
+ template<unsigned ndim>
+ bool operator>(const vector<ndim> &a, const vector<ndim> &b) {
+ return b < a;
+ }
+
+ template<unsigned ndim>
+ bool operator>=(const vector<ndim> &a, const vector<ndim> &b) {
+ return !(a < b);
+ }
+
+
+
+ template<unsigned ndim>
+ vector<ndim> abs(const vector<ndim> &a) {
+ vector<ndim> c(NOINIT);
+ for (unsigned i = 0; i < ndim; ++i) c[i] = fabs(a[i]);
+ return c;
+ }
+
+ template<unsigned ndim>
+ double distance2(const vector<ndim> &a, const vector<ndim> &b) {
+ return (b - a).length2();
+ }
+
+ template<unsigned ndim>
+ double distance(const vector<ndim> &a, const vector<ndim> &b) {
+ return (b - a).length();
+ }
+
+ template<unsigned ndim>
+ bool equal(const vector<ndim> &a, const vector<ndim> &b) {
+ return (b - a).isZero();
+ }
+
+ template<unsigned ndim>
+ int smallestAxis(const vector<ndim> &a) {
+ int x = 0;
+ double y = fabs(a[0]);
+ for (unsigned i = 1; i < ndim; ++i) {
+ double z = fabs(a[i]);
+ if (z <= y) { y = z; x = i; }
+ }
+ return x;
+ }
+
+ template<unsigned ndim>
+ int largestAxis(const vector<ndim> &a) {
+ int x = 0;
+ double y = fabs(a[0]);
+ for (unsigned i = 1; i < ndim; ++i) {
+ double z = fabs(a[i]);
+ if (z > y) { y = z; x = i; }
+ }
+ return x;
+ }
+
+ template<unsigned ndim>
+ vector<2> select(const vector<ndim> &a, int a1, int a2) {
+ vector<2> r(NOINIT);
+ r.v[0] = a.v[a1]; r.v[1] = a.v[a2];
+ return r;
+ }
+
+ template<unsigned ndim>
+ vector<3> select(const vector<ndim> &a, int a1, int a2, int a3) {
+ vector<3> r(NOINIT);
+ r.v[0] = a.v[a1]; r.v[1] = a.v[a2]; r.v[2] = a.v[a3];
+ return r;
+ }
+
+ template<unsigned ndim, typename assign_t, typename oper_t>
+ vector<ndim> &assign_op(vector<ndim> &a, const assign_t &t, oper_t op) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] = op(t[i]);
+ return a;
+ }
+
+ template<unsigned ndim, typename assign1_t, typename assign2_t, typename oper_t>
+ vector<ndim> &assign_op(vector<ndim> &a, const assign1_t &t1, const assign2_t &t2, oper_t op) {
+ for (unsigned i = 0; i < ndim; ++i) a[i] = op(t1[i], t2[i]);
+ return a;
+ }
+
+ template<unsigned ndim, typename iter_t>
+ void bounds(iter_t begin, iter_t end, vector<ndim> &min, vector<ndim> &max) {
+ if (begin == end) {
+ min.setZero();
+ max.setZero();
+ } else {
+ min = max = *begin;
+ while (++begin != end) {
+ vector<ndim> v = *begin;
+ assign_op(min, min, v, carve::util::min_functor());
+ assign_op(max, max, v, carve::util::max_functor());
+ }
+ }
+ }
+
+ template<unsigned ndim, typename iter_t, typename adapt_t>
+ void bounds(iter_t begin, iter_t end, adapt_t adapt, vector<ndim> &min, vector<ndim> &max) {
+ if (begin == end) {
+ min.setZero();
+ max.setZero();
+ } else {
+ min = max = adapt(*begin);
+ while (++begin != end) {
+ vector<ndim> v = adapt(*begin);
+ assign_op(min, min, v, carve::util::min_functor());
+ assign_op(max, max, v, carve::util::max_functor());
+ }
+ }
+ }
+
+ template<unsigned ndim, typename iter_t, typename adapt_t>
+ void centroid(iter_t begin, iter_t end, adapt_t adapt, vector<ndim> &c) {
+ c.setZero();
+ int n = 0;
+ while (begin != end) { c += adapt(*begin++); ++n; }
+ c /= double(n);
+ }
+
+ template<unsigned ndim, typename val_t>
+ double dot(const vector<ndim> &a, const val_t &b) {
+ double r = 0.0;
+ for (unsigned i = 0; i < ndim; ++i) r += a[i] * b[i];
+ return r;
+ }
+
+ static inline vector<3> cross(const vector<3> &a, const vector<3> &b) {
+ // Compute a x b
+ return VECTOR(+(a.y * b.z - a.z * b.y),
+ -(a.x * b.z - a.z * b.x),
+ +(a.x * b.y - a.y * b.x));
+ }
+
+ static inline double cross(const vector<2> &a, const vector<2> &b) {
+ // Compute a x b
+ return a.x * b.y - b.x * a.y;
+ }
+
+ static inline double dotcross(const vector<3> &a, const vector<3> &b, const vector<3> &c) {
+ // Compute a . (b x c)
+ return
+ (a.x * b.y * c.z + a.y * b.z * c.x + a.z * b.x * c.y) -
+ (a.x * b.z * c.y + a.y * b.x * c.z + a.z * b.y * c.x);
+ }
+
+
+
+ template<unsigned ndim>
+ double distance(const axis_pos &a, const vector<ndim> &b) {
+ return fabs(b[a.axis] - a.pos);
+ }
+
+ template<unsigned ndim>
+ double distance2(const axis_pos &a, const vector<ndim> &b) {
+ double r = fabs(b[a.axis] - a.pos);
+ return r * r;
+ }
+
+ template<unsigned ndim> bool operator<(const axis_pos &a, const vector<ndim> &b) { return a.pos < b[a.axis]; }
+ template<unsigned ndim> bool operator<(const vector<ndim> &a, const axis_pos &b) { return a[b.axis] < b.pos; }
+
+ template<unsigned ndim> bool operator<=(const axis_pos &a, const vector<ndim> &b) { return a.pos <= b[a.axis]; }
+ template<unsigned ndim> bool operator<=(const vector<ndim> &a, const axis_pos &b) { return a[b.axis] <= b.pos; }
+
+ template<unsigned ndim> bool operator>(const axis_pos &a, const vector<ndim> &b) { return a.pos > b[a.axis]; }
+ template<unsigned ndim> bool operator>(const vector<ndim> &a, const axis_pos &b) { return a[b.axis] > b.pos; }
+
+ template<unsigned ndim> bool operator>=(const axis_pos &a, const vector<ndim> &b) { return a.pos >= b[a.axis]; }
+ template<unsigned ndim> bool operator>=(const vector<ndim> &a, const axis_pos &b) { return a[b.axis] >= b.pos; }
+
+ template<unsigned ndim> bool operator==(const axis_pos &a, const vector<ndim> &b) { return a.pos == b[a.axis]; }
+ template<unsigned ndim> bool operator==(const vector<ndim> &a, const axis_pos &b) { return a[b.axis] == b.pos; }
+
+ template<unsigned ndim> bool operator!=(const axis_pos &a, const vector<ndim> &b) { return a.pos != b[a.axis]; }
+ template<unsigned ndim> bool operator!=(const vector<ndim> &a, const axis_pos &b) { return a[b.axis] != b.pos; }
+
+
+
+ template<unsigned ndim>
+ bool ray<ndim>::OK() const {
+ return !D.isZero();
+ }
+
+ template<unsigned ndim>
+ ray<ndim> rayThrough(const vector<ndim> &a, const vector<ndim> &b) {
+ return ray<ndim>(b - a, a);
+ }
+
+ static inline double distance2(const ray<3> &r, const vector<3> &v) {
+ return cross(r.D, v - r.v).length2() / r.D.length2();
+ }
+
+ static inline double distance(const ray<3> &r, const vector<3> &v) {
+ return sqrt(distance2(r, v));
+ }
+
+ static inline double distance2(const ray<2> &r, const vector<2> &v) {
+ double t = cross(r.D, v - r.v);
+ return (t * t) / r.D.length2();
+ }
+
+ static inline double distance(const ray<2> &r, const vector<2> &v) {
+ return sqrt(distance2(r, v));
+ }
+
+
+
+ template<unsigned ndim>
+ void linesegment<ndim>::update() {
+ midpoint = (v2 + v1) / 2.0;
+ half_length = (v2 - v1) / 2.0;
+ }
+
+ template<unsigned ndim>
+ bool linesegment<ndim>::OK() const {
+ return !half_length.isZero();
+ }
+
+ template<unsigned ndim>
+ void linesegment<ndim>::flip() {
+ std::swap(v1, v2);
+ half_length = (v2 - v1) / 2.0;
+ }
+
+ template<unsigned ndim>
+ aabb<ndim> linesegment<ndim>::getAABB() const {
+ aabb<ndim> r;
+ r.fit(v1, v2);
+ return r;
+ }
+
+ template<unsigned ndim>
+ linesegment<ndim>::linesegment(const vector_t &_v1, const vector_t &_v2) : v1(_v1), v2(_v2) {
+ update();
+ }
+
+
+
+ template<unsigned ndim>
+ double distance2(const linesegment<ndim> &l, const vector<ndim> &v) {
+ vector<ndim> D = l.v2 - l.v1;
+ double t = dot(v - l.v1, D) / dot(D, D);
+ if (t <= 0.0) return (v - l.v1).length2();
+ if (t >= 1.0) return (v - l.v2).length2();
+ vector<ndim> vc = D * t + l.v1;
+ return (v - vc).length2();
+ }
+
+ template<unsigned ndim>
+ double distance(const linesegment<ndim> &l, const vector<ndim> &v) {
+ return sqrt(distance2(l, v));
+ }
+
+
+ template<unsigned ndim>
+ void plane<ndim>::negate() {
+ N.negate();
+ d = -d;
+ }
+
+ template<unsigned ndim>
+ plane<ndim>::plane() {
+ N.setZero();
+ N[0] = 1.0;
+ d= 0.0;
+ }
+
+ template<unsigned ndim>
+ plane<ndim>::plane(const vector_t &_N, vector_t _p) : N(_N), d(-dot(_p, _N)) {
+ }
+
+ template<unsigned ndim>
+ plane<ndim>::plane(const vector_t &_N, double _d) : N(_N), d(_d) {
+ }
+
+
+
+ template<unsigned ndim>
+ plane<ndim> operator-(const plane<ndim> &p) {
+ return plane<ndim>(-p.N, -p.d);
+ }
+
+ template<unsigned ndim, typename val_t>
+ double distance(const plane<ndim> &plane, const val_t &point) {
+ return dot(plane.N, point) + plane.d;
+ }
+
+ template<unsigned ndim, typename val_t>
+ double distance2(const plane<ndim> &plane, const val_t &point) {
+ double d = distance(plane, point);
+ return d * d;
+ }
+
+ template<unsigned ndim>
+ vector<ndim> closestPoint(const plane<ndim> &p, const vector<ndim> &v) {
+ return v - p.N * (p.d + dot(p.N, v)) / dot(p.N, p.N);
+ }
+
+
+
+ template<unsigned ndim>
+ aabb<ndim> sphere<ndim>::getAABB() const {
+ aabb<ndim> r;
+ r.fit(C - r, C + r);
+ }
+
+ template<unsigned ndim>
+ sphere<ndim>::sphere() {
+ C.setZero();
+ r = 1.0;
+ }
+
+ template<unsigned ndim>
+ sphere<ndim>::sphere(const vector_t &_C, double _r) : C(_C), r(_r) {
+ }
+
+
+
+ template<unsigned ndim, typename val_t>
+ double distance(const sphere<ndim> &sphere, const val_t &point) {
+ return std::max(0.0, distance(sphere.C, point) - sphere.r);
+ }
+
+ template<unsigned ndim, typename val_t>
+ double distance2(const sphere<ndim> &sphere, const val_t &point) {
+ return std::max(0.0, distance2(sphere.C, point) - sphere.r * sphere.r);
+ }
+
+ template<unsigned ndim>
+ vector<ndim> closestPoint(const sphere<ndim> &sphere, const vector<ndim> &point) {
+ return (point - sphere.C).normalized() * sphere.r;
+ }
+
+
+
+ template<unsigned ndim>
+ aabb<ndim> tri<ndim>::getAABB() const {
+ aabb<ndim> aabb;
+ aabb.fit(v[0], v[1], v[2]);
+ return aabb;
+ }
+
+ template<unsigned ndim>
+ tri<ndim>::tri(vector_t _v[3]) {
+ std::copy(v, v+3, _v);
+ }
+
+ template<unsigned ndim>
+ tri<ndim>::tri(const vector_t &a, const vector_t &b, const vector_t &c) {
+ v[0] = a;
+ v[1] = b;
+ v[2] = c;
+ }
+
+
+
+ template<unsigned ndim>
+ std::ostream &operator<<(std::ostream &o, const vector<ndim> &v) {
+ o << v.asStr();
+ return o;
+ }
+
+ template<unsigned ndim>
+ std::ostream &operator<<(std::ostream &o, const carve::geom::plane<ndim> &p) {
+ o << p.N << ";" << p.d;
+ return o;
+ }
+
+ template<unsigned ndim>
+ std::ostream &operator<<(std::ostream &o, const carve::geom::sphere<ndim> &sphere) {
+ o << "{sphere " << sphere.C << ";" << sphere.r << "}";
+ return o;
+ }
+
+ template<unsigned ndim>
+ std::ostream &operator<<(std::ostream &o, const carve::geom::tri<ndim> &tri) {
+ o << "{tri " << tri.v[0] << ";" << tri.v[1] << ";" << tri.v[2] << "}";
+ return o;
+ }
+
+
+
+ template<unsigned ndim>
+ double distance(const tri<ndim> &tri, const vector<ndim> &pt) {
+ return distance(closestPoint(tri, pt), pt);
+ }
+
+
+
+ template<unsigned ndim>
+ double distance2(const tri<ndim> &tri, const vector<ndim> &pt) {
+ return distance2(closestPoint(tri, pt), pt);
+ }
+ }
+}
--- /dev/null
+// Copyright 2006 Tobias Sargeant (toby@permuted.net)
+// All rights reserved.
+
+#pragma once
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+#pragma once
+namespace carve {
+ namespace heap {
+ namespace detail {
+
+
+
+ struct ignore_position_t {
+ template<typename value_t>
+ void operator()(value_t &val, size_t idx) const {}
+ };
+
+
+
+ template<typename random_access_iter_t,
+ typename distance_t,
+ typename value_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void _adjust_heap(random_access_iter_t begin,
+ distance_t pos,
+ distance_t len,
+ value_t val,
+ pred_t pred,
+ pos_notifier_t notify) {
+ const distance_t top = pos;
+
+ distance_t child = pos * 2 + 2;
+ while (child < len) {
+ if (pred(begin[child], begin[child - 1])) child--;
+
+ begin[pos] = begin[child];
+ notify(begin[pos], pos);
+ pos = child;
+ child = pos * 2 + 2;
+ }
+
+ if (child == len) {
+ child--;
+ begin[pos] = begin[child];
+ notify(begin[pos], pos);
+ pos = child;
+ }
+
+ distance_t parent = (pos - 1) / 2;
+ while (pos > top && pred(begin[parent], val)) {
+ begin[pos] = begin[parent];
+ notify(begin[pos], pos);
+ pos = parent;
+ parent = (pos - 1) / 2;
+ }
+
+ begin[pos] = val;
+ notify(begin[pos], pos);
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename distance_t,
+ typename value_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void _push_heap(random_access_iter_t begin,
+ distance_t pos,
+ value_t val,
+ pred_t pred,
+ pos_notifier_t notify) {
+ distance_t parent = (pos - 1) / 2;
+ while (pos > 0 && pred(begin[parent], val)) {
+ begin[pos] = begin[parent];
+ notify(begin[pos], pos);
+ pos = parent;
+ parent = (pos - 1) / 2;
+ }
+ begin[pos] = val;
+ notify(begin[pos], pos);
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename distance_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void _remove_heap(random_access_iter_t begin,
+ distance_t pos,
+ distance_t len,
+ pred_t pred,
+ pos_notifier_t notify) {
+ --len;
+ if (pos != len) {
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+ value_t removed = begin[pos];
+ _adjust_heap(begin, pos, len, begin[len], pred, notify);
+ begin[len] = removed;
+ notify(begin[len], len);
+ }
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename distance_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void _make_heap(random_access_iter_t begin,
+ distance_t len,
+ pred_t pred,
+ pos_notifier_t notify) {
+ for (distance_t pos = len / 2; pos > 0; ) {
+ --pos;
+ _adjust_heap(begin, pos, len, begin[pos], pred, ignore_position_t());
+ }
+ for (distance_t pos = 0; pos < len; ++pos) {
+ notify(begin[pos], pos);
+ }
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename distance_t,
+ typename pred_t>
+ void _make_heap(random_access_iter_t begin,
+ distance_t len,
+ pred_t pred,
+ ignore_position_t) {
+ for (distance_t pos = len / 2; pos > 0; ) {
+ --pos;
+ _adjust_heap(begin, pos, len, begin[pos], pred, ignore_position_t());
+ }
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename distance_t,
+ typename pred_t>
+ bool _is_heap(random_access_iter_t begin,
+ distance_t len,
+ pred_t pred) {
+ distance_t parent = 0;
+
+ for (distance_t child = 1; child < len; ++child) {
+ if (pred(begin[parent], begin[child])) {
+ return false;
+ }
+ if (++child == len) break;
+ if (pred(begin[parent], begin[child])) {
+ return false;
+ }
+ ++parent;
+ }
+
+ return true;
+ }
+
+
+
+ }
+
+
+
+ template<typename random_access_iter_t>
+ void adjust_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ random_access_iter_t pos) {
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+
+ detail::_adjust_heap(begin, pos - begin, end - begin, *pos, std::less<value_t>());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t>
+ void adjust_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ random_access_iter_t pos,
+ pred_t pred) {
+ detail::_adjust_heap(begin, pos - begin, end - begin, *pos, pred);
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void adjust_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ random_access_iter_t pos,
+ pred_t pred,
+ pos_notifier_t notify) {
+ detail::_adjust_heap(begin, pos - begin, end - begin, *pos, pred, notify);
+ }
+
+
+
+ template<typename random_access_iter_t>
+ void remove_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ random_access_iter_t pos) {
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+
+ detail::_remove_heap(begin, pos - begin, end - begin, std::less<value_t>(), detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t>
+ void remove_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ random_access_iter_t pos,
+ pred_t pred) {
+ detail::_remove_heap(begin, pos - begin, end - begin, pred, detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void remove_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ random_access_iter_t pos,
+ pred_t pred,
+ pos_notifier_t notify) {
+ detail::_remove_heap(begin, pos - begin, end - begin, pred, notify);
+ }
+
+
+
+ template<typename random_access_iter_t>
+ void pop_heap(random_access_iter_t begin,
+ random_access_iter_t end) {
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ detail::_remove_heap(begin, distance_t(0), end - begin, std::less<value_t>(), detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t>
+ void pop_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred) {
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ detail::_remove_heap(begin, distance_t(0), end - begin, pred, detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void pop_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred,
+ pos_notifier_t notify) {
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ detail::_remove_heap(begin, distance_t(0), end - begin, pred, notify);
+ }
+
+
+
+ template<typename random_access_iter_t>
+ void push_heap(random_access_iter_t begin,
+ random_access_iter_t end) {
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ distance_t pos = end - begin - 1;
+ detail::_push_heap(begin, pos, begin[pos], std::less<value_t>(), detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t>
+ void push_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred) {
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ distance_t pos = end - begin - 1;
+ detail::_push_heap(begin, pos, begin[pos], pred, detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void push_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred,
+ pos_notifier_t notify) {
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ distance_t pos = end - begin - 1;
+ detail::_push_heap(begin, pos, begin[pos], pred, notify);
+ }
+
+
+
+ template<typename random_access_iter_t>
+ void make_heap(random_access_iter_t begin,
+ random_access_iter_t end) {
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+
+ detail::_make_heap(begin, end - begin, std::less<value_t>(), detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t>
+ void make_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred) {
+ detail::_make_heap(begin, end - begin, pred, detail::ignore_position_t());
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void make_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred,
+ pos_notifier_t notify) {
+ detail::_make_heap(begin, end - begin, pred, notify);
+ }
+
+
+
+ template<typename random_access_iter_t>
+ bool is_heap(random_access_iter_t begin,
+ random_access_iter_t end) {
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+
+ return detail::_is_heap(begin, end - begin, std::less<value_t>());
+ }
+
+
+
+ template<typename random_access_iter_t, typename pred_t>
+ bool is_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred) {
+ return detail::_is_heap(begin, end - begin, pred);
+ }
+
+
+
+ template<typename random_access_iter_t>
+ void sort_heap(random_access_iter_t begin,
+ random_access_iter_t end) {
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+ typedef typename std::iterator_traits<random_access_iter_t>::value_type value_t;
+
+ for (distance_t len = end - begin; len > 1; --len) {
+ detail::_remove_heap(begin, distance_t(0), len, std::less<value_t>(), detail::ignore_position_t());
+ }
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t>
+ void sort_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred) {
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ for (distance_t len = end - begin; len > 1; --len) {
+ detail::_remove_heap(begin, distance_t(0), len, pred, detail::ignore_position_t());
+ }
+ }
+
+
+
+ template<typename random_access_iter_t,
+ typename pred_t,
+ typename pos_notifier_t>
+ void sort_heap(random_access_iter_t begin,
+ random_access_iter_t end,
+ pred_t pred,
+ pos_notifier_t notify) {
+ typedef typename std::iterator_traits<random_access_iter_t>::difference_type distance_t;
+
+ for (distance_t len = end - begin; len > 1; --len) {
+ detail::_remove_heap(begin, distance_t(0), len, pred, detail::ignore_position_t());
+ notify(begin[len], len);
+ }
+ notify(begin[0], 0);
+ }
+
+
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/poly.hpp>
+#include <carve/mesh.hpp>
+#include <carve/polyline.hpp>
+#include <carve/pointset.hpp>
+
+
+
+namespace carve {
+ namespace input {
+
+ struct Data {
+ Data() {
+ }
+
+ virtual ~Data() {
+ }
+
+ virtual void transform(const carve::math::Matrix & /* transform */) {
+ }
+ };
+
+
+
+ struct VertexData : public Data {
+ std::vector<carve::geom3d::Vector> points;
+
+ VertexData() : Data() {
+ }
+
+ virtual ~VertexData() {
+ }
+
+ virtual void transform(const carve::math::Matrix &transform) {
+ for (size_t i = 0; i < points.size(); ++i) {
+ points[i] *= transform;
+ }
+ }
+
+ size_t addVertex(carve::geom3d::Vector point) {
+ size_t index = points.size();
+ points.push_back(point);
+ return index;
+ }
+
+ inline void reserveVertices(int count) {
+ points.reserve(count);
+ }
+
+ size_t getVertexCount() const {
+ return points.size();
+ }
+
+ const carve::geom3d::Vector &getVertex(int index) const {
+ return points[index];
+ }
+ };
+
+
+
+ struct PolyhedronData : public VertexData {
+ std::vector<int> faceIndices;
+ int faceCount;
+
+ PolyhedronData() : VertexData(), faceIndices(), faceCount(0) {
+ }
+
+ virtual ~PolyhedronData() {
+ }
+
+ void reserveFaces(int count, int avgFaceSize) {
+ faceIndices.reserve(faceIndices.size() + count * (1 + avgFaceSize));
+ }
+
+ int getFaceCount() const {
+ return faceCount;
+ }
+
+ template <typename Iter>
+ void addFace(Iter begin, Iter end) {
+ size_t n = std::distance(begin, end);
+ faceIndices.reserve(faceIndices.size() + n + 1);
+ faceIndices.push_back(n);
+ std::copy(begin, end, std::back_inserter(faceIndices));
+ ++faceCount;
+ }
+
+ void addFace(int a, int b, int c) {
+ faceIndices.push_back(3);
+ faceIndices.push_back(a);
+ faceIndices.push_back(b);
+ faceIndices.push_back(c);
+ ++faceCount;
+ }
+
+ void addFace(int a, int b, int c, int d) {
+ faceIndices.push_back(4);
+ faceIndices.push_back(a);
+ faceIndices.push_back(b);
+ faceIndices.push_back(c);
+ faceIndices.push_back(d);
+ ++faceCount;
+ }
+
+ void clearFaces() {
+ faceIndices.clear();
+ faceCount = 0;
+ }
+
+ carve::poly::Polyhedron *create() const {
+ return new carve::poly::Polyhedron(points, faceCount, faceIndices);
+ }
+
+ carve::mesh::MeshSet<3> *createMesh() const {
+ return new carve::mesh::MeshSet<3>(points, faceCount, faceIndices);
+ }
+ };
+
+
+
+ struct PolylineSetData : public VertexData {
+ typedef std::pair<bool, std::vector<int> > polyline_data_t;
+ std::list<polyline_data_t> polylines;
+
+ PolylineSetData() : VertexData(), polylines() {
+ }
+
+ virtual ~PolylineSetData() {
+ }
+
+ void beginPolyline(bool closed = false) {
+ polylines.push_back(std::make_pair(closed, std::vector<int>()));
+ }
+
+ void reservePolyline(size_t len) {
+ polylines.back().second.reserve(len);
+ }
+
+ void addPolylineIndex(int idx) {
+ polylines.back().second.push_back(idx);
+ }
+
+ carve::line::PolylineSet *create() const {
+ carve::line::PolylineSet *p = new carve::line::PolylineSet(points);
+
+ for (std::list<polyline_data_t>::const_iterator i = polylines.begin();
+ i != polylines.end();
+ ++i) {
+ p->addPolyline((*i).first, (*i).second.begin(), (*i).second.end());
+ }
+ return p;
+ }
+ };
+
+
+
+ struct PointSetData : public VertexData {
+
+ PointSetData() : VertexData() {
+ }
+
+ virtual ~PointSetData() {
+ }
+
+ carve::point::PointSet *create() const {
+ carve::point::PointSet *p = new carve::point::PointSet(points);
+ return p;
+ }
+ };
+
+
+
+ class Input {
+ public:
+ std::list<Data *> input;
+
+ Input() {
+ }
+
+ ~Input() {
+ for (std::list<Data *>::iterator i = input.begin(); i != input.end(); ++i) {
+ delete (*i);
+ }
+ }
+
+ void addDataBlock(Data *data) {
+ input.push_back(data);
+ }
+
+ void transform(const carve::math::Matrix &transform) {
+ if (transform == carve::math::Matrix::IDENT()) return;
+ for (std::list<Data *>::iterator i = input.begin(); i != input.end(); ++i) {
+ (*i)->transform(transform);
+ }
+ }
+
+ template<typename T>
+ static inline T *create(Data *d) {
+ return NULL;
+ }
+ };
+
+ template<>
+ inline carve::mesh::MeshSet<3> *Input::create(Data *d) {
+ PolyhedronData *p = dynamic_cast<PolyhedronData *>(d);
+ if (p == NULL) return NULL;
+ return p->createMesh();
+ }
+
+ template<>
+ inline carve::poly::Polyhedron *Input::create(Data *d) {
+ PolyhedronData *p = dynamic_cast<PolyhedronData *>(d);
+ if (p == NULL) return NULL;
+ return p->create();
+ }
+
+ template<>
+ inline carve::line::PolylineSet *Input::create(Data *d) {
+ PolylineSetData *p = dynamic_cast<PolylineSetData *>(d);
+ if (p == NULL) return NULL;
+ return p->create();
+ }
+
+ template<>
+ inline carve::point::PointSet *Input::create(Data *d) {
+ PointSetData *p = dynamic_cast<PointSetData *>(d);
+ if (p == NULL) return NULL;
+ return p->create();
+ }
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/geom2d.hpp>
+#include <carve/poly.hpp>
+#include <carve/mesh.hpp>
+#include <carve/csg.hpp>
+
+namespace carve {
+ namespace interpolate {
+
+ static inline std::vector<double> polyInterpolate(const std::vector<carve::geom2d::P2> &s,
+ const carve::geom2d::P2 &v) {
+ // see hormann et al. 2006
+ const size_t SZ = s.size();
+ std::vector<double> r;
+ std::vector<double> A;
+ std::vector<double> D;
+
+ std::vector<double> result;
+
+ r.resize(SZ);
+ A.resize(SZ);
+ D.resize(SZ);
+
+ result.resize(SZ, 0.0);
+
+ for (size_t i = 0; i < SZ; ++i) {
+ size_t i2 = (i + 1) % SZ;
+ carve::geom2d::P2 si = s[i] - v;
+ carve::geom2d::P2 si2 = s[i2] - v;
+
+ r[i] = sqrt(dot(si, si));
+ A[i] = cross(si, si2) / 2.0;
+ D[i] = dot(si, si2);
+ if (fabs(r[i]) < 1e-16) {
+ result[i] = 1.0;
+ return result;
+ } else if (fabs(A[i]) < 1e-16 && D[i] < 0.0) {
+ double r2 = sqrt(dot(si2, si2));
+ result[i2] = r[i] / (r[i] + r2);
+ result[i] = r2 / (r[i] + r2);
+ return result;
+ }
+ }
+
+ double w_sum = 0.0;
+
+ for (size_t i = 0; i < SZ; ++i) {
+ size_t i_m = (i + SZ - 1) % SZ;
+ size_t i_p = (i + 1) % SZ;
+
+ double w = 0.0;
+ if (fabs(A[i_m]) > 1e-16)
+ w += (r[i_m] - D[i_m] / r[i]) / A[i_m];
+ if (fabs(A[i]) > 1e-16)
+ w += (r[i_p] - D[i] / r[i]) / A[i];
+
+ result[i] = w;
+ w_sum += w;
+ }
+
+ for (size_t i = 0; i < SZ; ++i) {
+ result[i] /= w_sum;
+ }
+
+// carve::geom2d::P2 test;
+// for (size_t i = 0; i < SZ; ++i) {
+// test = test + result[i] * s[i];
+// }
+
+ return result;
+ }
+
+ template<typename iter_t,
+ typename adapt_t,
+ typename val_t,
+ typename mod_t>
+ val_t interp(iter_t begin,
+ iter_t end,
+ adapt_t adapt,
+ const std::vector<val_t> &vals,
+ double x,
+ double y,
+ mod_t mod = mod_t()) {
+ std::vector<carve::geom2d::P2> s;
+ s.reserve(std::distance(begin, end));
+ std::transform(begin, end, std::back_inserter(s), adapt);
+ std::vector<double> weight = polyInterpolate(s, carve::geom::VECTOR(x, y));
+
+ val_t v;
+ for (size_t z = 0; z < weight.size(); z++) {
+ v += weight[z] * vals[z];
+ }
+
+ return mod(v);
+ }
+
+ template<typename iter_t,
+ typename adapt_t,
+ typename val_t>
+ val_t interp(iter_t begin,
+ iter_t end,
+ adapt_t adapt,
+ const std::vector<val_t> &vals,
+ double x,
+ double y) {
+ return interp(begin, end, adapt, vals, x, y, identity_t<val_t>());
+ }
+
+ template<typename vertex_t,
+ typename adapt_t,
+ typename val_t,
+ typename mod_t>
+ val_t interp(const std::vector<vertex_t> &poly,
+ adapt_t adapt,
+ const std::vector<val_t> &vals,
+ double x,
+ double y,
+ mod_t mod = mod_t()) {
+ return interp(poly.begin(), poly.end(), adapt, vals, x, y, mod);
+ }
+
+ template<typename vertex_t,
+ typename adapt_t,
+ typename val_t>
+ val_t interp(const std::vector<vertex_t> &poly,
+ adapt_t adapt,
+ const std::vector<val_t> &vals,
+ double x,
+ double y) {
+ return interp(poly.begin(), poly.end(), adapt, vals, x, y, identity_t<val_t>());
+ }
+
+ template<typename val_t,
+ typename mod_t>
+ val_t interp(const std::vector<carve::geom2d::P2> &poly,
+ const std::vector<val_t> &vals,
+ double x,
+ double y,
+ mod_t mod = mod_t()) {
+ std::vector<double> weight = polyInterpolate(poly, carve::geom::VECTOR(x, y));
+
+ val_t v;
+ for (size_t z = 0; z < weight.size(); z++) {
+ v += weight[z] * vals[z];
+ }
+
+ return mod(v);
+ }
+
+ template<typename val_t>
+ val_t interp(const std::vector<carve::geom2d::P2> &poly,
+ const std::vector<val_t> &vals,
+ double x,
+ double y) {
+ return interp(poly, vals, x, y, identity_t<val_t>());
+ }
+
+ class Interpolator {
+ public:
+ virtual void interpolate(const carve::mesh::MeshSet<3>::face_t *new_face,
+ const carve::mesh::MeshSet<3>::face_t *orig_face,
+ bool flipped) =0;
+
+ Interpolator() {
+ }
+
+ virtual ~Interpolator() {
+ }
+
+ class Hook : public carve::csg::CSG::Hook {
+ Interpolator *interpolator;
+ public:
+ virtual void resultFace(const carve::mesh::MeshSet<3>::face_t *new_face,
+ const carve::mesh::MeshSet<3>::face_t *orig_face,
+ bool flipped) {
+ interpolator->interpolate(new_face, orig_face, flipped);
+ }
+
+ Hook(Interpolator *_interpolator) : interpolator(_interpolator) {
+ }
+
+ virtual ~Hook() {
+ }
+ };
+
+ void installHooks(carve::csg::CSG &csg) {
+ csg.hooks.registerHook(new Hook(this), carve::csg::CSG::Hooks::RESULT_FACE_BIT);
+ }
+ };
+
+ template<typename attr_t>
+ class FaceVertexAttr : public Interpolator {
+
+ protected:
+ struct fv_hash {
+ size_t operator()(const std::pair<const carve::mesh::MeshSet<3>::face_t *, unsigned> &v) const {
+ return size_t(v.first) ^ size_t(v.second);
+ }
+ };
+
+ typedef std::unordered_map<const carve::mesh::MeshSet<3>::vertex_t *, attr_t> attrvmap_t;
+ typedef std::unordered_map<std::pair<const carve::mesh::MeshSet<3>::face_t *, unsigned>, attr_t, fv_hash> attrmap_t;
+
+ attrmap_t attrs;
+
+ public:
+ bool hasAttribute(const carve::mesh::MeshSet<3>::face_t *f, unsigned v) {
+ return attrs.find(std::make_pair(f, v)) != attrs.end();
+ }
+
+ attr_t getAttribute(const carve::mesh::MeshSet<3>::face_t *f, unsigned v, const attr_t &def = attr_t()) {
+ typename attrmap_t::const_iterator fv = attrs.find(std::make_pair(f, v));
+ if (fv != attrs.end()) {
+ return (*fv).second;
+ }
+ return def;
+ }
+
+ void setAttribute(const carve::mesh::MeshSet<3>::face_t *f, unsigned v, const attr_t &attr) {
+ attrs[std::make_pair(f, v)] = attr;
+ }
+
+ virtual void interpolate(const carve::mesh::MeshSet<3>::face_t *new_face,
+ const carve::mesh::MeshSet<3>::face_t *orig_face,
+ bool flipped) {
+ std::vector<attr_t> vertex_attrs;
+ attrvmap_t base_attrs;
+ vertex_attrs.reserve(orig_face->nVertices());
+
+ for (carve::mesh::MeshSet<3>::face_t::const_edge_iter_t e = orig_face->begin(); e != orig_face->end(); ++e) {
+ typename attrmap_t::const_iterator a = attrs.find(std::make_pair(orig_face, e.idx()));
+ if (a == attrs.end()) return;
+ vertex_attrs.push_back((*a).second);
+ base_attrs[e->vert] = vertex_attrs.back();
+ }
+
+ for (carve::mesh::MeshSet<3>::face_t::const_edge_iter_t e = new_face->begin(); e != new_face->end(); ++e) {
+ const carve::mesh::MeshSet<3>::vertex_t *vertex = e->vert;
+ typename attrvmap_t::const_iterator b = base_attrs.find(vertex);
+ if (b != base_attrs.end()) {
+ attrs[std::make_pair(new_face, e.idx())] = (*b).second;
+ } else {
+ carve::geom2d::P2 p = orig_face->project(e->vert->v);
+ attr_t attr = interp(orig_face->begin(),
+ orig_face->end(),
+ orig_face->projector(),
+ vertex_attrs,
+ p.x,
+ p.y);
+ attrs[std::make_pair(new_face, e.idx())] = attr;
+ }
+ }
+ }
+
+ FaceVertexAttr() : Interpolator() {
+ }
+
+ virtual ~FaceVertexAttr() {
+ }
+
+ };
+
+
+ template<typename attr_t>
+ class FaceAttr : public Interpolator {
+
+ protected:
+ struct f_hash {
+ size_t operator()(const carve::mesh::MeshSet<3>::face_t * const &f) const {
+ return size_t(f);
+ }
+ };
+
+ typedef std::unordered_map<const carve::mesh::MeshSet<3>::face_t *, attr_t, f_hash> attrmap_t;
+
+ attrmap_t attrs;
+
+ public:
+ bool hasAttribute(const carve::mesh::MeshSet<3>::face_t *f) {
+ return attrs.find(f) != attrs.end();
+ }
+
+ attr_t getAttribute(const carve::mesh::MeshSet<3>::face_t *f, const attr_t &def = attr_t()) {
+ typename attrmap_t::const_iterator i = attrs.find(f);
+ if (i != attrs.end()) {
+ return (*i).second;
+ }
+ return def;
+ }
+
+ void setAttribute(const carve::mesh::MeshSet<3>::face_t *f, const attr_t &attr) {
+ attrs[f] = attr;
+ }
+
+ virtual void interpolate(const carve::mesh::MeshSet<3>::face_t *new_face,
+ const carve::mesh::MeshSet<3>::face_t *orig_face,
+ bool flipped) {
+ typename attrmap_t::const_iterator i = attrs.find(orig_face);
+ if (i != attrs.end()) {
+ attrs[new_face] = (*i).second;
+ }
+ }
+
+ FaceAttr() : Interpolator() {
+ }
+
+ virtual ~FaceAttr() {
+ }
+
+ };
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/collection_types.hpp>
+#include <carve/iobj.hpp>
+
+namespace carve {
+ namespace csg {
+
+ /**
+ * \class Intersections
+ * \brief Storage for computed intersections between vertices, edges and faces.
+ *
+ */
+ struct Intersections : public std::unordered_map<IObj, IObjVMapSmall, IObj_hash> {
+ typedef carve::mesh::MeshSet<3>::vertex_t vertex_t;
+ typedef carve::mesh::MeshSet<3>::edge_t edge_t;
+ typedef carve::mesh::MeshSet<3>::face_t face_t;
+
+ typedef std::unordered_map<IObj, IObjVMapSmall, IObj_hash> super;
+
+ ~Intersections() {
+ }
+
+ /**
+ * \brief Record the position of intersection between a pair of intersection objects.
+ *
+ * @param a The first intersecting object.
+ * @param b The second intersecting object.
+ * @param p The point of intersection.
+ */
+ void record(IObj a, IObj b, vertex_t *p) {
+ if (a > b) std::swap(a, b);
+ (*this)[a][b] = p;
+ (*this)[b][a] = p;
+ }
+
+ /**
+ * \brief Test whether vertex \a v intersects face \a f.
+ *
+ * @param v The vertex to test.
+ * @param f The face to test.
+ *
+ * @return true, if \a v intersects \a f.
+ */
+ bool intersectsFace(vertex_t *v, face_t *f) const;
+
+ /**
+ * \brief Collect sets of vertices, edges and faces that intersect \a obj
+ *
+ * @param[in] obj The intersection object to search for intersections.
+ * @param[out] collect_v A vector of vertices intersecting \a obj.
+ * @param[out] collect_e A vector of edges intersecting \a obj.
+ * @param[out] collect_f A vector of faces intersecting \a obj.
+ */
+ void collect(const IObj &obj,
+ std::vector<vertex_t *> *collect_v,
+ std::vector<edge_t *> *collect_e,
+ std::vector<face_t *> *collect_f) const;
+
+
+ /**
+ * \brief Determine whether two intersection objects intersect.
+ *
+ * @param a The first intersection object.
+ * @param b The second intersection object.
+ *
+ * @return true, if \a a and \a b intersect.
+ */
+ bool intersectsExactly(const IObj &a, const IObj &b) {
+ Intersections::const_iterator i = find(a);
+ if (i == end()) return false;
+ return i->second.find(b) != i->second.end();
+ }
+
+ /**
+ * \brief Determine whether an intersection object intersects a vertex.
+ *
+ * @param a The intersection object.
+ * @param v The vertex.
+ *
+ * @return true, if \a a and \a v intersect.
+ */
+ bool intersects(const IObj &a, vertex_t *v) {
+ Intersections::const_iterator i = find(a);
+ if (i == end()) return false;
+ if (i->second.find(v) != i->second.end()) return true;
+ return false;
+ }
+
+ /**
+ * \brief Determine whether an intersection object intersects an edge.
+ *
+ * @param a The intersection object.
+ * @param e The edge.
+ *
+ * @return true, if \a a and \a e intersect (either on the edge,
+ * or at either endpoint).
+ */
+ bool intersects(const IObj &a, edge_t *e) {
+ Intersections::const_iterator i = find(a);
+ if (i == end()) return false;
+ for (super::data_type::const_iterator j = i->second.begin(); j != i->second.end(); ++j) {
+ const IObj &obj = j->first;
+ switch (obj.obtype) {
+ case IObj::OBTYPE_VERTEX:
+ if (obj.vertex == e->v1() || obj.vertex == e->v2()) return true;
+ break;
+ case IObj::OBTYPE_EDGE:
+ if (obj.edge == e) return true;
+ break;
+ default:
+ break;
+ }
+ }
+ return false;
+ }
+
+ /**
+ * \brief Determine whether an intersection object intersects a face.
+ *
+ * @param a The intersection object.
+ * @param f The face.
+ *
+ * @return true, if \a a and \a f intersect (either on the face,
+ * or at any associated edge or vertex).
+ */
+ bool intersects(const IObj &a, face_t *f) {
+ Intersections::const_iterator i = find(a);
+ if (i == end()) return false;
+ if (i->second.find(f) != i->second.end()) return true;
+ edge_t *e = f->edge;
+ do {
+ if (i->second.find(e) != i->second.end()) return true;
+ if (i->second.find(e->vert) != i->second.end()) return true;
+ e = e->next;
+ } while (e != f->edge);
+ return false;
+ }
+
+ /**
+ * \brief Determine whether an edge intersects another edge.
+ *
+ * @param e The edge.
+ * @param f The face.
+ *
+ * @return true, if \a e and \a f intersect.
+ */
+ bool intersects(edge_t *e1, edge_t *e2) {
+ if (intersects(e1->v1(), e2) || intersects(e1->v2(), e2) || intersects(IObj(e1), e2)) return true;
+ return false;
+ }
+
+ /**
+ * \brief Determine whether an edge intersects a face.
+ *
+ * @param e The edge.
+ * @param f The face.
+ *
+ * @return true, if \a e and \a f intersect.
+ */
+ bool intersects(edge_t *e, face_t *f) {
+ if (intersects(e->v1(), f) || intersects(e->v2(), f) || intersects(IObj(e), f)) return true;
+ return false;
+ }
+
+ /**
+ * \brief Determine the faces intersected by an edge.
+ *
+ * @tparam face_set_t A collection type holding face_t *
+ * @param[in] e The edge.
+ * @param[out] f The resulting set of faces.
+ */
+ template<typename face_set_t>
+ void intersectedFaces(edge_t *e, face_set_t &f) const {
+ std::vector<face_t *> intersected_faces;
+ std::vector<edge_t *> intersected_edges;
+ std::vector<vertex_t *> intersected_vertices;
+
+ collect(e, &intersected_vertices, &intersected_edges, &intersected_faces);
+
+ for (unsigned i = 0; i < intersected_vertices.size(); ++i) {
+ facesForVertex(intersected_vertices[i], f);
+ }
+ for (unsigned i = 0; i < intersected_edges.size(); ++i) {
+ facesForEdge(intersected_edges[i], f);
+ }
+ f.insert(intersected_faces.begin(), intersected_faces.end());
+ }
+
+ /**
+ * \brief Determine the faces intersected by a vertex.
+ *
+ * @tparam face_set_t A collection type holding face_t *
+ * @param[in] v The vertex.
+ * @param[out] f The resulting set of faces.
+ */
+ template<typename face_set_t>
+ void intersectedFaces(vertex_t *v, face_set_t &f) const {
+ std::vector<face_t *> intersected_faces;
+ std::vector<edge_t *> intersected_edges;
+ std::vector<vertex_t *> intersected_vertices;
+
+ collect(v, &intersected_vertices, &intersected_edges, &intersected_faces);
+
+ for (unsigned i = 0; i < intersected_vertices.size(); ++i) {
+ facesForVertex(intersected_vertices[i], f);
+ }
+ for (unsigned i = 0; i < intersected_edges.size(); ++i) {
+ facesForEdge(intersected_edges[i], f);
+ }
+ f.insert(intersected_faces.begin(), intersected_faces.end());
+ }
+
+ /**
+ * \brief Collect the set of faces that contain all vertices in \a verts.
+ *
+ * @tparam vertex_set_t A collection type holding vertex_t *
+ * @tparam face_set_t A collection type holding face_t *
+ * @param[in] verts A set of vertices.
+ * @param[out] result The resulting set of faces.
+ */
+ template<typename vertex_set_t, typename face_set_t>
+ void commonFaces(const vertex_set_t &verts, face_set_t &result) {
+
+ std::set<face_t *> ifaces, temp, out;
+ typename vertex_set_t::const_iterator i = verts.begin();
+ if (i == verts.end()) return;
+ intersectedFaces((*i), ifaces);
+ while (++i != verts.end()) {
+ temp.clear();
+ intersectedFaces((*i), temp);
+
+ out.clear();
+ std::set_intersection(temp.begin(), temp.end(),
+ ifaces.begin(), ifaces.end(),
+ set_inserter(out));
+ ifaces.swap(out);
+ }
+ std::copy(ifaces.begin(), ifaces.end(), set_inserter(result));
+ }
+
+ void clear() {
+ super::clear();
+ }
+
+ };
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+namespace carve {
+ namespace csg {
+ struct IObj {
+ enum {
+ OBTYPE_NONE = 0,
+ OBTYPE_VERTEX = 1,
+ OBTYPE_EDGE = 2,
+ OBTYPE_FACE = 4
+ } obtype;
+
+ union {
+ carve::mesh::MeshSet<3>::vertex_t *vertex;
+ carve::mesh::MeshSet<3>::edge_t *edge;
+ carve::mesh::MeshSet<3>::face_t *face;
+ intptr_t val;
+ };
+
+ IObj() : obtype(OBTYPE_NONE), val(0) { }
+ IObj(carve::mesh::MeshSet<3>::vertex_t *v) : obtype(OBTYPE_VERTEX), vertex(v) { }
+ IObj(carve::mesh::MeshSet<3>::edge_t *e) : obtype(OBTYPE_EDGE), edge(e) { }
+ IObj(carve::mesh::MeshSet<3>::face_t *f) : obtype(OBTYPE_FACE), face(f) { }
+ char typeChar() const { return "NVExF"[obtype]; }
+ };
+
+
+
+ struct IObj_hash {
+ inline size_t operator()(const IObj &i) const {
+ return (size_t)i.val;
+ }
+ inline size_t operator()(const std::pair<const IObj, const IObj> &i) const {
+ return (size_t)i.first.val ^ (size_t)i.second.val;
+ }
+ };
+
+
+
+ typedef std::unordered_set<std::pair<const IObj, const IObj>, IObj_hash> IObjPairSet;
+
+ typedef std::unordered_map<IObj, carve::mesh::MeshSet<3>::vertex_t *, IObj_hash> IObjVMap;
+ typedef std::map<IObj, carve::mesh::MeshSet<3>::vertex_t *> IObjVMapSmall;
+
+ class VertexIntersections :
+ public std::unordered_map<carve::mesh::MeshSet<3>::vertex_t *, IObjPairSet> {
+ };
+
+
+
+ static inline bool operator==(const carve::csg::IObj &a, const carve::csg::IObj &b) {
+ return a.obtype == b.obtype && a.val == b.val;
+ }
+
+ static inline bool operator!=(const carve::csg::IObj &a, const carve::csg::IObj &b) {
+ return a.obtype != b.obtype || a.val != b.val;
+ }
+
+ static inline bool operator<(const carve::csg::IObj &a, const carve::csg::IObj &b) {
+ return a.obtype < b.obtype || (a.obtype == b.obtype && a.val < b.val);
+ }
+
+ static inline bool operator<=(const carve::csg::IObj &a, const carve::csg::IObj &b) {
+ return a.obtype < b.obtype || (a.obtype == b.obtype && a.val <= b.val);
+ }
+
+ static inline bool operator>(const carve::csg::IObj &a, const carve::csg::IObj &b) {
+ return a.obtype > b.obtype || (a.obtype == b.obtype && a.val > b.val);
+ }
+
+ static inline bool operator>=(const carve::csg::IObj &a, const carve::csg::IObj &b) {
+ return a.obtype > b.obtype || (a.obtype == b.obtype && a.val >= b.val);
+ }
+
+ static inline std::ostream &operator<<(std::ostream &o, const carve::csg::IObj &a) {
+ switch (a.obtype) {
+ case carve::csg::IObj::OBTYPE_NONE: o << "NONE{}"; break;
+ case carve::csg::IObj::OBTYPE_VERTEX: o << "VERT{" << a.vertex << "}"; break;
+ case carve::csg::IObj::OBTYPE_EDGE: o << "EDGE{" << a.edge << "}"; break;
+ case carve::csg::IObj::OBTYPE_FACE: o << "FACE{" << a.face << "}"; break;
+ }
+ return o;
+ }
+
+ }
+}
+
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+#include <carve/geom.hpp>
+#include <carve/aabb.hpp>
+
+#include <queue>
+#include <list>
+#include <limits>
+
+namespace carve {
+ namespace geom {
+ template<unsigned ndim,
+ typename data_t,
+ typename inserter_t,
+ typename aabb_calc_t>
+ class kd_node {
+ kd_node(const kd_node &);
+ kd_node &operator=(const kd_node &);
+
+ public:
+ kd_node *c_neg;
+ kd_node *c_pos;
+ kd_node *parent;
+ axis_pos splitpos;
+
+ typedef vector<ndim> vector_t;
+ typedef std::list<data_t> container_t;
+
+ container_t data;
+
+ kd_node(kd_node *_parent = NULL) : c_neg(NULL), c_pos(NULL), parent(_parent), splitpos(0, 0.0) {
+ }
+
+ ~kd_node() {
+ if (c_neg) delete c_neg;
+ if (c_pos) delete c_pos;
+ }
+
+ template<typename visitor_t>
+ void closeNodes(const vector_t &p, double d, visitor_t &visit) const {
+ if (c_neg) {
+ double delta = splitpos.pos - p[splitpos.axis];
+ if (delta <= d) c_neg->closeNodes(p, d, visit);
+ if (delta >= -d) c_pos->closeNodes(p, d, visit);
+ } else {
+ visit(this);
+ }
+ }
+
+ void removeData(const data_t &d) {
+ typename container_t::iterator i = std::find(data.begin(), data.end(), d);
+
+ if (i != data.end()) {
+ data.erase(i);
+ }
+ }
+
+ void addData(const data_t &d) {
+ data.push_back(d);
+ }
+
+ void insert(const data_t &data, inserter_t &inserter) {
+ inserter.insert(this, data);
+ }
+
+ void insert(const data_t &data) {
+ inserter_t inserter;
+ insert(data, inserter);
+ }
+
+ void remove(const data_t &data, inserter_t &inserter) {
+ inserter.remove(this, data);
+ }
+
+ void remove(const data_t &data) {
+ inserter_t inserter;
+ remove(data, inserter);
+ }
+
+ carve::geom::aabb<ndim> nodeAABB() const {
+ carve::geom::aabb<ndim> aabb;
+ if (c_neg) {
+ aabb = c_neg->nodeAABB();
+ aabb.unionAABB(c_pos->nodeAABB());
+ } else {
+ if (data.size()) {
+ typename container_t::const_iterator i = data.begin();
+ aabb = aabb_calc_t()(*i);
+ while (i != data.end()) {
+ aabb.unionAABB(aabb_calc_t()(*i));
+ ++i;
+ }
+ }
+ }
+ return aabb;
+ }
+
+ bool split(axis_pos split_at, inserter_t &inserter) {
+ if (c_neg) {
+ // already split
+ return false;
+ }
+
+ c_neg = new kd_node(this);
+ c_pos = new kd_node(this);
+
+ // choose an axis and split point.
+ splitpos = split_at;
+
+ carve::geom::aabb<ndim> aabb;
+
+ if (splitpos.axis < 0 ||
+ splitpos.axis >= ndim ||
+ splitpos.pos == std::numeric_limits<double>::max()) {
+ // need an aabb
+ if (data.size()) {
+ typename container_t::const_iterator i = data.begin();
+ aabb = aabb_calc_t()(*i);
+ while (i != data.end()) {
+ aabb.unionAABB(aabb_calc_t()(*i));
+ ++i;
+ }
+ }
+ }
+
+ if (splitpos.axis < 0 || splitpos.axis >= ndim) {
+
+ // choose an axis;
+
+ // if no axis was specified, force calculation of the split position.
+ splitpos.pos = std::numeric_limits<double>::max();
+
+ // choose the axis of the AABB with the biggest extent.
+ splitpos.axis = largestAxis(aabb.extent);
+
+ if (parent && splitpos.axis == parent->splitpos.axis) {
+ // but don't choose the same axis as the parent node;
+ // choose the axis with the second greatest AABB extent.
+ double e = -1.0;
+ int a = -1;
+ for (unsigned i = 0; i < ndim; ++i) {
+ if (i == splitpos.axis) continue;
+ if (e < aabb.extent[i]) { a = i; e = aabb.extent[i]; }
+ }
+ if (a != -1) {
+ splitpos.axis = a;
+ }
+ }
+ }
+
+ if (splitpos.pos == std::numeric_limits<double>::max()) {
+ carve::geom::vector<ndim> min = aabb.min();
+ carve::geom::vector<ndim> max = aabb.max();
+ splitpos.pos = aabb.pos.v[splitpos.axis];
+ }
+
+ inserter.propagate(this);
+
+ return true;
+ }
+
+ bool split(axis_pos split_at = axis_pos(-1, std::numeric_limits<double>::max())) {
+ inserter_t inserter;
+ return split(split_at, inserter);
+ }
+
+ void splitn(int num, inserter_t &inserter) {
+ if (num <= 0) return;
+ if (!c_neg) {
+ split(inserter);
+ }
+ if (c_pos) c_pos->splitn(num-1, inserter);
+ if (c_neg) c_neg->splitn(num-1, inserter);
+ }
+
+ void splitn(int num) {
+ inserter_t inserter;
+ splitn(num, inserter);
+ }
+
+ template<typename split_t>
+ void splitn(int num, split_t splitter, inserter_t &inserter) {
+ if (num <= 0) return;
+ if (!c_neg) {
+ split(inserter, splitter(this));
+ }
+ if (c_pos) c_pos->splitn(num-1, inserter, splitter);
+ if (c_neg) c_neg->splitn(num-1, inserter, splitter);
+ }
+
+ template<typename split_t>
+ void splitn(int num, split_t splitter) {
+ inserter_t inserter;
+ splitn(num, splitter, inserter);
+ }
+
+ template<typename pred_t>
+ void splitpred(pred_t pred, inserter_t &inserter, int depth = 0) {
+ if (!c_neg) {
+ axis_pos splitpos(-1, std::numeric_limits<double>::max());
+ if (!pred(this, depth, splitpos)) return;
+ split(splitpos, inserter);
+ }
+ if (c_pos) c_pos->splitpred(pred, inserter, depth + 1);
+ if (c_neg) c_neg->splitpred(pred, inserter, depth + 1);
+ }
+
+ template<typename pred_t>
+ void splitpred(pred_t pred, int depth = 0) {
+ inserter_t inserter;
+ splitpred(pred, inserter, depth);
+ }
+
+ // distance_t must provide:
+ // double operator()(kd_node::data_t, vector<ndim>);
+ // double operator()(axis_pos, vector<ndim>);
+ template<typename distance_t>
+ struct near_point_query {
+
+ // q_t - the priority queue value type.
+ // q_t.first: distance from object to query point.
+ // q_t.second: pointer to object
+ typedef std::pair<double, const typename kd_node::data_t *> q_t;
+
+ // the queue priority should sort from smallest distance to largest, and on equal distance, by object pointer.
+ struct pcmp {
+ bool operator()(const q_t &a, const q_t &b) {
+ return (a.first > b.first) || ((a.first == b.first) && (a.second < b.second));
+ }
+ };
+
+ vector<ndim> point;
+ const kd_node *node;
+ std::priority_queue<q_t, std::vector<q_t>, pcmp> pq;
+
+ distance_t dist;
+ double dist_to_parent_split;
+
+ void addToPQ(kd_node *node) {
+ if (node->c_neg) {
+ addToPQ(node->c_neg);
+ addToPQ(node->c_pos);
+ } else {
+ for (typename kd_node::container_t::const_iterator i = node->data.begin(); i != node->data.end(); ++i) {
+ double d = dist((*i), point);
+ pq.push(std::make_pair(d, &(*i)));
+ }
+ }
+ }
+
+ const typename kd_node::data_t *next() {
+ while (1) {
+ if (pq.size()) {
+ q_t t = pq.top();
+ if (!node->parent || t.first < dist_to_parent_split) {
+ pq.pop();
+ return t.second;
+ }
+ }
+
+ if (!node->parent) return NULL;
+
+ if (node->parent->c_neg == node) {
+ addToPQ(node->parent->c_pos);
+ } else {
+ addToPQ(node->parent->c_neg);
+ }
+
+ node = node->parent;
+ dist_to_parent_split = dist(node->splitpos, point);
+ }
+ }
+
+ near_point_query(const vector<ndim> _point, const kd_node *_node) : point(_point), node(_node), pq(), dist() {
+ while (node->c_neg) {
+ node = (point[node->axis] < node->pos) ? node->c_neg : node->c_pos;
+ }
+ if (node->parent) {
+ dist_to_parent_split = dist(node->parent->splitpos, point);
+ } else {
+ dist_to_parent_split = HUGE_VAL;
+ }
+ addToPQ(node);
+ }
+ };
+
+ };
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/math_constants.hpp>
+
+#include <math.h>
+
+namespace carve {
+ namespace geom {
+ template<unsigned ndim> struct vector;
+ }
+}
+
+namespace carve {
+ namespace math {
+ struct Matrix3;
+ int cubic_roots(double c3, double c2, double c1, double c0, double *roots);
+
+ void eigSolveSymmetric(const Matrix3 &m,
+ double &l1, carve::geom::vector<3> &e1,
+ double &l2, carve::geom::vector<3> &e2,
+ double &l3, carve::geom::vector<3> &e3);
+
+ void eigSolve(const Matrix3 &m, double &l1, double &l2, double &l3);
+
+ static inline bool ZERO(double x) { return fabs(x) < carve::EPSILON; }
+
+ static inline double radians(double deg) { return deg * M_PI / 180.0; }
+ static inline double degrees(double rad) { return rad * 180.0 / M_PI; }
+
+ static inline double ANG(double x) {
+ return (x < 0) ? x + M_TWOPI : x;
+ }
+
+ template<typename T>
+ static inline const T &clamp(const T &val, const T &min, const T &max) {
+ if (val < min) return min;
+ if (val > max) return max;
+ return val;
+ }
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <math.h>
+
+#ifndef M_SQRT_3
+#define M_SQRT_3 1.73205080756887729352
+#endif
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+#ifndef M_TWOPI
+#define M_TWOPI (M_PI + M_PI)
+#endif
+
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <cstring>
+
+#include <carve/carve.hpp>
+
+#include <carve/math.hpp>
+#include <carve/geom.hpp>
+
+namespace carve {
+ namespace math {
+
+ struct Quaternion {
+ double x, y, z, w;
+
+ Quaternion(double _x, double _y, double _z, double _w) : x(_x), y(_y), z(_z), w(_w) {
+ }
+
+ Quaternion(double angle, const carve::geom::vector<3> &axis) {
+ double s = axis.length();
+ if (!carve::math::ZERO(s)) {
+ double c = 1.0 / s;
+ double omega = -0.5 * angle;
+ s = sin(omega);
+ x = axis.x * c * s;
+ y = axis.y * c * s;
+ z = axis.z * c * s;
+ w = cos(omega);
+ normalize();
+ } else {
+ x = y = z = 0.0;
+ w = 1.0;
+ }
+ }
+
+ double lengthSquared() const {
+ return x * x + y * y + z * z + w * w;
+ }
+
+ double length() const {
+ return sqrt(lengthSquared());
+ }
+
+ Quaternion normalized() const {
+ return Quaternion(*this).normalize();
+ }
+
+ Quaternion &normalize() {
+ double l = length();
+ if (l == 0.0) {
+ x = 1.0; y = 0.0; z = 0.0; w = 0.0;
+ } else {
+ x /= l; y /= l; z /= l; w /= l;
+ }
+ return *this;
+ }
+ };
+
+ struct Matrix3 {
+ // access: .m[col][row], .v[col * 4 + row], ._cr
+ union {
+ double m[3][3];
+ double v[9];
+ struct {
+ // transposed
+ double _11, _12, _13;
+ double _21, _22, _23;
+ double _31, _32, _33;
+ };
+ };
+ Matrix3(double __11, double __21, double __31,
+ double __12, double __22, double __32,
+ double __13, double __23, double __33) {
+ // nb, args are row major, storage is column major.
+ _11 = __11; _12 = __12; _13 = __13;
+ _21 = __21; _22 = __22; _23 = __23;
+ _31 = __31; _32 = __32; _33 = __33;
+ }
+ Matrix3(double _m[3][3]) {
+ std::memcpy(m, _m, sizeof(m));
+ }
+ Matrix3(double _v[9]) {
+ std::memcpy(v, _v, sizeof(v));
+ }
+ Matrix3() {
+ _11 = 1.00; _12 = 0.00; _13 = 0.00;
+ _21 = 0.00; _22 = 1.00; _23 = 0.00;
+ _31 = 0.00; _32 = 0.00; _33 = 1.00;
+ }
+ };
+
+ struct Matrix {
+ // access: .m[col][row], .v[col * 4 + row], ._cr
+ union {
+ double m[4][4];
+ double v[16];
+ struct {
+ // transposed
+ double _11, _12, _13, _14;
+ double _21, _22, _23, _24;
+ double _31, _32, _33, _34;
+ double _41, _42 ,_43, _44;
+ };
+ };
+ Matrix(double __11, double __21, double __31, double __41,
+ double __12, double __22, double __32, double __42,
+ double __13, double __23, double __33, double __43,
+ double __14, double __24, double __34, double __44) {
+ // nb, args are row major, storage is column major.
+ _11 = __11; _12 = __12; _13 = __13; _14 = __14;
+ _21 = __21; _22 = __22; _23 = __23; _24 = __24;
+ _31 = __31; _32 = __32; _33 = __33; _34 = __34;
+ _41 = __41; _42 = __42; _43 = __43; _44 = __44;
+ }
+ Matrix(double _m[4][4]) {
+ std::memcpy(m, _m, sizeof(m));
+ }
+ Matrix(double _v[16]) {
+ std::memcpy(v, _v, sizeof(v));
+ }
+ Matrix() {
+ _11 = 1.00; _12 = 0.00; _13 = 0.00; _14 = 0.00;
+ _21 = 0.00; _22 = 1.00; _23 = 0.00; _24 = 0.00;
+ _31 = 0.00; _32 = 0.00; _33 = 1.00; _34 = 0.00;
+ _41 = 0.00; _42 = 0.00; _43 = 0.00; _44 = 1.00;
+ }
+
+ static Matrix ROT(const Quaternion &q) {
+ const double w = q.w;
+ const double x = q.x;
+ const double y = q.y;
+ const double z = q.z;
+ return Matrix(1 - 2*y*y - 2*z*z, 2*x*y - 2*z*w, 2*x*z + 2*y*w, 0.0,
+ 2*x*y + 2*z*w, 1 - 2*x*x - 2*z*z, 2*y*z - 2*x*w, 0.0,
+ 2*x*z - 2*y*w, 2*y*z + 2*x*w, 1 - 2*x*x - 2*y*y, 0.0,
+ 0.0, 0.0, 0.0, 1.0);
+ }
+ static Matrix ROT(double angle, const carve::geom::vector<3> &axis) {
+ return ROT(Quaternion(angle, axis));
+ }
+ static Matrix ROT(double angle, double x, double y, double z) {
+ return ROT(Quaternion(angle, carve::geom::VECTOR(x, y, z)));
+ }
+ static Matrix TRANS(double x, double y, double z) {
+ return Matrix(1.0, 0.0, 0.0, x,
+ 0.0, 1.0, 0.0, y,
+ 0.0, 0.0, 1.0, z,
+ 0.0, 0.0, 0.0, 1.0);
+ }
+ static Matrix TRANS(const carve::geom::vector<3> &v) {
+ return TRANS(v.x, v.y, v.z);
+ }
+ static Matrix SCALE(double x, double y, double z) {
+ return Matrix(x, 0.0, 0.0, 0.0,
+ 0.0, y, 0.0, 0.0,
+ 0.0, 0.0, z, 0.0,
+ 0.0, 0.0, 0.0, 1.0);
+ }
+ static Matrix SCALE(const carve::geom::vector<3> &v) {
+ return SCALE(v.x, v.y, v.z);
+ }
+ static Matrix IDENT() {
+ return Matrix(1.0, 0.0, 0.0, 0.0,
+ 0.0, 1.0, 0.0, 0.0,
+ 0.0, 0.0, 1.0, 0.0,
+ 0.0, 0.0, 0.0, 1.0);
+ }
+ };
+
+ static inline bool operator==(const Matrix &A, const Matrix &B) {
+ for (size_t i = 0; i < 16; ++i) if (A.v[i] != B.v[i]) return false;
+ return true;
+ }
+ static inline bool operator!=(const Matrix &A, const Matrix &B) {
+ return !(A == B);
+ }
+ static inline carve::geom::vector<3> operator*(const Matrix &A, const carve::geom::vector<3> &b) {
+ return carve::geom::VECTOR(
+ A._11 * b.x + A._21 * b.y + A._31 * b.z + A._41,
+ A._12 * b.x + A._22 * b.y + A._32 * b.z + A._42,
+ A._13 * b.x + A._23 * b.y + A._33 * b.z + A._43
+ );
+ }
+
+ static inline carve::geom::vector<3> &operator*=(carve::geom::vector<3> &b, const Matrix &A) {
+ b = A * b;
+ return b;
+ }
+
+ static inline carve::geom::vector<3> operator*(const Matrix3 &A, const carve::geom::vector<3> &b) {
+ return carve::geom::VECTOR(
+ A._11 * b.x + A._21 * b.y + A._31 * b.z,
+ A._12 * b.x + A._22 * b.y + A._32 * b.z,
+ A._13 * b.x + A._23 * b.y + A._33 * b.z
+ );
+ }
+
+ static inline carve::geom::vector<3> &operator*=(carve::geom::vector<3> &b, const Matrix3 &A) {
+ b = A * b;
+ return b;
+ }
+
+ static inline Matrix operator*(const Matrix &A, const Matrix &B) {
+ Matrix c;
+ for (int i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ c.m[i][j] = 0.0;
+ for (int k = 0; k < 4; k++) {
+ c.m[i][j] += A.m[k][j] * B.m[i][k];
+ }
+ }
+ }
+ return c;
+ }
+
+ static inline Matrix3 operator*(const Matrix3 &A, const Matrix3 &B) {
+ Matrix3 c;
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ c.m[i][j] = 0.0;
+ for (int k = 0; k < 3; k++) {
+ c.m[i][j] += A.m[k][j] * B.m[i][k];
+ }
+ }
+ }
+ return c;
+ }
+
+
+
+ struct matrix_transformation {
+ Matrix matrix;
+
+ matrix_transformation(const Matrix &_matrix) : matrix(_matrix) {
+ }
+
+ carve::geom::vector<3> operator()(const carve::geom::vector<3> &vector) const {
+ return matrix * vector;
+ }
+ };
+
+
+
+ }
+}
--- /dev/null
+// Begin License:
+// Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
+// All rights reserved.
+//
+// This file is part of the Carve CSG Library (http://carve-csg.com/)
+//
+// This file may be used under the terms of the GNU General Public
+// License version 2.0 as published by the Free Software Foundation
+// and appearing in the file LICENSE.GPL2 included in the packaging of
+// this file.
+//
+// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
+// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE.
+// End:
+
+
+#pragma once
+
+#include <carve/carve.hpp>
+
+#include <carve/geom.hpp>
+#include <carve/geom3d.hpp>
+#include <carve/tag.hpp>
+#include <carve/djset.hpp>
+#include <carve/aabb.hpp>
+#include <carve/rtree.hpp>
+
+#include <iostream>
+
+namespace carve {
+ namespace poly {
+ class Polyhedron;
+ }
+
+ namespace mesh {
+
+
+ template<unsigned ndim> class Edge;
+ template<unsigned ndim> class Face;
+ template<unsigned ndim> class Mesh;
+ template<unsigned ndim> class MeshSet;
+
+
+
+ // A Vertex may participate in several meshes. If the Mesh belongs
+ // to a MeshSet, then the vertices come from the vertex_storage
+ // member of the MeshSet. This allows one to construct one or more
+ // Meshes out of sets of connected faces (possibly using vertices
+ // from a variety of MeshSets, and other storage), then create a
+ // MeshSet from the Mesh(es), causing the vertices to be
+ // collected, cloned and repointed into the MeshSet.
+
+ // Normally, in a half-edge structure, Vertex would have a member
+ // pointing to an incident edge, allowing the enumeration of
+ // adjacent faces and edges. Because we want to support vertex
+ // sharing between Meshes and groups of Faces, this is made more
+ // complex. If Vertex contained a list of incident edges, one from
+ // each disjoint face set, then this could be done (with the
+ // caveat that you'd need to pass in a Mesh pointer to the
+ // adjacency queries). However, it seems that this would
+ // unavoidably complicate the process of incorporating or removing
+ // a vertex into an edge.
+
+ // In most cases it is expected that a vertex will be arrived at
+ // via an edge or face in the mesh (implicit or explicit) of
+ // interest, so not storing this information will not hurt,
+ // overly.
+ template<unsigned ndim>
+ class Vertex : public tagable {
+ public:
+ typedef carve::geom::vector<ndim> vector_t;
+ typedef MeshSet<ndim> owner_t;
+ typedef carve::geom::aabb<ndim> aabb_t;
+
+ carve::geom::vector<ndim> v;
+
+ Vertex(const vector_t &_v) : tagable(), v(_v) {
+ }
+
+ Vertex() : tagable(), v() {
+ }
+
+ aabb_t getAABB() const {
+ return aabb_t(v, carve::geom::vector<ndim>::ZERO());
+ }
+ };
+
+
+
+ struct hash_vertex_pair {
+ template<unsigned ndim>
+ size_t operator()(const std::pair<Vertex<ndim> *, Vertex<ndim> *> &pair) const {
+ size_t r = (size_t)pair.first;
+ size_t s = (size_t)pair.second;
+ return r ^ ((s >> 16) | (s << 16));
+ &nb
projects / blender-staging.git / commitdiff