Revision: 77311
http://sourceforge.net/p/brlcad/code/77311
Author: starseeker
Date: 2020-10-01 20:42:02 +0000 (Thu, 01 Oct 2020)
Log Message:
-----------
Upstream has moved quite a bit from this version of SPSR, we've added the
library call to it which still isn't in upstream, and experiments with the
latest PoissonRecon project code indicated it now takes a huge amount of RAM to
compile out of the box. Strip this down somewhat and move it into a subdir of
libbg, since a system version is currently guaranteed not to be suitable... No
particular point in paying the third party build maintenance overhead right now.
Modified Paths:
--------------
brlcad/trunk/INSTALL
brlcad/trunk/configure
brlcad/trunk/doc/legal/embedded/SPSR.txt
brlcad/trunk/misc/CMake/BRLCAD_Targets.cmake
brlcad/trunk/src/external/Creo/CMakeLists.txt
brlcad/trunk/src/libbg/CMakeLists.txt
brlcad/trunk/src/libbg/spsr.c
brlcad/trunk/src/other/CMakeLists.txt
Added Paths:
-----------
brlcad/trunk/src/libbg/spsr/
brlcad/trunk/src/libbg/spsr/Allocator.h
brlcad/trunk/src/libbg/spsr/Array.h
brlcad/trunk/src/libbg/spsr/Array.inl
brlcad/trunk/src/libbg/spsr/BSplineData.h
brlcad/trunk/src/libbg/spsr/BSplineData.inl
brlcad/trunk/src/libbg/spsr/BinaryNode.h
brlcad/trunk/src/libbg/spsr/Factor.cpp
brlcad/trunk/src/libbg/spsr/Factor.h
brlcad/trunk/src/libbg/spsr/FunctionData.h
brlcad/trunk/src/libbg/spsr/FunctionData.inl
brlcad/trunk/src/libbg/spsr/Geometry.cpp
brlcad/trunk/src/libbg/spsr/Geometry.h
brlcad/trunk/src/libbg/spsr/Geometry.inl
brlcad/trunk/src/libbg/spsr/Hash.h
brlcad/trunk/src/libbg/spsr/LICENSE
brlcad/trunk/src/libbg/spsr/MAT.h
brlcad/trunk/src/libbg/spsr/MAT.inl
brlcad/trunk/src/libbg/spsr/MarchingCubes.cpp
brlcad/trunk/src/libbg/spsr/MarchingCubes.h
brlcad/trunk/src/libbg/spsr/MemoryUsage.h
brlcad/trunk/src/libbg/spsr/MultiGridOctreeData.IsoSurface.inl
brlcad/trunk/src/libbg/spsr/MultiGridOctreeData.SortedTreeNodes.inl
brlcad/trunk/src/libbg/spsr/MultiGridOctreeData.h
brlcad/trunk/src/libbg/spsr/MultiGridOctreeData.inl
brlcad/trunk/src/libbg/spsr/Octree.h
brlcad/trunk/src/libbg/spsr/Octree.inl
brlcad/trunk/src/libbg/spsr/PPolynomial.h
brlcad/trunk/src/libbg/spsr/PPolynomial.inl
brlcad/trunk/src/libbg/spsr/PlyVertexMini.h
brlcad/trunk/src/libbg/spsr/PointStream.h
brlcad/trunk/src/libbg/spsr/PointStream.inl
brlcad/trunk/src/libbg/spsr/Polynomial.h
brlcad/trunk/src/libbg/spsr/Polynomial.inl
brlcad/trunk/src/libbg/spsr/SPSR.cpp
brlcad/trunk/src/libbg/spsr/SPSR.h
brlcad/trunk/src/libbg/spsr/SparseMatrix.h
brlcad/trunk/src/libbg/spsr/SparseMatrix.inl
brlcad/trunk/src/libbg/spsr/Vector.h
brlcad/trunk/src/libbg/spsr/Vector.inl
brlcad/trunk/src/libbg/spsr/cvertex.h
Removed Paths:
-------------
brlcad/trunk/src/other/libspsr/
brlcad/trunk/src/other/libspsr.dist
Modified: brlcad/trunk/INSTALL
===================================================================
--- brlcad/trunk/INSTALL 2020-10-01 19:56:41 UTC (rev 77310)
+++ brlcad/trunk/INSTALL 2020-10-01 20:42:02 UTC (rev 77311)
@@ -661,16 +661,6 @@
Aliases: ENABLE_OPENNURBS
---- BRLCAD_SPSR ---
-
-Option for enabling and disabling compilation of the Screened Poisson
-Surface Reconstruction library provided with BRL-CAD's source code.
-Default is AUTO, responsive to the toplevel BRLCAD_BUNDLED_LIBS option
-and testing first for a system version if BRLCAD_BUNDLED_LIBS is also AUTO.
-
-Aliases: ENABLE_SPSR
-
-
--- BRLCAD_SC ---
Option for enabling and disabling compilation of the NIST Step Class
Modified: brlcad/trunk/configure
===================================================================
--- brlcad/trunk/configure 2020-10-01 19:56:41 UTC (rev 77310)
+++ brlcad/trunk/configure 2020-10-01 20:42:02 UTC (rev 77311)
@@ -162,10 +162,6 @@
shift;;
--disable-opennurbs) options="$options
-DBRLCAD_OPENNURBS=SYSTEM";
shift;;
- --enable-spsr) options="$options -DBRLCAD_SPSR=BUNDLED";
- shift;;
- --disable-spsr) options="$options -DBRLCAD_SPSR=SYSTEM";
- shift;;
--enable-scl) options="$options -DBRLCAD_SC=BUNDLED";
shift;;
--disable-scl) options="$options -DBRLCAD_SC=SYSTEM";
Modified: brlcad/trunk/doc/legal/embedded/SPSR.txt
===================================================================
--- brlcad/trunk/doc/legal/embedded/SPSR.txt 2020-10-01 19:56:41 UTC (rev
77310)
+++ brlcad/trunk/doc/legal/embedded/SPSR.txt 2020-10-01 20:42:02 UTC (rev
77311)
@@ -85,49 +85,43 @@
reason for such modification.
file:/include/bg/spsr.h
-file:/src/other/libspsr/Src/Allocator.h
-file:/src/other/libspsr/Src/Array.h
-file:/src/other/libspsr/Src/Array.inl
-file:/src/other/libspsr/Src/BSplineData.h
-file:/src/other/libspsr/Src/BSplineData.inl
-file:/src/other/libspsr/Src/BinaryNode.h
-file:/src/other/libspsr/Src/CmdLineParser.cpp
-file:/src/other/libspsr/Src/CmdLineParser.h
-file:/src/other/libspsr/Src/CmdLineParser.inl
-file:/src/other/libspsr/Src/Factor.cpp
-file:/src/other/libspsr/Src/Factor.h
-file:/src/other/libspsr/Src/FunctionData.h
-file:/src/other/libspsr/Src/FunctionData.inl
-file:/src/other/libspsr/Src/Geometry.cpp
-file:/src/other/libspsr/Src/Geometry.h
-file:/src/other/libspsr/Src/Geometry.inl
-file:/src/other/libspsr/Src/Hash.h
-file:/src/other/libspsr/Src/MAT.h
-file:/src/other/libspsr/Src/MAT.inl
-file:/src/other/libspsr/Src/MarchingCubes.cpp
-file:/src/other/libspsr/Src/MarchingCubes.h
-file:/src/other/libspsr/Src/MemoryUsage.h
-file:/src/other/libspsr/Src/MultiGridOctreeData.IsoSurface.inl
-file:/src/other/libspsr/Src/MultiGridOctreeData.SortedTreeNodes.inl
-file:/src/other/libspsr/Src/MultiGridOctreeData.h
-file:/src/other/libspsr/Src/MultiGridOctreeData.inl
-file:/src/other/libspsr/Src/MyTime.h
-file:/src/other/libspsr/Src/Octree.h
-file:/src/other/libspsr/Src/Octree.inl
-file:/src/other/libspsr/Src/PPolynomial.h
-file:/src/other/libspsr/Src/PPolynomial.inl
-file:/src/other/libspsr/Src/PlyVertexMini.h
-file:/src/other/libspsr/Src/PointStream.h
-file:/src/other/libspsr/Src/PointStream.inl
-file:/src/other/libspsr/Src/PoissonRecon.cpp
-file:/src/other/libspsr/Src/Polynomial.h
-file:/src/other/libspsr/Src/Polynomial.inl
-file:/src/other/libspsr/Src/SPSR.cpp
-file:/src/other/libspsr/Src/SPSR.h
-file:/src/other/libspsr/Src/SparseMatrix.h
-file:/src/other/libspsr/Src/SparseMatrix.inl
-file:/src/other/libspsr/Src/SurfaceTrimmer.cpp
-file:/src/other/libspsr/Src/Vector.h
-file:/src/other/libspsr/Src/Vector.inl
-file:/src/other/libspsr/Src/cvertex.h
+file:/src/libbg/spsr/Allocator.h
+file:/src/libbg/spsr/Array.h
+file:/src/libbg/spsr/Array.inl
+file:/src/libbg/spsr/BSplineData.h
+file:/src/libbg/spsr/BSplineData.inl
+file:/src/libbg/spsr/BinaryNode.h
+file:/src/libbg/spsr/Factor.cpp
+file:/src/libbg/spsr/Factor.h
+file:/src/libbg/spsr/FunctionData.h
+file:/src/libbg/spsr/FunctionData.inl
+file:/src/libbg/spsr/Geometry.cpp
+file:/src/libbg/spsr/Geometry.h
+file:/src/libbg/spsr/Geometry.inl
+file:/src/libbg/spsr/Hash.h
+file:/src/libbg/spsr/MAT.h
+file:/src/libbg/spsr/MAT.inl
+file:/src/libbg/spsr/MarchingCubes.cpp
+file:/src/libbg/spsr/MarchingCubes.h
+file:/src/libbg/spsr/MemoryUsage.h
+file:/src/libbg/spsr/MultiGridOctreeData.IsoSurface.inl
+file:/src/libbg/spsr/MultiGridOctreeData.SortedTreeNodes.inl
+file:/src/libbg/spsr/MultiGridOctreeData.h
+file:/src/libbg/spsr/MultiGridOctreeData.inl
+file:/src/libbg/spsr/Octree.h
+file:/src/libbg/spsr/Octree.inl
+file:/src/libbg/spsr/PPolynomial.h
+file:/src/libbg/spsr/PPolynomial.inl
+file:/src/libbg/spsr/PlyVertexMini.h
+file:/src/libbg/spsr/PointStream.h
+file:/src/libbg/spsr/PointStream.inl
+file:/src/libbg/spsr/Polynomial.h
+file:/src/libbg/spsr/Polynomial.inl
+file:/src/libbg/spsr/SPSR.cpp
+file:/src/libbg/spsr/SPSR.h
+file:/src/libbg/spsr/SparseMatrix.h
+file:/src/libbg/spsr/SparseMatrix.inl
+file:/src/libbg/spsr/Vector.h
+file:/src/libbg/spsr/Vector.inl
+file:/src/libbg/spsr/cvertex.h
Modified: brlcad/trunk/misc/CMake/BRLCAD_Targets.cmake
===================================================================
--- brlcad/trunk/misc/CMake/BRLCAD_Targets.cmake 2020-10-01 19:56:41 UTC
(rev 77310)
+++ brlcad/trunk/misc/CMake/BRLCAD_Targets.cmake 2020-10-01 20:42:02 UTC
(rev 77311)
@@ -121,9 +121,9 @@
endwhile(NOT "${f_ext}" STREQUAL "")
# C++
- if(${srcfile_ext} STREQUAL ".cxx" OR ${srcfile_ext} STREQUAL ".cpp" OR
${srcfile_ext} STREQUAL ".cc")
+ if(${srcfile_ext} STREQUAL ".cxx" OR ${srcfile_ext} STREQUAL ".cpp" OR
${srcfile_ext} STREQUAL ".cc" OR ${srcfile_ext} STREQUAL ".inl")
set(file_language CXX)
- endif(${srcfile_ext} STREQUAL ".cxx" OR ${srcfile_ext} STREQUAL ".cpp" OR
${srcfile_ext} STREQUAL ".cc")
+ endif(${srcfile_ext} STREQUAL ".cxx" OR ${srcfile_ext} STREQUAL ".cpp" OR
${srcfile_ext} STREQUAL ".cc" OR ${srcfile_ext} STREQUAL ".inl")
if(${srcfile_ext} STREQUAL ".hxx" OR ${srcfile_ext} STREQUAL ".hpp" OR
${srcfile_ext} STREQUAL ".hh")
set(file_language CXX)
endif(${srcfile_ext} STREQUAL ".hxx" OR ${srcfile_ext} STREQUAL ".hpp" OR
${srcfile_ext} STREQUAL ".hh")
Modified: brlcad/trunk/src/external/Creo/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/external/Creo/CMakeLists.txt 2020-10-01 19:56:41 UTC
(rev 77310)
+++ brlcad/trunk/src/external/Creo/CMakeLists.txt 2020-10-01 20:42:02 UTC
(rev 77311)
@@ -207,7 +207,6 @@
libnmg
librt
libwdb
- SPSR
gdiam
openNURBS
poly2tri
@@ -295,7 +294,6 @@
configure_file(creo-brl.dat.in
${CMAKE_CURRENT_BINARY_DIR}/${CFG_TYPE}/creo-brl-debug.dat @ONLY)
endforeach(CFG_TYPE ${CMAKE_CONFIGURATION_TYPES})
set_property(TARGET creo-brl-exe APPEND PROPERTY COMPILE_DEFINITIONS
"TINYCTHREAD_DLL_IMPORTS")
- set_property(TARGET creo-brl-exe APPEND PROPERTY COMPILE_DEFINITIONS
"SPSR_DLL_IMPORTS")
endif(DEFINED CREO_EXEC_PLUGIN)
# Install the resource files
Modified: brlcad/trunk/src/libbg/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/libbg/CMakeLists.txt 2020-10-01 19:56:41 UTC (rev
77310)
+++ brlcad/trunk/src/libbg/CMakeLists.txt 2020-10-01 20:42:02 UTC (rev
77311)
@@ -8,6 +8,47 @@
)
BRLCAD_LIB_INCLUDE_DIRS(bg BG_INCLUDE_DIRS "")
+set(SPSR_srcs
+ spsr/Allocator.h
+ spsr/Array.h
+ spsr/Array.inl
+ spsr/BinaryNode.h
+ spsr/BSplineData.h
+ spsr/BSplineData.inl
+ spsr/cvertex.h
+ spsr/Factor.cpp
+ spsr/Factor.h
+ spsr/FunctionData.h
+ spsr/FunctionData.inl
+ spsr/Geometry.cpp
+ spsr/Geometry.h
+ spsr/Geometry.inl
+ spsr/Hash.h
+ spsr/MarchingCubes.cpp
+ spsr/MarchingCubes.h
+ spsr/MAT.h
+ spsr/MAT.inl
+ spsr/MultiGridOctreeData.h
+ spsr/MultiGridOctreeData.inl
+ spsr/MultiGridOctreeData.IsoSurface.inl
+ spsr/MultiGridOctreeData.SortedTreeNodes.inl
+ spsr/Octree.h
+ spsr/Octree.inl
+ spsr/PlyVertexMini.h
+ spsr/PointStream.h
+ spsr/PointStream.inl
+ spsr/Polynomial.h
+ spsr/Polynomial.inl
+ spsr/PPolynomial.h
+ spsr/PPolynomial.inl
+ spsr/SparseMatrix.h
+ spsr/SparseMatrix.inl
+ spsr/SPSR.cpp
+ spsr/SPSR.h
+ spsr/Vector.h
+ spsr/Vector.inl
+ )
+
set(LIBBG_SOURCES
aabb_ray.c
chull.c
@@ -23,6 +64,7 @@
polygon_triangulate.cpp
polygon_point_in.c
spsr.c
+ ${SPSR_srcs}
tri_pt.c
tri_ray.c
tri_tri.c
@@ -31,15 +73,13 @@
util.c
)
-BRLCAD_ADDLIB(libbg "${LIBBG_SOURCES}" "libbn;libbu;SPSR;${P2T_LIBRARY}")
+BRLCAD_ADDLIB(libbg "${LIBBG_SOURCES}" "libbn;libbu;${P2T_LIBRARY}")
set_target_properties(libbg PROPERTIES VERSION 20.0.1 SOVERSION 20)
if (HIDE_INTERNAL_SYMBOLS)
- set_property(TARGET libbg APPEND PROPERTY COMPILE_DEFINITIONS
"SPSR_DLL_IMPORTS")
if (TARGET poly2tri OR HIDE_INTERNAL_SYMBOLS_EXT)
set_property(TARGET libbg APPEND PROPERTY COMPILE_DEFINITIONS
"P2T_DLL_IMPORTS")
endif (TARGET poly2tri OR HIDE_INTERNAL_SYMBOLS_EXT)
if (TARGET libbg-obj)
- set_property(TARGET libbg-obj APPEND PROPERTY COMPILE_DEFINITIONS
"SPSR_DLL_IMPORTS")
if (TARGET poly2tri OR HIDE_INTERNAL_SYMBOLS_EXT)
set_property(TARGET libbg-obj APPEND PROPERTY COMPILE_DEFINITIONS
"P2T_DLL_IMPORTS")
endif (TARGET poly2tri OR HIDE_INTERNAL_SYMBOLS_EXT)
@@ -55,6 +95,7 @@
delaunator.hpp
earcut.hpp
pointgen.c
+ spsr/LICENSE
)
CMAKEFILES(${bg_ignore})
@@ -66,3 +107,4 @@
# indent-tabs-mode: t
# End:
# ex: shiftwidth=2 tabstop=8
+
Added: brlcad/trunk/src/libbg/spsr/Allocator.h
===================================================================
--- brlcad/trunk/src/libbg/spsr/Allocator.h (rev 0)
+++ brlcad/trunk/src/libbg/spsr/Allocator.h 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,163 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+#ifndef ALLOCATOR_INCLUDED
+#define ALLOCATOR_INCLUDED
+#include <vector>
+
+class AllocatorState{
+public:
+ int index,remains;
+};
+/** This templated class assists in memory allocation and is well suited for
instances
+ * when it is known that the sequence of memory allocations is performed in a
stack-based
+ * manner, so that memory allocated last is released first. It also
preallocates memory
+ * in chunks so that multiple requests for small chunks of memory do not
require separate
+ * system calls to the memory manager.
+ * The allocator is templated off of the class of objects that we would like
it to allocate,
+ * ensuring that appropriate constructors and destructors are called as
necessary.
+ */
+template<class T>
+class Allocator
+{
+ int blockSize;
+ int index , remains;
+ std::vector< T* > memory;
+public:
+ Allocator( void ){ blockSize = index = remains = 0; }
+ ~Allocator( void ){ reset(); }
+
+ /** This method is the allocators destructor. It frees up any of the
memory that
+ * it has allocated. */
+ void reset( void )
+ {
+ for(size_t i=0;i<memory.size();i++){delete[] memory[i];}
+ memory.clear();
+ blockSize=index=remains=0;
+ }
+ /** This method returns the memory state of the allocator. */
+ AllocatorState getState( void ) const
+ {
+ AllocatorState s;
+ s.index=index;
+ s.remains=remains;
+ return s;
+ }
+
+
+ /** This method rolls back the allocator so that it makes all of the
memory previously
+ * allocated available for re-allocation. Note that it does it not
call the constructor
+ * again, so after this method has been called, assumptions about the
state of the values
+ * in memory are no longer valid. */
+ void rollBack(void)
+ {
+ if( memory.size() )
+ {
+ for( size_t i=0 ; i<memory.size() ; i++ )
+ {
+ for( int j=0 ; j<blockSize ; j++ )
+ {
+ memory[i][j].~T();
+ new(&memory[i][j]) T();
+ }
+ }
+ index=0;
+ remains=blockSize;
+ }
+ }
+ /** This method rolls back the allocator to the previous memory state
and makes all of the memory previously
+ * allocated available for re-allocation. Note that it does it not
call the constructor
+ * again, so after this method has been called, assumptions about the
state of the values
+ * in memory are no longer valid. */
+ void rollBack(const AllocatorState& state){
+ if(state.index<index || (state.index==index &&
state.remains<remains)){
+ if(state.index<index){
+ for(int j=state.remains;j<blockSize;j++){
+ memory[state.index][j].~T();
+ new(&memory[state.index][j]) T();
+ }
+ for(int i=state.index+1;i<index-1;i++){
+ for(int j=0;j<blockSize;j++){
+ memory[i][j].~T();
+ new(&memory[i][j]) T();
+ }
+ }
+ for(int j=0;j<remains;j++){
+ memory[index][j].~T();
+ new(&memory[index][j]) T();
+ }
+ index=state.index;
+ remains=state.remains;
+ }
+ else{
+ for(int j=0;j<state.remains;j<remains){
+ memory[index][j].~T();
+ new(&memory[index][j]) T();
+ }
+ remains=state.remains;
+ }
+ }
+ }
+
+ /** This method initializes the constructor and the blockSize variable
specifies the
+ * the number of objects that should be pre-allocated at a time. */
+ void set( int blockSize)
+ {
+ reset();
+ this->blockSize = blockSize;
+ index=-1;
+ remains=0;
+ }
+
+ /** This method returns a pointer to an array of elements objects. If
there is left over pre-allocated
+ * memory, this method simply returns a pointer to the next free piece
of memory, otherwise it pre-allocates
+ * more memory. Note that if the number of objects requested is larger
than the value blockSize with which
+ * the allocator was initialized, the request for memory will fail.
+ */
+ T* newElements( int elements=1 )
+ {
+ T* mem;
+ if( !elements ) return NULL;
+ if( elements>blockSize ) fprintf( stderr , "[ERROR] Allocator:
elements bigger than block-size: %d>%d\n" , elements , blockSize ) , exit( 0 );
+ if( remains<elements )
+ {
+ if( index==memory.size()-1 )
+ {
+ mem = new T[blockSize];
+ if( !mem ) fprintf( stderr , "[ERROR] Failed to
allocate memory\n" ) , exit(0);
+ memory.push_back( mem );
+ }
+ index++;
+ remains=blockSize;
+ }
+ mem = &(memory[index][blockSize-remains]);
+ remains -= elements;
+ return mem;
+ }
+};
+#endif // ALLOCATOR_INCLUDE
Property changes on: brlcad/trunk/src/libbg/spsr/Allocator.h
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/Array.h
===================================================================
--- brlcad/trunk/src/libbg/spsr/Array.h (rev 0)
+++ brlcad/trunk/src/libbg/spsr/Array.h 2020-10-01 20:42:02 UTC (rev 77311)
@@ -0,0 +1,97 @@
+/*
+Copyright (c) 2011, Michael Kazhdan and Ming Chuang
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+#ifndef ARRAY_INCLUDED
+#define ARRAY_INCLUDED
+
+#include <vector>
+
+#define ARRAY_DEBUG 0
+
+// Code from http://stackoverflow.com
+void* aligned_malloc( size_t size , size_t align )
+{
+ // Align enough for the data, the alignment padding, and room to store
a pointer to the actual start of the memory
+ void* mem = malloc( size + align + sizeof( void* ) );
+ // The position at which we could potentially start addressing
+ char* amem = ( (char*)mem ) + sizeof( void* );
+ // Add align-1 to the start of the address and then zero out at most of
the first align-1 bits.
+ amem = ( char* )( ( (size_t)( ( (char*)amem ) + (align-1) ) ) & ~(
align-1 ) );
+ // Pre-write the actual address
+ ( ( void** ) amem )[-1] = mem;
+ return amem;
+}
+void aligned_free( void* mem ) { free( ( ( void** )mem )[-1] ); }
+
+#if ARRAY_DEBUG
+#pragma message ( "[WARNING] Array debugging is enabled" )
+#include "Array.inl"
+#define Pointer( ... ) Array< __VA_ARGS__ >
+#define ConstPointer( ... ) ConstArray< __VA_ARGS__ >
+template< class C > void FreePointer( Array< C >& a ){ a.Free( ); }
+template< class C > void AlignedFreePointer( Array< C >& a ){ a.Free( ); }
+template< class C > void VFreePointer( Array< C >& a ){ a.Free( ); }
+template< class C > void DeletePointer( Array< C >& a ){ a.Delete( ); }
+
+template< class C > Array< C > NewPointer( size_t
size , const char* name=NULL ){ return Array< C >::New
( size , name ); }
+template< class C > Array< C > AllocPointer( size_t
size , const char* name=NULL ){ return Array< C >::Alloc
( size , false , name ); }
+template< class C > Array< C > AlignedAllocPointer( size_t
size , size_t alignment , const char* name=NULL ){ return Array< C
>::AlignedAlloc( size , alignment , false , name ); }
+template< class C > Array< C > ReAllocPointer( Array< C >& a , size_t
size , const char* name=NULL ){ return Array< C >::ReAlloc
( a , size , false , name ); }
+
+template< class C > Array< C > NullPointer( void ){ return Array< C >( ); }
+
+template< class C > C* PointerAddress( Array< C >& a ) { return
a.pointer(); }
+template< class C > const C* PointerAddress( ConstArray< C >& a ) { return
a.pointer(); }
+template< class C > Array< C > GetPointer( C& c ) { return
Array< C >::FromPointer( &c , 1 ); }
+template< class C > ConstArray< C > GetPointer( const C& c ) { return
ConstArray< C >::FromPointer( &c , 1 ); }
+template< class C > Array< C > GetPointer( std::vector< C >& v ){
return Array< C >::FromPointer( &v[0] , v.size() ); }
+template< class C > ConstArray< C > GetPointer( const std::vector< C >& v ){
return ConstArray< C >::FromPointer( &v[0] , v.size() ); }
+
+#else // !ARRAY_DEBUG
+#define Pointer( ... ) __VA_ARGS__*
+#define ConstPointer( ... ) const __VA_ARGS__*
+
+#define FreePointer( ... ) { if( __VA_ARGS__ ) free(
__VA_ARGS__ ) , __VA_ARGS__ = NULL; }
+#define AlignedFreePointer( ... ) { if( __VA_ARGS__ ) aligned_free(
__VA_ARGS__ ) , __VA_ARGS__ = NULL; }
+#define DeletePointer( ... ) { if( __VA_ARGS__ ) delete[]
__VA_ARGS__ , __VA_ARGS__ = NULL; }
+
+template< class C > C* NewPointer( size_t size ,
const char* name=NULL ){ return new C[size]; }
+template< class C > C* AllocPointer( size_t size ,
const char* name=NULL ){ return (C*) malloc( sizeof(C) *
size ); }
+template< class C > C* AlignedAllocPointer( size_t size , size_t
alignment , const char* name=NULL ){ return (C*)aligned_malloc(
sizeof(C) * size , alignment ); }
+template< class C > C* ReAllocPointer( C* c , size_t size ,
const char* name=NULL ){ return (C*) realloc( c , sizeof(C) *
size ); }
+
+template< class C > C* NullPointer( void ){ return NULL; }
+
+template< class C > C* PointerAddress( C* c ){ return c; }
+template< class C > const C* PointerAddress( const C* c ){ return c; }
+template< class C > C* GetPointer( C& c ){ return &c; }
+template< class C > const C* GetPointer( const C& c ){ return &c; }
+template< class C > C* GetPointer( std::vector< C >& v ){ return
&v[0]; }
+template< class C > const C* GetPointer( const std::vector< C >& v ){ return
&v[0]; }
+#endif // ARRAY_DEBUG
+#endif // ARRAY_INCLUDED
Property changes on: brlcad/trunk/src/libbg/spsr/Array.h
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/Array.inl
===================================================================
--- brlcad/trunk/src/libbg/spsr/Array.inl (rev 0)
+++ brlcad/trunk/src/libbg/spsr/Array.inl 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,655 @@
+/*
+Copyright (c) 2011, Michael Kazhdan and Ming Chuang
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+#define FULL_ARRAY_DEBUG 0 // Note that this is not thread-safe
+
+#include <stdio.h>
+#include <emmintrin.h>
+#include <vector>
+#include <stddef.h>
+
+inline bool isfinitef( float fp ){ float f=fp; return ((*(unsigned
*)&f)&0x7f800000)!=0x7f800000; }
+
+
+template< class C > bool IsValid( const C& c );
+#if _DEBUG
+template< > inline bool IsValid< float >( const float& f ) { return
isfinitef( f ) && ( f==0.f || abs(f)>1e-31f ); }
+#else // !_DEBUG
+template< > inline bool IsValid< float >( const float& f ) { return
isfinitef( f ); }
+#endif // _DEBUG
+template< > inline bool IsValid< __m128 >( const __m128& m )
+{
+ const __m128* addr = &m;
+ if( size_t(addr) & 15 ) return false;
+ else return true;
+}
+template< class C > inline bool IsValid( const C& c ){ return true; }
+
+
+#if FULL_ARRAY_DEBUG
+class DebugMemoryInfo
+{
+public:
+ const void* address;
+ char name[512];
+};
+static std::vector< DebugMemoryInfo > memoryInfo;
+#endif // FULL_ARRAY_DEBUG
+
+template< class C >
+class Array
+{
+ void _assertBounds( long long idx ) const
+ {
+ if( idx<min || idx>=max )
+ {
+ fprintf( stderr , "Array index out-of-bounds: %lld <=
%lld < %lld\n" , min , idx , max );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+protected:
+ C *data , *_data;
+ long long min , max;
+#if FULL_ARRAY_DEBUG
+ static void _AddMemoryInfo( const void* ptr , const char* name )
+ {
+ size_t sz = memoryInfo.size();
+ memoryInfo.resize( sz + 1 );
+ memoryInfo[sz].address = ptr;
+ if( name ) strcpy( memoryInfo[sz].name , name );
+ else memoryInfo[sz].name[0] = 0;
+ }
+ static void _RemoveMemoryInfo( const void* ptr )
+ {
+ {
+ size_t idx;
+ for( idx=0 ; idx<memoryInfo.size( ) ; idx++ ) if(
memoryInfo[idx].address==ptr ) break;
+ if( idx==memoryInfo.size() )
+ {
+ fprintf( stderr , "Could not find memory in
address table\n" );
+ ASSERT( 0 );
+ }
+ else
+ {
+ memoryInfo[idx] =
memoryInfo[memoryInfo.size()-1];
+ memoryInfo.pop_back( );
+ }
+ }
+ }
+#endif // FULL_ARRAY_DEBUG
+
+public:
+ long long minimum( void ) const { return min; }
+ long long maximum( void ) const { return max; }
+
+ static inline Array New( size_t size , const char* name=NULL )
+ {
+ Array a;
+ a._data = a.data = new C[size];
+ a.min = 0;
+#pragma message( "[WARNING] Casting unsigned to signed" )
+ a.max = ( long long ) size;
+#if FULL_ARRAY_DEBUG
+ _AddMemoryInfo( a._data , name );
+#endif // FULL_ARRAY_DEBUG
+ return a;
+ }
+ static inline Array Alloc( size_t size , bool clear , const char*
name=NULL )
+ {
+ Array a;
+ a._data = a.data = ( C* ) malloc( size * sizeof( C ) );
+ if( clear ) memset( a.data , 0 , size * sizeof( C ) );
+// else memset( a.data , -1 , size * sizeof( C ) );
+ a.min = 0;
+#pragma message( "[WARNING] Casting unsigned to signed" )
+ a.max = ( long long ) size;
+#if FULL_ARRAY_DEBUG
+ _AddMemoryInfo( a._data , name );
+#endif // FULL_ARRAY_DEBUG
+ return a;
+ }
+ static inline Array AlignedAlloc( size_t size , size_t alignment , bool
clear , const char* name=NULL )
+ {
+ Array a;
+ a.data = ( C* ) aligned_malloc( sizeof(C) * size , alignment );
+ a._data = ( C* )( ( ( void** )a.data )[-1] );
+ if( clear ) memset( a.data , 0 , size * sizeof( C ) );
+// else memset( a.data , -1 , size * sizeof( C ) );
+ a.min = 0;
+#pragma message( "[WARNING] Casting unsigned to signed" )
+ a.max = ( long long ) size;
+#if FULL_ARRAY_DEBUG
+ _AddMemoryInfo( a._data , name );
+#endif // FULL_ARRAY_DEBUG
+ return a;
+ }
+ static inline Array ReAlloc( Array& a , size_t size , bool clear ,
const char* name=NULL )
+ {
+ Array _a;
+ _a._data = _a.data = ( C* ) realloc( a.data , size * sizeof( C
) );
+ if( clear ) memset( _a.data , 0 , size * sizeof( C ) );
+#if FULL_ARRAY_DEBUG
+ _RemoveMemoryInfo( a._data );
+#endif // FULL_ARRAY_DEBUG
+ a._data = NULL;
+ _a.min = 0;
+#pragma message( "[WARNING] Casting unsigned to signed" )
+ _a.max = ( long long ) size;
+#if FULL_ARRAY_DEBUG
+ _AddMemoryInfo( _a._data , name );
+#endif // FULL_ARRAY_DEBUG
+ return _a;
+ }
+
+ Array( void )
+ {
+ data = _data = NULL;
+ min = max = 0;
+ }
+ template< class D >
+ Array( Array< D >& a )
+ {
+ _data = NULL;
+ if( !a )
+ {
+ data = NULL;
+ min = max = 0;
+ }
+ else
+ {
+ // [WARNING] Changing szC and szD to size_t causes some
really strange behavior.
+ long long szC = sizeof( C );
+ long long szD = sizeof( D );
+ data = (C*)&a[0];
+ min = ( a.minimum() * szD ) / szC;
+ max = ( a.maximum() * szD ) / szC;
+ if( min*szC!=a.minimum()*szD ||
max*szC!=a.maximum()*szD )
+ {
+ fprintf( stderr , "Could not convert array [
%lld , %lld ] * %lld => [ %lld , %lld ] * %lld\n" , a.minimum() , a.maximum() ,
szD , min , max , szC );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+ }
+ static Array FromPointer( C* data , long long max )
+ {
+ Array a;
+ a._data = NULL;
+ a.data = data;
+ a.min = 0;
+ a.max = max;
+ return a;
+ }
+ static Array FromPointer( C* data , long long min , long long max )
+ {
+ Array a;
+ a._data = NULL;
+ a.data = data;
+ a.min = min;
+ a.max = max;
+ return a;
+ }
+ inline bool operator == ( const Array< C >& a ) const { return
data==a.data; }
+ inline bool operator != ( const Array< C >& a ) const { return
data!=a.data; }
+ inline bool operator == ( const C* c ) const { return data==c; }
+ inline bool operator != ( const C* c ) const { return data!=c; }
+ inline C* operator -> ( void )
+ {
+ _assertBounds( 0 );
+ return data;
+ }
+ inline const C* operator -> ( ) const
+ {
+ _assertBounds( 0 );
+ return data;
+ }
+ inline C& operator[]( long long idx )
+ {
+ _assertBounds( idx );
+ return data[idx];
+ }
+ inline const C& operator[]( long long idx ) const
+ {
+ _assertBounds( idx );
+ return data[idx];
+ }
+ inline Array operator + ( int idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array operator + ( long long idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array operator + ( unsigned int idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array operator + ( unsigned long long idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array& operator += ( int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator += ( long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator += ( unsigned int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator += ( unsigned long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator ++ ( void ) { return (*this) += 1; }
+ Array operator - ( int idx ) const { return (*this) + (-idx); }
+ Array operator - ( long long idx ) const { return (*this) + (-idx); }
+ Array operator - ( unsigned int idx ) const { return (*this) +
(-idx); }
+ Array operator - ( unsigned long long idx ) const { return (*this) +
(-idx); }
+ Array& operator -= ( int idx ) { return (*this) += (-idx); }
+ Array& operator -= ( long long idx ) { return (*this) += (-idx); }
+ Array& operator -= ( unsigned int idx ) { return (*this) += (-idx); }
+ Array& operator -= ( unsigned long long idx ) { return (*this) +=
(-idx); }
+ Array& operator -- ( void ) { return (*this) -= 1; }
+ long long operator - ( const Array& a ) const { return ( long long )(
data - a.data ); }
+
+ void Free( void )
+ {
+ if( _data )
+ {
+ free( _data );
+#if FULL_ARRAY_DEBUG
+ _RemoveMemoryInfo( _data );
+#endif // FULL_ARRAY_DEBUG
+ }
+ (*this) = Array( );
+ }
+ void Delete( void )
+ {
+ if( _data )
+ {
+ delete[] _data;
+#if FULL_ARRAY_DEBUG
+ _RemoveMemoryInfo( _data );
+#endif // FULL_ARRAY_DEBUG
+ }
+ (*this) = Array( );
+ }
+ C* pointer( void ){ return data; }
+ const C* pointer( void ) const { return data; }
+ bool operator !( void ) const { return data==NULL; }
+ operator bool( ) const { return data!=NULL; }
+};
+
+template< class C >
+class ConstArray
+{
+ void _assertBounds( long long idx ) const
+ {
+ if( idx<min || idx>=max )
+ {
+ fprintf( stderr , "ConstArray index out-of-bounds: %lld
<= %lld < %lld\n" , min , idx , max );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+protected:
+ const C *data;
+ long long min , max;
+public:
+ long long minimum( void ) const { return min; }
+ long long maximum( void ) const { return max; }
+
+ inline ConstArray( void )
+ {
+ data = NULL;
+ min = max = 0;
+ }
+ inline ConstArray( const Array< C >& a )
+ {
+ // [WARNING] Changing szC and szD to size_t causes some really
strange behavior.
+ data = ( const C* )a.pointer( );
+ min = a.minimum();
+ max = a.maximum();
+ }
+ template< class D >
+ inline ConstArray( const Array< D >& a )
+ {
+ // [WARNING] Changing szC and szD to size_t causes some really
strange behavior.
+ long long szC = ( long long ) sizeof( C );
+ long long szD = ( long long ) sizeof( D );
+ data = ( const C* )a.pointer( );
+ min = ( a.minimum() * szD ) / szC;
+ max = ( a.maximum() * szD ) / szC;
+ if( min*szC!=a.minimum()*szD || max*szC!=a.maximum()*szD )
+ {
+// fprintf( stderr , "Could not convert const array [ %lld
, %lld ] * %lld => [ %lld , %lld ] * %lld\n" , a.minimum() , a.maximum() , szD
, min , max , szC );
+ fprintf( stderr , "Could not convert const array [ %lld
, %lld ] * %lld => [ %lld , %lld ] * %lld\n %lld %lld %lld\n" , a.minimum() ,
a.maximum() , szD , min , max , szC , a.minimum() , a.minimum()*szD ,
(a.minimum()*szD)/szC );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+ template< class D >
+ inline ConstArray( const ConstArray< D >& a )
+ {
+ // [WARNING] Changing szC and szD to size_t causes some really
strange behavior.
+ long long szC = sizeof( C );
+ long long szD = sizeof( D );
+ data = ( const C*)a.pointer( );
+ min = ( a.minimum() * szD ) / szC;
+ max = ( a.maximum() * szD ) / szC;
+ if( min*szC!=a.minimum()*szD || max*szC!=a.maximum()*szD )
+ {
+ fprintf( stderr , "Could not convert array [ %lld ,
%lld ] * %lld => [ %lld , %lld ] * %lld\n" , a.minimum() , a.maximum() , szD ,
min , max , szC );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+ static ConstArray FromPointer( const C* data , long long max )
+ {
+ ConstArray a;
+ a.data = data;
+ a.min = 0;
+ a.max = max;
+ return a;
+ }
+ static ConstArray FromPointer( const C* data , long long min , long
long max )
+ {
+ ConstArray a;
+ a.data = data;
+ a.min = min;
+ a.max = max;
+ return a;
+ }
+
+ inline bool operator == ( const ConstArray< C >& a ) const { return
data==a.data; }
+ inline bool operator != ( const ConstArray< C >& a ) const { return
data!=a.data; }
+ inline bool operator == ( const C* c ) const { return data==c; }
+ inline bool operator != ( const C* c ) const { return data!=c; }
+ inline const C* operator -> ( void )
+ {
+ _assertBounds( 0 );
+ return data;
+ }
+ inline const C& operator[]( long long idx ) const
+ {
+ _assertBounds( idx );
+ return data[idx];
+ }
+ inline ConstArray operator + ( int idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray operator + ( long long idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray operator + ( unsigned int idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray operator + ( unsigned long long idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray& operator += ( int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator += ( long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator += ( unsigned int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator += ( unsigned long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator ++ ( void ) { return (*this) += 1; }
+ ConstArray operator - ( int idx ) const { return (*this) + (-idx); }
+ ConstArray operator - ( long long idx ) const { return (*this) +
(-idx); }
+ ConstArray operator - ( unsigned int idx ) const { return (*this) +
(-idx); }
+ ConstArray operator - ( unsigned long long idx ) const { return
(*this) + (-idx); }
+ ConstArray& operator -= ( int idx ) { return (*this) += (-idx); }
+ ConstArray& operator -= ( long long idx ) { return (*this) +=
(-idx); }
+ ConstArray& operator -= ( unsigned int idx ) { return (*this) +=
(-idx); }
+ ConstArray& operator -= ( unsigned long long idx ) { return (*this)
+= (-idx); }
+ ConstArray& operator -- ( void ) { return (*this) -= 1; }
+ long long operator - ( const ConstArray& a ) const { return ( long long
)( data - a.data ); }
+ long long operator - ( const Array< C >& a ) const { return ( long long
)( data - a.pointer() ); }
+
+ const C* pointer( void ) const { return data; }
+ bool operator !( void ) { return data==NULL; }
+ operator bool( ) { return data!=NULL; }
+};
+
+#if FULL_ARRAY_DEBUG
+inline void PrintMemoryInfo( void ){ for( size_t i=0 ; i<memoryInfo.size() ;
i++ ) printf( "%d] %s\n" , i , memoryInfo[i].name ); }
+#endif // FULL_ARRAY_DEBUG
+template< class C >
+Array< C > memcpy( Array< C > destination , const void* source , size_t size )
+{
+ if( size>destination.maximum()*sizeof(C) )
+ {
+ fprintf( stderr , "Size of copy exceeds destination maximum:
%lld > %lld\n" , ( long long )( size ) , ( long long )(
destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( &destination[0] , source , size );
+ return destination;
+}
+template< class C , class D >
+Array< C > memcpy( Array< C > destination , Array< D > source , size_t size )
+{
+ if( size>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of copy exceeds destination maximum:
%lld > %lld\n" , ( long long )( size ) , ( long long )(
destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld >
%lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) )
);
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( &destination[0] , &source[0] , size );
+ return destination;
+}
+template< class C , class D >
+Array< C > memcpy( Array< C > destination , ConstArray< D > source , size_t
size )
+{
+ if( size>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of copy exceeds destination maximum:
%lld > %lld\n" , ( long long )( size ) , ( long long )(
destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld >
%lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) )
);
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( &destination[0] , &source[0] , size );
+ return destination;
+}
+template< class D >
+void* memcpy( void* destination , Array< D > source , size_t size )
+{
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld >
%lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) )
);
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( destination , &source[0] , size );
+ return destination;
+}
+template< class D >
+void* memcpy( void* destination , ConstArray< D > source , size_t size )
+{
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld >
%lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) )
);
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( destination , &source[0] , size );
+ return destination;
+}
+template< class C >
+Array< C > memset( Array< C > destination , int value , size_t size )
+{
+ if( size>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of set exceeds destination maximum:
%lld > %lld\n" , ( long long )( size ) , ( long long )(
destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memset( &destination[0] , value , size );
+ return destination;
+}
+
+template< class C >
+size_t fread( Array< C > destination , size_t eSize , size_t count , FILE* fp )
+{
+ if( count*eSize>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of read exceeds source maximum: %lld >
%lld\n" , ( long long )( count*eSize ) , ( long long )(
destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ return fread( &destination[0] , eSize , count , fp );
+}
+template< class C >
+size_t fwrite( Array< C > source , size_t eSize , size_t count , FILE* fp )
+{
+ if( count*eSize>source.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of write exceeds source maximum: %lld >
%lld\n" , ( long long )( count*eSize ) , ( long long )(
source.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ return fwrite( &source[0] , eSize , count , fp );
+}
+template< class C >
+size_t fwrite( ConstArray< C > source , size_t eSize , size_t count , FILE* fp
)
+{
+ if( count*eSize>source.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of write exceeds source maximum: %lld >
%lld\n" , ( long long )( count*eSize ) , ( long long )(
source.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ return fwrite( &source[0] , eSize , count , fp );
+}
+template< class C >
+void qsort( Array< C > base , size_t numElements , size_t elementSize , int
(*compareFunction)( const void* , const void* ) )
+{
+ if( sizeof(C)!=elementSize )
+ {
+ fprintf( stderr , "Element sizes differ: %lld != %lld\n" , (
long long )( sizeof(C) ) , ( long long )( elementSize ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( base.minimum()>0 || base.maximum()<numElements )
+ {
+ fprintf( stderr , "Array access out of bounds: %lld <= 0 <=
%lld <= %lld\n" , base.minimum() , base.maximum() , ( long long )( numElements
) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ qsort( base.pointer() , numElements , elementSize , compareFunction );
+}
Property changes on: brlcad/trunk/src/libbg/spsr/Array.inl
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/BSplineData.h
===================================================================
--- brlcad/trunk/src/libbg/spsr/BSplineData.h (rev 0)
+++ brlcad/trunk/src/libbg/spsr/BSplineData.h 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,191 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+#ifndef BSPLINE_DATA_INCLUDED
+#define BSPLINE_DATA_INCLUDED
+
+#include "PPolynomial.h"
+#include "Array.h"
+
+template< int Degree >
+struct BSplineElementCoefficients
+{
+ int coeffs[Degree+1];
+ BSplineElementCoefficients( void ){ memset( coeffs , 0 , sizeof( int )
* ( Degree+1 ) ); }
+ int& operator[]( int idx ){ return coeffs[idx]; }
+ const int& operator[]( int idx ) const { return coeffs[idx]; }
+};
+template< int Degree >
+struct BSplineElements : public std::vector< BSplineElementCoefficients<
Degree > >
+{
+ static const int _off = (Degree+1)/2;
+ void _addLeft ( int offset , int boundary );
+ void _addRight( int offset , int boundary );
+public:
+ enum
+ {
+ NONE = 0,
+ DIRICHLET = -1,
+ NEUMANN = 1
+ };
+ // Coefficients are ordered as "/" "-" "\"
+ int denominator;
+
+ BSplineElements( void ) { denominator = 1; }
+ BSplineElements( int res , int offset , int boundary=NONE , int inset=0
);
+
+ void upSample( BSplineElements& high ) const;
+ void differentiate( BSplineElements< Degree-1 >& d ) const;
+
+ void print( FILE* fp=stdout ) const
+ {
+ for( int i=0 ; i<std::vector< BSplineElementCoefficients<
Degree > >::size() ; i++ )
+ {
+ printf( "%d]" , i );
+ for( int j=0 ; j<=Degree ; j++ ) printf( " %d" ,
(*this)[i][j] );
+ printf( " (%d)\n" , denominator );
+ }
+ }
+};
+
+template< int Degree >
+class BSplineData
+{
+ int _boundaryType;
+ double _vvIntegrals[Degree+1][Degree+1];
+ double _vdIntegrals[Degree+1][Degree ];
+ double _dvIntegrals[Degree ][Degree+1];
+ double _ddIntegrals[Degree ][Degree ];
+
+public:
+ struct Integrator
+ {
+ struct IntegralTables
+ {
+ double vv_ccIntegrals[2*Degree+1][2*Degree+1] ,
vv_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ double dv_ccIntegrals[2*Degree+1][2*Degree+1] ,
dv_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ double vd_ccIntegrals[2*Degree+1][2*Degree+1] ,
vd_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ double dd_ccIntegrals[2*Degree+1][2*Degree+1] ,
dd_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ };
+ std::vector< IntegralTables > iTables;
+ double dot( int depth , int off1 , int off2 , bool d1 , bool d2
, bool childParent=false ) const;
+ };
+ double dot( int depth1 , int off1 , int depth2 , int off2 , bool d1 ,
bool d2 , bool inset=false ) const;
+ void setIntegrator( Integrator& integrator , bool inset , bool
useDotRatios=false ) const;
+ template< int Radius >
+ struct CenterEvaluator
+ {
+ struct ValueTables
+ {
+ double vValues[2*Degree+1][ 3*(2*Radius+1) ];
+ double dValues[2*Degree+1][ 3*(2*Radius+1) ];
+ };
+ std::vector< ValueTables > vTables;
+ double value( int depth , int off1 , int off2 , bool d , bool
childParent=false ) const;
+ };
+ template< int Radius >
+ void setCenterEvaluator( CenterEvaluator< Radius >& evaluator , double
smoothingRadius , double dSmoothingRadius, bool inset ) const;
+ double value( int depth , int off , double smoothingRadius , double s ,
bool d , bool inset=false ) const;
+ template< int Radius >
+ struct CornerEvaluator
+ {
+ struct ValueTables
+ {
+ double vValues[2*Degree+1][4*Radius+3];
+ double dValues[2*Degree+1][4*Radius+3];
+ };
+ std::vector< ValueTables > vTables;
+ double value( int depth , int off1 , int c1 , int off2 , bool d
, bool childParent=false ) const;
+ };
+ template< int Radius >
+ void setCornerEvaluator( CornerEvaluator< Radius >& evaluator , double
smoothingRadius , double dSmoothingRadius, bool inset ) const;
+
+ struct BSplineComponents
+ {
+ Polynomial< Degree > polys[Degree+1];
+ Polynomial< Degree >& operator[] ( int idx ) { return
polys[idx]; }
+ const Polynomial< Degree >& operator[] ( int idx ) const {
return polys[idx]; }
+ void printnl( void ) const { for( int d=0 ; d<=Degree ; d++ )
polys[d].printnl(); }
+ BSplineComponents scale( double s ) const { BSplineComponents b
; for( int d=0 ; d<=Degree ; d++ ) b[d] = polys[d].scale(s) ; return b; }
+ BSplineComponents shift( double s ) const { BSplineComponents b
; for( int d=0 ; d<=Degree ; d++ ) b[d] = polys[d].shift(s) ; return b; }
+ };
+
+ int depth;
+ size_t functionCount , sampleCount;
+ PPolynomial< Degree > baseFunction , leftBaseFunction ,
rightBaseFunction , leftRightBaseFunction;
+ PPolynomial< Degree-1 > dBaseFunction , dLeftBaseFunction ,
dRightBaseFunction , dLeftRightBaseFunction;
+ BSplineComponents baseBSpline , leftBSpline , rightBSpline ,
leftRightBSpline;
+ Pointer( PPolynomial< Degree > ) baseFunctions;
+ Pointer( BSplineComponents ) baseBSplines;
+
+ BSplineData( void );
+
+ const static int VV_DOT_FLAG = 1;
+ const static int DV_DOT_FLAG = 2;
+ const static int DD_DOT_FLAG = 4;
+ const static int VALUE_FLAG = 1;
+ const static int D_VALUE_FLAG = 2;
+ template< class Real >
+ struct DotTables
+ {
+ size_t functionCount;
+ Pointer( Real ) vvDotTable;
+ Pointer( Real ) dvDotTable;
+ Pointer( Real ) ddDotTable;
+
+ DotTables( void );
+ ~DotTables( void );
+
+ inline size_t Index( int i1 , int i2 ) const;
+ static inline size_t SymmetricIndex( int i1 , int i2 );
+ static inline int SymmetricIndex( int i1 , int i2 , size_t&
index );
+ };
+ template< class Real >
+ struct ValueTables
+ {
+ size_t functionCount , sampleCount;
+ Pointer( Real ) valueTable;
+ Pointer( Real ) dValueTable;
+
+ ValueTables( void );
+ ~ValueTables( void );
+
+ inline size_t Index( int i1 , int i2 ) const;
+ void setSampleSpan( int idx , int& start , int& end , double
smooth=0 ) const;
+ };
+ void set( int maxDepth , int boundaryType=BSplineElements< Degree
>::NONE );
+ template< class Real >
+ typename BSplineData< Degree >::template DotTables< Real >
getDotTables( int flags , bool useDotRatios=true , bool inset=false ) const;
+ template< class Real >
+ typename BSplineData< Degree >::template ValueTables< Real >
getValueTables( int flags , double valueSmooth=0 , double normalSmooth=0 )
const;
+};
+
+template< int Degree1 , int Degree2 > void SetBSplineElementIntegrals( double
integrals[Degree1+1][Degree2+1] );
+
+#include "BSplineData.inl"
+#endif // BSPLINE_DATA_INCLUDED
Property changes on: brlcad/trunk/src/libbg/spsr/BSplineData.h
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/BSplineData.inl
===================================================================
--- brlcad/trunk/src/libbg/spsr/BSplineData.inl (rev 0)
+++ brlcad/trunk/src/libbg/spsr/BSplineData.inl 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,733 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+/////////////////
+// BSplineData //
+/////////////////
+// Support[i]:
+// Odd: i +/- 0.5 * ( 1 + Degree )
+// i - 0.5 * ( 1 + Degree ) < 0
+// <=> i < 0.5 * ( 1 + Degree )
+// i + 0.5 * ( 1 + Degree ) > 0
+// <=> i > - 0.5 * ( 1 + Degree )
+// i + 0.5 * ( 1 + Degree ) > r
+// <=> i > r - 0.5 * ( 1 + Degree )
+// i - 0.5 * ( 1 + Degree ) < r
+// <=> i < r + 0.5 * ( 1 + Degree )
+// Even: i + 0.5 +/- 0.5 * ( 1 + Degree )
+// i - 0.5 * Degree < 0
+// <=> i < 0.5 * Degree
+// i + 1 + 0.5 * Degree > 0
+// <=> i > -1 - 0.5 * Degree
+// i + 1 + 0.5 * Degree > r
+// <=> i > r - 1 - 0.5 * Degree
+// i - 0.5 * Degree < r
+// <=> i < r + 0.5 * Degree
+template< int Degree > inline bool LeftOverlap( unsigned int depth , int
offset )
+{
+ offset <<= 1;
+ if( Degree & 1 ) return (offset < 1+Degree) && (offset > -1-Degree );
+ else return (offset < Degree) && (offset > -2-Degree );
+}
+template< int Degree > inline bool RightOverlap( unsigned int depth , int
offset )
+{
+ offset <<= 1;
+ int r = 1<<(depth+1);
+ if( Degree & 1 ) return (offset > 2-1-Degree) && (offset < 2+1+Degree );
+ else return (offset > 2-2-Degree) && (offset < 2+ Degree );
+}
+template< int Degree > inline int ReflectLeft( unsigned int depth , int offset
)
+{
+ if( Degree&1 ) return -offset;
+ else return -1-offset;
+}
+template< int Degree > inline int ReflectRight( unsigned int depth , int
offset )
+{
+ int r = 1<<(depth+1);
+ if( Degree&1 ) return r -offset;
+ else return r-1-offset;
+}
+
+template< int Degree >
+BSplineData< Degree >::BSplineData( void )
+{
+ functionCount = sampleCount = 0;
+ SetBSplineElementIntegrals< Degree , Degree >( _vvIntegrals );
+ SetBSplineElementIntegrals< Degree , Degree-1 >( _vdIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree >( _dvIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree-1 >( _ddIntegrals );
+}
+
+template< int Degree >
+double BSplineData< Degree >::Integrator::dot( int depth , int off1 , int off2
, bool d1 , bool d2 , bool childParent ) const
+{
+ if( depth<0 || depth>=int( iTables.size() ) ) return 0.;
+ const typename Integrator::IntegralTables& iTable = iTables[depth];
+ if( childParent )
+ {
+ int c = off1&1;
+ off1 >>= 1 , depth--;
+ int ii , d = off2-off1 , res = (1<<depth);
+ if( depth<0 || off1<0 || off2<0 || off1>=res || off2>=res ||
d<-Degree || d>Degree ) return 0;
+ if ( off1< Degree ) ii = off1;
+ else if( off1>=res-Degree ) ii = 2*Degree + off1 - (res-1);
+ else ii = Degree;
+ if ( d1 && d2 ) return
iTable.dd_cpIntegrals[2*ii+c][d+Degree];
+ else if( d1 ) return
iTable.dv_cpIntegrals[2*ii+c][d+Degree];
+ else if( d2 ) return
iTable.vd_cpIntegrals[2*ii+c][d+Degree];
+ else return
iTable.vv_cpIntegrals[2*ii+c][d+Degree];
+ }
+ else
+ {
+ int ii , d = off2-off1 , res = (1<<depth);
+ if( off1<0 || off2<0 || off1>=res || off2>=res || d<-Degree ||
d>Degree ) return 0;
+ if ( off1< Degree ) ii = off1;
+ else if( off1>=res-Degree ) ii = 2*Degree + off1 - (res-1);
+ else ii = Degree;
+ if ( d1 && d2 ) return iTable.dd_ccIntegrals[ii][d+Degree];
+ else if( d1 ) return iTable.dv_ccIntegrals[ii][d+Degree];
+ else if( d2 ) return iTable.vd_ccIntegrals[ii][d+Degree];
+ else return iTable.vv_ccIntegrals[ii][d+Degree];
+ }
+}
+template< int Degree >
+template< int Radius >
+double BSplineData< Degree >::CenterEvaluator< Radius >::value( int depth ,
int off1 , int off2 , bool d , bool childParent ) const
+{
+ if( depth<0 || depth>=int( vTables.size() ) ) return 0.;
+ if( childParent )
+ {
+ int c = off1&1;
+ off1 >>= 1 , depth--;
+ const typename CenterEvaluator::ValueTables& vTable =
vTables[depth];
+ int ii , dd = off1-off2 , res = (1<<depth);
+ if( depth<0 || off1<0 || off2<0 || off1>=res || off2>=res ||
dd<-Radius || dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*3+2*c];
+ else return vTable.vValues[ii][(dd+Radius)*3+2*c];
+ }
+ else
+ {
+ const typename CenterEvaluator::ValueTables& vTable =
vTables[depth];
+ int ii , dd = off1-off2 , res = (1<<depth);
+ if( off1<0 || off2<0 || off1>=res || off2>=res || dd<-Radius ||
dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*3+1];
+ else return vTable.vValues[ii][(dd+Radius)*3+1];
+ }
+}
+template< int Degree >
+template< int Radius >
+double BSplineData< Degree >::CornerEvaluator< Radius >::value( int depth ,
int off1 , int c1 , int off2 , bool d , bool childParent ) const
+{
+ if( c1<0 || c1>=2 )
+ {
+ fprintf( stderr , "[WARNING] Clamping corner to {0,1}\n" );
+ c1 = std::max< int >( 0 , std::min< int >( c1 , 1 ) );
+ }
+ if( depth<0 || depth>=int( vTables.size() ) ) return 0.;
+ if( childParent )
+ {
+ int c = off1&1;
+ off1 >>= 1 , depth--;
+ const typename CornerEvaluator::ValueTables& vTable =
vTables[depth];
+ int ii , dd = off1-off2 , res = (1<<depth);
+ if( depth<0 || off1<0 || off2<0 || off1>=res || off2>=res ||
dd<-Radius || dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*2+c+c1];
+ else return vTable.vValues[ii][(dd+Radius)*2+c+c1];
+ }
+ else
+ {
+ const typename CornerEvaluator::ValueTables& vTable =
vTables[depth];
+ int ii , dd = off1-off2 , res = (1<<depth);
+ if( off1<0 || off2<0 || off1>=res || off2>=res || dd<-Radius ||
dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*2+2*c1];
+ else return vTable.vValues[ii][(dd+Radius)*2+2*c1];
+ }
+}
+template< int Degree >
+void BSplineData< Degree >::set( int maxDepth , int boundaryType )
+{
+ _boundaryType = boundaryType;
+
+ depth = maxDepth;
+ // [Warning] This assumes that the functions spacing is dual
+ functionCount = BinaryNode::CumulativeCenterCount( depth );
+ sampleCount = BinaryNode::CenterCount( depth ) +
BinaryNode::CornerCount( depth );
+ baseFunctions = NewPointer< PPolynomial< Degree > >( functionCount );
+ baseBSplines = NewPointer< BSplineComponents >( functionCount );
+
+ baseFunction = PPolynomial< Degree >::BSpline();
+ for( int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree
>::BSplineComponent( i ).shift( double(-(Degree+1)/2) + i - 0.5 );
+ dBaseFunction = baseFunction.derivative();
+ StartingPolynomial< Degree > sPolys[Degree+4];
+
+ for( int i=0 ; i<Degree+3 ; i++ )
+ {
+ sPolys[i].start = double(-(Degree+1)/2) + i - 1.5;
+ sPolys[i].p *= 0;
+ if( i<=Degree ) sPolys[i].p += baseBSpline[i
].shift( -1 ) * _boundaryType;
+ if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1];
+ for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
+ }
+ leftBaseFunction.set( sPolys , Degree+3 );
+ for( int i=0 ; i<Degree+3 ; i++ )
+ {
+ sPolys[i].start = double(-(Degree+1)/2) + i - 0.5;
+ sPolys[i].p *= 0;
+ if( i<=Degree ) sPolys[i].p += baseBSpline[i ];
+ if( i>=1 && i<=Degree+1 ) sPolys[i].p +=
baseBSpline[i-1].shift( 1 ) * _boundaryType;
+ for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
+ }
+ rightBaseFunction.set( sPolys , Degree+3 );
+ for( int i=0 ; i<Degree+4 ; i++ )
+ {
+ sPolys[i].start = double(-(Degree+1)/2) + i - 1.5;
+ sPolys[i].p *= 0;
+ if( i<=Degree ) sPolys[i].p += baseBSpline[i
].shift( -1 ) * _boundaryType; // The left-shifted B-spline
+ if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1];
// The centered B-Spline
+ if( i>=2 && i<=Degree+2 ) sPolys[i].p +=
baseBSpline[i-2].shift( 1 ) * _boundaryType; // The right-shifted B-spline
+ for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
+ }
+ leftRightBaseFunction.set( sPolys , Degree+4 );
+
+ dLeftBaseFunction = leftBaseFunction.derivative();
+ dRightBaseFunction = rightBaseFunction.derivative();
+ dLeftRightBaseFunction = leftRightBaseFunction.derivative();
+ leftRightBSpline = leftBSpline = rightBSpline = baseBSpline;
+ leftBSpline [1] += leftBSpline[2].shift( -1 ) , leftBSpline[0] *= 0;
+ rightBSpline[1] += rightBSpline[0].shift( 1 ) , rightBSpline[2] *= 0;
+ leftRightBSpline[1] += leftRightBSpline[2].shift( -1 ) +
leftRightBSpline[0].shift( 1 ) , leftRightBSpline[0] *= 0 , leftRightBSpline[2]
*= 0 ;
+
+ double c , w;
+ for( size_t i=0 ; i<functionCount ; i++ )
+ {
+ BinaryNode::CenterAndWidth( int(i) , c , w );
+ baseFunctions[i] = baseFunction.scale(w).shift(c);
+ baseBSplines[i] = baseBSpline.scale(w).shift(c);
+ if( _boundaryType )
+ {
+ int d , off , r;
+ BinaryNode::DepthAndOffset( int(i) , d , off );
+ r = 1<<d;
+ if ( off==0 && off==r-1 ) baseFunctions[i] =
leftRightBaseFunction.scale(w).shift(c);
+ else if( off==0 ) baseFunctions[i] =
leftBaseFunction.scale(w).shift(c);
+ else if( off==r-1 ) baseFunctions[i] =
rightBaseFunction.scale(w).shift(c);
+ if ( off==0 && off==r-1 ) baseBSplines [i] =
leftRightBSpline.scale(w).shift(c);
+ else if( off==0 ) baseBSplines [i] =
leftBSpline.scale(w).shift(c);
+ else if( off==r-1 ) baseBSplines [i] =
rightBSpline.scale(w).shift(c);
+ }
+ }
+}
+template< int Degree >
+double BSplineData< Degree >::dot( int depth1 , int off1 , int depth2 , int
off2 , bool d1 , bool d2 , bool inset ) const
+{
+ const int _Degree1 = (d1 ? (Degree-1) : Degree) , _Degree2 = (d2 ?
(Degree-1) : Degree);
+ int sums[ Degree+1 ][ Degree+1 ];
+
+ int depth = std::max< int >( depth1 , depth2 );
+
+ BSplineElements< Degree > b1( 1<<depth1 , off1 , _boundaryType , inset
? ( 1<<(depth1-2) ) : 0 ) , b2( 1<<depth2 , off2 , _boundaryType , inset ? (
1<<(depth2-2) ) : 0 );
+
+ BSplineElements< Degree > b;
+ while( depth1<depth ) b=b1 , b.upSample( b1 ) , depth1++;
+ while( depth2<depth ) b=b2 , b.upSample( b2 ) , depth2++;
+
+ BSplineElements< Degree-1 > db1 , db2;
+ b1.differentiate( db1 ) , b2.differentiate( db2 );
+
+ int start1=-1 , end1=-1 , start2=-1 , end2=-1;
+ for( int i=0 ; i<int( b1.size() ) ; i++ ) for( int j=0 ; j<=Degree ;
j++ )
+ {
+ if( b1[i][j] && start1==-1 ) start1 = i;
+ if( b1[i][j] ) end1 = i+1;
+ if( b2[i][j] && start2==-1 ) start2 = i;
+ if( b2[i][j] ) end2 = i+1;
+ }
+ if( start1==end1 || start2==end2 || start1>=end2 || start2>=end1 )
return 0.;
+ int start = std::max< int >( start1 , start2 ) , end = std::min< int >(
end1 , end2 );
+ memset( sums , 0 , sizeof( sums ) );
+ for( int i=start ; i<end ; i++ ) for( int j=0 ; j<=_Degree1 ; j++ )
for( int k=0 ; k<=_Degree2 ; k++ ) sums[j][k] += ( d1 ? db1[i][j] : b1[i][j] )
* ( d2 ? db2[i][k] : b2[i][k] );
+ double _dot = 0;
+ if ( d1 && d2 ) for( int j=0 ; j<=_Degree1 ; j++ ) for( int k=0 ;
k<=_Degree2 ; k++ ) _dot += _ddIntegrals[j][k] * sums[j][k];
+ else if( d1 ) for( int j=0 ; j<=_Degree1 ; j++ ) for( int k=0 ;
k<=_Degree2 ; k++ ) _dot += _dvIntegrals[j][k] * sums[j][k];
+ else if( d2 ) for( int j=0 ; j<=_Degree1 ; j++ ) for( int k=0 ;
k<=_Degree2 ; k++ ) _dot += _vdIntegrals[j][k] * sums[j][k];
+ else for( int j=0 ; j<=_Degree1 ; j++ ) for( int k=0 ;
k<=_Degree2 ; k++ ) _dot += _vvIntegrals[j][k] * sums[j][k];
+ _dot /= b1.denominator;
+ _dot /= b2.denominator;
+ if ( d1 && d2 ) return _dot * (1<<depth);
+ else if( d1 || d2 ) return _dot;
+ else return _dot / (1<<depth);
+}
+template< int Degree >
+double BSplineData< Degree >::value( int depth , int off , double
smoothingRadius , double s , bool d , bool inset ) const
+{
+ PPolynomial< Degree+1 > function;
+ PPolynomial< Degree > dFunction;
+
+ if( off<0 || off>=(1<<depth) ) return 0;
+ int idx = BinaryNode::CenterIndex( depth , off );
+
+ if( smoothingRadius>0 ) function = baseFunctions[idx].MovingAverage(
smoothingRadius );
+ else function = baseFunctions[idx];
+ dFunction = function.derivative();
+
+ if( d ) return dFunction(s);
+ else return function(s);
+}
+template< int Degree >
+void BSplineData< Degree >::setIntegrator( Integrator& integrator , bool inset
, bool useDotRatios ) const
+{
+ integrator.iTables.resize( depth+1 );
+ for( int d=0 ; d<=depth ; d++ ) for( int i=0 ; i<=2*Degree ; i++ ) for(
int j=-Degree ; j<=Degree ; j++ )
+ {
+ int res = 1<<d , ii = (i<=Degree ? i : i+res-1 - 2*Degree );
+ integrator.iTables[d].vv_ccIntegrals[i][j+Degree] = dot( d , ii
, d , ii+j , false , false , inset );
+ integrator.iTables[d].dv_ccIntegrals[i][j+Degree] = dot( d , ii
, d , ii+j , true , false , inset );
+ integrator.iTables[d].vd_ccIntegrals[i][j+Degree] = dot( d , ii
, d , ii+j , false , true , inset );
+ integrator.iTables[d].dd_ccIntegrals[i][j+Degree] = dot( d , ii
, d , ii+j , true , true , inset );
+ if( useDotRatios )
+ {
+ integrator.iTables[d].dv_ccIntegrals[i][j+Degree] /=
integrator.iTables[d].vv_ccIntegrals[i][j+Degree];
+ integrator.iTables[d].vd_ccIntegrals[i][j+Degree] /=
integrator.iTables[d].vv_ccIntegrals[i][j+Degree];
+ integrator.iTables[d].dd_ccIntegrals[i][j+Degree] /=
integrator.iTables[d].vv_ccIntegrals[i][j+Degree];
+ }
+ }
+ for( int d=1 ; d<=depth ; d++ ) for( int i=0 ; i<=2*Degree ; i++ ) for(
int j=-Degree ; j<=Degree ; j++ )
+ {
+ int res = 1<<d , ii = (i<=Degree ? i : i+(res/2)-1 - 2*Degree );
+ for( int c=0 ; c<2 ; c++ )
+ {
+ integrator.iTables[d].vv_cpIntegrals[2*i+c][j+Degree] =
dot( d , 2*ii+c , d-1 , ii+j , false , false , inset );
+ integrator.iTables[d].dv_cpIntegrals[2*i+c][j+Degree] =
dot( d , 2*ii+c , d-1 , ii+j , true , false , inset );
+ integrator.iTables[d].vd_cpIntegrals[2*i+c][j+Degree] =
dot( d , 2*ii+c , d-1 , ii+j , false , true , inset );
+ integrator.iTables[d].dd_cpIntegrals[2*i+c][j+Degree] =
dot( d , 2*ii+c , d-1 , ii+j , true , true , inset );
+ if( useDotRatios )
+ {
+
integrator.iTables[d].dv_cpIntegrals[2*i+c][j+Degree] /=
integrator.iTables[d].vv_cpIntegrals[2*i+c][j+Degree];
+
integrator.iTables[d].vd_cpIntegrals[2*i+c][j+Degree] /=
integrator.iTables[d].vv_cpIntegrals[2*i+c][j+Degree];
+
integrator.iTables[d].dd_cpIntegrals[2*i+c][j+Degree] /=
integrator.iTables[d].vv_cpIntegrals[2*i+c][j+Degree];
+ }
+ }
+ }
+}
+template< int Degree >
+template< int Radius >
+void BSplineData< Degree >::setCenterEvaluator( CenterEvaluator< Radius >&
evaluator , double smoothingRadius , double dSmoothingRadius , bool inset )
const
+{
+ evaluator.vTables.resize( depth+1 );
+ for( int d=0 ; d<=depth ; d++ ) for( int i=0 ; i<=2*Degree ; i++ ) for(
int j=-Radius ; j<=Radius ; j++ ) for( int k=-1 ; k<=1 ; k++ )
+ {
+ int res = 1<<d , ii = (i<=Degree ? i : i+res-1 - 2*Degree );
+ double s = 0.5+ii+j+0.25*k;
+ evaluator.vTables[d].vValues[i][(j+Radius)*3+(k+1)] = value( d
, ii , smoothingRadius , s/res , false , inset );
+ evaluator.vTables[d].dValues[i][(j+Radius)*3+(k+1)] = value( d
, ii , dSmoothingRadius , s/res , true , inset );
+ }
+}
+template< int Degree >
+template< int Radius >
+void BSplineData< Degree >::setCornerEvaluator( CornerEvaluator< Radius >&
evaluator , double smoothingRadius , double dSmoothingRadius , bool inset )
const
+{
+ evaluator.vTables.resize( depth+1 );
+ for( int d=0 ; d<=depth ; d++ ) for( int i=0 ; i<=2*Degree ; i++ ) for(
int j=-Radius ; j<=Radius ; j++ ) for( int k=0 ; k<=2 ; k++ )
+ {
+ int res = 1<<d , ii = (i<=Degree ? i : i+res-1 - 2*Degree );
+ double s = ii+j+0.5*k;
+ evaluator.vTables[d].vValues[i][(j+Radius)*2+k] = value( d , ii
, smoothingRadius , s/res , false , inset );
+ evaluator.vTables[d].dValues[i][(j+Radius)*2+k] = value( d , ii
, dSmoothingRadius , s/res , true , inset );
+ }
+}
+
+
+template< int Degree >
+template< class Real >
+BSplineData< Degree >::DotTables< Real >::DotTables( void )
+{
+ vvDotTable = NullPointer< Real >();
+ dvDotTable = NullPointer< Real >();
+ ddDotTable = NullPointer< Real >();
+}
+template< int Degree >
+template< class Real >
+BSplineData< Degree >::DotTables< Real >::~DotTables( void )
+{
+ DeletePointer( vvDotTable );
+ DeletePointer( dvDotTable );
+ DeletePointer( ddDotTable );
+}
+template< int Degree >
+template< class Real >
+inline size_t BSplineData< Degree >::DotTables< Real >::Index( int i1 , int i2
) const { return size_t(i1)*functionCount + size_t(i2); }
+template< int Degree >
+template< class Real >
+inline size_t BSplineData< Degree >::DotTables< Real >::SymmetricIndex( int i1
, int i2 )
+{
+ size_t _i1 = i1 , _i2 = i2;
+ if( i1>i2 ) return ((_i1*_i1+i1)>>1)+_i2;
+ else return ((_i2*_i2+i2)>>1)+_i1;
+}
+template< int Degree >
+template< class Real >
+inline int BSplineData< Degree >::DotTables< Real >::SymmetricIndex( int i1 ,
int i2 , size_t& index )
+{
+ size_t _i1 = i1 , _i2 = i2;
+ if( i1<i2 )
+ {
+ index = ((_i2*_i2+_i2)>>1)+_i1;
+ return 1;
+ }
+ else
+ {
+ index = ((_i1*_i1+_i1)>>1)+_i2;
+ return 0;
+ }
+}
+template< int Degree >
+template< class Real >
+typename BSplineData< Degree >::template DotTables< Real > BSplineData< Degree
>::getDotTables( int flags , bool useDotRatios , bool inset ) const
+{
+ typename BSplineData< Degree >::template DotTables< Real > dTables;
+ dTables.functionCount = functionCount;
+
+ size_t size = ( functionCount*functionCount + functionCount )>>1;
+ size_t fullSize = functionCount*functionCount;
+ if( flags & VV_DOT_FLAG )
+ {
+ dTables.vvDotTable = NewPointer< Real >( size );
+ memset( dTables.vvDotTable , 0 , sizeof(Real)*size );
+ }
+ if( flags & DV_DOT_FLAG )
+ {
+ dTables.dvDotTable = NewPointer< Real >( fullSize );
+ memset( dTables.dvDotTable , 0 , sizeof(Real)*fullSize );
+ }
+ if( flags & DD_DOT_FLAG )
+ {
+ dTables.ddDotTable = NewPointer< Real >( size );
+ memset( dTables.ddDotTable , 0 , sizeof(Real)*size );
+ }
+ int vvSums[Degree+1][Degree+1];
+ int vdSums[Degree+1][Degree ];
+ int dvSums[Degree ][Degree+1];
+ int ddSums[Degree ][Degree ];
+ double vvIntegrals[Degree+1][Degree+1];
+ double vdIntegrals[Degree+1][Degree ];
+ double dvIntegrals[Degree ][Degree+1];
+ double ddIntegrals[Degree ][Degree ];
+ SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
+ SetBSplineElementIntegrals< Degree , Degree-1 >( vdIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree >( dvIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree-1 >( ddIntegrals );
+
+ for( int d1=0 ; d1<=depth ; d1++ ) for( int off1=0 ; off1<(1<<d1) ;
off1++ )
+ {
+ int ii = BinaryNode::CenterIndex( d1 , off1 );
+ BSplineElements< Degree > b1( 1<<d1 , off1 , _boundaryType ,
inset ? ( 1<<(d1-2) ) : 0 );
+ BSplineElements< Degree-1 > db1;
+ b1.differentiate( db1 );
+ int start1 , end1;
+
+ start1 = -1 , end1 = -1;
+ for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ;
j<=Degree ; j++ )
+ {
+ if( b1[i][j] && start1==-1 ) start1 = i;
+ if( b1[i][j] ) end1 = i+1;
+ }
+ if( start1==end1 ) continue;
+ for( int d2=d1 ; d2<=depth ; d2++ )
+ {
+ for( int off2=0 ; off2<(1<<d2) ; off2++ )
+ {
+ int start2 = off2-Degree;
+ int end2 = off2+Degree+1;
+ if( start2>=end1 || start1>=end2 ) continue;
+ start2 = std::max< int >( start1 , start2 );
+ end2 = std::min< int >( end1 , end2 );
+ if( d1==d2 && off2<off1 ) continue;
+ int jj = BinaryNode::CenterIndex( d2 , off2 );
+ BSplineElements< Degree > b2( 1<<d2 , off2 ,
_boundaryType , inset ? ( 1<<(d2-2) ) : 0 );
+ BSplineElements< Degree-1 > db2;
+ b2.differentiate( db2 );
+
+ size_t idx = DotTables< Real >::SymmetricIndex(
ii , jj );
+ size_t idx1 = DotTables< Real >::Index( ii , jj
) , idx2 = DotTables< Real >::Index( jj , ii );
+
+ memset( vvSums , 0 , sizeof( int ) * ( Degree+1
) * ( Degree+1 ) );
+ memset( vdSums , 0 , sizeof( int ) * ( Degree+1
) * ( Degree ) );
+ memset( dvSums , 0 , sizeof( int ) * ( Degree
) * ( Degree+1 ) );
+ memset( ddSums , 0 , sizeof( int ) * ( Degree
) * ( Degree ) );
+ for( int i=start2 ; i<end2 ; i++ )
+ {
+ for( int j=0 ; j<=Degree ; j++ ) for(
int k=0 ; k<=Degree ; k++ ) vvSums[j][k] += b1[i][j] * b2[i][k];
+ for( int j=0 ; j<=Degree ; j++ ) for(
int k=0 ; k< Degree ; k++ ) vdSums[j][k] += b1[i][j] * db2[i][k];
+ for( int j=0 ; j< Degree ; j++ ) for(
int k=0 ; k<=Degree ; k++ ) dvSums[j][k] += db1[i][j] * b2[i][k];
+ for( int j=0 ; j< Degree ; j++ ) for(
int k=0 ; k< Degree ; k++ ) ddSums[j][k] += db1[i][j] * db2[i][k];
+ }
+ double vvDot = 0 , dvDot = 0 , vdDot = 0 ,
ddDot = 0;
+ for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ;
k<=Degree ; k++ ) vvDot += vvIntegrals[j][k] * vvSums[j][k];
+ for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ;
k< Degree ; k++ ) vdDot += vdIntegrals[j][k] * vdSums[j][k];
+ for( int j=0 ; j< Degree ; j++ ) for( int k=0 ;
k<=Degree ; k++ ) dvDot += dvIntegrals[j][k] * dvSums[j][k];
+ for( int j=0 ; j< Degree ; j++ ) for( int k=0 ;
k< Degree ; k++ ) ddDot += ddIntegrals[j][k] * ddSums[j][k];
+ vvDot /= (1<<d2);
+ ddDot *= (1<<d2);
+ vvDot /= ( b1.denominator * b2.denominator );
+ dvDot /= ( b1.denominator * b2.denominator );
+ vdDot /= ( b1.denominator * b2.denominator );
+ ddDot /= ( b1.denominator * b2.denominator );
+ if( fabs(vvDot)<1e-15 ) continue;
+ if( flags & VV_DOT_FLAG )
dTables.vvDotTable[idx] = Real( vvDot );
+ if( flags & DV_DOT_FLAG )
dTables.dvDotTable[idx1] = Real( dvDot );
+ if( flags & DV_DOT_FLAG )
dTables.dvDotTable[idx2] = Real( dvDot );
+ if( flags & DD_DOT_FLAG )
dTables.ddDotTable[idx ] = Real( ddDot );
+ }
+ BSplineElements< Degree > b;
+ b = b1;
+ b.upSample( b1 );
+ b1.differentiate( db1 );
+ start1 = -1;
+ for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ;
j<=Degree ; j++ )
+ {
+ if( b1[i][j] && start1==-1 ) start1 = i;
+ if( b1[i][j] ) end1 = i+1;
+ }
+ }
+ }
+ return dTables;
+}
+template< int Degree >
+template< class Real >
+BSplineData< Degree >::ValueTables< Real >::ValueTables( void )
+{
+ valueTable = NullPointer< Real >();
+ dValueTable = NullPointer< Real >();
+}
+template< int Degree >
+template< class Real >
+BSplineData< Degree >::ValueTables< Real >::~ValueTables( void )
+{
+ DeletePointer( valueTable );
+ DeletePointer( dValueTable );
+}
+template< int Degree >
+template< class Real >
+inline size_t BSplineData< Degree >::ValueTables< Real >::Index( int i1 , int
i2 ) const { return size_t(i1)*functionCount + size_t(i2); }
+template< int Degree >
+template< class Real >
+typename BSplineData< Degree >::template ValueTables< Real > BSplineData<
Degree >::getValueTables( int flags , double valueSmooth , double
derivativeSmooth ) const
+{
+ typename BSplineData< Degree >::template ValueTables< Real > vTables;
+ vTables.functionCount = functionCount;
+ vTables.sampleCount = sampleCount;
+
+ if( flags & VALUE_FLAG ) vTables.valueTable = NewPointer< Real >(
functionCount*sampleCount );
+ if( flags & D_VALUE_FLAG ) vTables.dValueTable = NewPointer< Real >(
functionCount*sampleCount );
+ PPolynomial< Degree+1 > function;
+ PPolynomial< Degree > dFunction;
+ for( size_t i=0 ; i<functionCount ; i++ )
+ {
+ if( valueSmooth>0 )
function=baseFunctions[i].MovingAverage( valueSmooth );
+ else function=baseFunctions[i];
+ if( derivativeSmooth>0 )
dFunction=baseFunctions[i].derivative().MovingAverage( derivativeSmooth );
+ else
dFunction=baseFunctions[i].derivative();
+
+ for( size_t j=0 ; j<sampleCount ; j++ )
+ {
+ double x=double(j)/(sampleCount-1);
+ if( flags & VALUE_FLAG )
vTables.valueTable[j*functionCount+i] = Real( function(x));
+ if( flags & D_VALUE_FLAG )
vTables.dValueTable[j*functionCount+i] = Real(dFunction(x));
+ }
+ }
+ return vTables;
+}
+template< int Degree >
+template< class Real >
+void BSplineData< Degree >::ValueTables< Real >::setSampleSpan( int idx , int&
start , int& end , double smooth ) const
+{
+ int d , off , res;
+ BinaryNode::DepthAndOffset( idx , d , off );
+ res = 1<<d;
+ double _start = ( off + 0.5 - 0.5*(Degree+1) ) / res - smooth;
+ double _end = ( off + 0.5 + 0.5*(Degree+1) ) / res + smooth;
+ // (start)/(sampleCount-1) >_start &&
(start-1)/(sampleCount-1)<=_start
+ // => start > _start * (sampleCount-1 ) && start <=
_start*(sampleCount-1) + 1
+ // => _start * (sampleCount-1) + 1 >= start > _start * (sampleCount-1)
+ start = int( floor( _start * (sampleCount-1) + 1 ) );
+ if( start<0 ) start = 0;
+ // (end)/(sampleCount-1)<_end && (end+1)/(sampleCount-1)>=_end
+ // => end < _end * (sampleCount-1 ) && end >= _end*(sampleCount-1) - 1
+ // => _end * (sampleCount-1) > end >= _end * (sampleCount-1) - 1
+ end = int( ceil( _end * (sampleCount-1) - 1 ) );
+ if( end>=int(sampleCount) ) end = int(sampleCount)-1;
+}
+
+
+/////////////////////
+// BSplineElements //
+/////////////////////
+template< int Degree >
+BSplineElements< Degree >::BSplineElements( int res , int offset , int
boundary , int inset )
+{
+ denominator = 1;
+ std::vector< BSplineElementCoefficients< Degree > >::resize( res ,
BSplineElementCoefficients< Degree >() );
+
+ for( int i=0 ; i<=Degree ; i++ )
+ {
+ int idx = -_off + offset + i;
+ if( idx>=0 && idx<res ) (*this)[idx][i] = 1;
+ }
+ if( boundary!=0 )
+ {
+ _addLeft( offset-2*res , boundary ) , _addRight( offset+2*res ,
boundary );
+ if( Degree&1 ) _addLeft( offset-res , boundary ) , _addRight(
offset+res , boundary );
+ else _addLeft( -offset-1 , boundary ) , _addRight(
-offset-1+2*res , boundary );
+ }
+ if( inset ) for( int i=0 ; i<inset && i<res ; i++ ) for( int j=0 ;
j<=Degree ; j++ ) (*this)[i][j] = (*this)[res-1-i][j] = 0;
+}
+template< int Degree >
+void BSplineElements< Degree >::_addLeft( int offset , int boundary )
+{
+ int res = int( std::vector< BSplineElementCoefficients< Degree >
>::size() );
+ bool set = false;
+ for( int i=0 ; i<=Degree ; i++ )
+ {
+ int idx = -_off + offset + i;
+ if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set =
true;
+ }
+ if( set ) _addLeft( offset-2*res , boundary );
+}
+template< int Degree >
+void BSplineElements< Degree >::_addRight( int offset , int boundary )
+{
+ int res = int( std::vector< BSplineElementCoefficients< Degree >
>::size() );
+ bool set = false;
+ for( int i=0 ; i<=Degree ; i++ )
+ {
+ int idx = -_off + offset + i;
+ if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set =
true;
+ }
+ if( set ) _addRight( offset+2*res , boundary );
+}
+template< int Degree >
+void BSplineElements< Degree >::upSample( BSplineElements< Degree >& high )
const
+{
+ fprintf( stderr , "[ERROR] B-spline up-sampling not supported for
degree %d\n" , Degree );
+ exit( 0 );
+}
+template<>
+void BSplineElements< 1 >::upSample( BSplineElements< 1 >& high ) const
+{
+ high.resize( size()*2 );
+ high.assign( high.size() , BSplineElementCoefficients<1>() );
+ for( int i=0 ; i<int(size()) ; i++ )
+ {
+ high[2*i+0][0] += 1 * (*this)[i][0];
+ high[2*i+0][1] += 0 * (*this)[i][0];
+ high[2*i+1][0] += 2 * (*this)[i][0];
+ high[2*i+1][1] += 1 * (*this)[i][0];
+
+ high[2*i+0][0] += 1 * (*this)[i][1];
+ high[2*i+0][1] += 2 * (*this)[i][1];
+ high[2*i+1][0] += 0 * (*this)[i][1];
+ high[2*i+1][1] += 1 * (*this)[i][1];
+ }
+ high.denominator = denominator * 2;
+}
+template<>
+void BSplineElements< 2 >::upSample( BSplineElements< 2 >& high ) const
+{
+ /* /----\
*/
+ /* / \
*/
+ /* / \ = 1 /--\ +3 /--\ +3 /--\ +1
/--\ */
+ /* / \ / \ / \ / \
/ \ */
+ /* |----------| |----------| |----------| |----------|
|----------| */
+
+ high.resize( size()*2 );
+ high.assign( high.size() , BSplineElementCoefficients<2>() );
+ for( int i=0 ; i<int(size()) ; i++ )
+ {
+ high[2*i+0][0] += 1 * (*this)[i][0];
+ high[2*i+0][1] += 0 * (*this)[i][0];
+ high[2*i+0][2] += 0 * (*this)[i][0];
+ high[2*i+1][0] += 3 * (*this)[i][0];
+ high[2*i+1][1] += 1 * (*this)[i][0];
+ high[2*i+1][2] += 0 * (*this)[i][0];
+
+ high[2*i+0][0] += 3 * (*this)[i][1];
+ high[2*i+0][1] += 3 * (*this)[i][1];
+ high[2*i+0][2] += 1 * (*this)[i][1];
+ high[2*i+1][0] += 1 * (*this)[i][1];
+ high[2*i+1][1] += 3 * (*this)[i][1];
+ high[2*i+1][2] += 3 * (*this)[i][1];
+
+ high[2*i+0][0] += 0 * (*this)[i][2];
+ high[2*i+0][1] += 1 * (*this)[i][2];
+ high[2*i+0][2] += 3 * (*this)[i][2];
+ high[2*i+1][0] += 0 * (*this)[i][2];
+ high[2*i+1][1] += 0 * (*this)[i][2];
+ high[2*i+1][2] += 1 * (*this)[i][2];
+ }
+ high.denominator = denominator * 4;
+}
+
+template< int Degree >
+void BSplineElements< Degree >::differentiate( BSplineElements< Degree-1 >& d
) const
+{
+ d.resize( std::vector< BSplineElementCoefficients< Degree > >::size() );
+ d.assign( d.size() , BSplineElementCoefficients< Degree-1 >() );
+ for( int i=0 ; i<int(std::vector< BSplineElementCoefficients< Degree >
>::size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
+ {
+ if( j-1>=0 ) d[i][j-1] -= (*this)[i][j];
+ if( j<Degree ) d[i][j ] += (*this)[i][j];
+ }
+ d.denominator = denominator;
+}
+// If we were really good, we would implement this integral table to store
+// rational values to improve precision...
+template< int Degree1 , int Degree2 >
+void SetBSplineElementIntegrals( double integrals[Degree1+1][Degree2+1] )
+{
+ for( int i=0 ; i<=Degree1 ; i++ )
+ {
+ Polynomial< Degree1 > p1 = Polynomial< Degree1
>::BSplineComponent( i );
+ for( int j=0 ; j<=Degree2 ; j++ )
+ {
+ Polynomial< Degree2 > p2 = Polynomial< Degree2
>::BSplineComponent( j );
+ integrals[i][j] = ( p1 * p2 ).integral( 0 , 1 );
+ }
+ }
+}
Property changes on: brlcad/trunk/src/libbg/spsr/BSplineData.inl
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/BinaryNode.h
===================================================================
--- brlcad/trunk/src/libbg/spsr/BinaryNode.h (rev 0)
+++ brlcad/trunk/src/libbg/spsr/BinaryNode.h 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,78 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+#ifndef BINARY_NODE_INCLUDED
+#define BINARY_NODE_INCLUDED
+
+#define MSVC_2010_FIX 1
+
+
+class BinaryNode
+{
+public:
+ static inline int CenterCount( int depth ) { return 1<<depth; }
+ static inline int CornerCount( int depth ) { return (1<<depth)+1; }
+ static inline int CumulativeCenterCount( int maxDepth ) { return
(1<<(maxDepth+1))-1; }
+ static inline int CumulativeCornerCount( int maxDepth ) { return
(1<<(maxDepth+1))+maxDepth; }
+ static inline int CenterIndex( int depth , int offSet ) { return
(1<<depth)+offSet-1; }
+ static inline int CornerIndex( int depth , int offSet ) { return
(1<<depth)+offSet+depth; }
+
+ static inline int CornerIndex( int maxDepth , int depth , int offSet ,
int forwardCorner ){ return (offSet+forwardCorner)<<(maxDepth-depth); }
+ template< class Real > static inline Real CornerIndexPosition(int
index,int maxDepth){ return Real(index)/(1<<maxDepth); }
+ template< class Real > static inline Real Width(int depth){ return
Real(1.0/(1<<depth)); }
+ template< class Real > static inline void CenterAndWidth( int depth ,
int offset , Real& center , Real& width )
+ {
+ width=Real (1.0/(1<<depth) );
+ center=Real((0.5+offset)*width);
+ }
+ template< class Real > static inline void CenterAndWidth( int idx ,
Real& center , Real& width )
+ {
+ int depth , offset;
+ DepthAndOffset( idx , depth , offset );
+ CenterAndWidth( depth , offset , center , width );
+ }
+ static inline void DepthAndOffset( int idx , int& depth , int& offset )
+ {
+ int i=idx+1;
+#if MSVC_2010_FIX
+ depth = 0;
+#else // !MSVC_2010_FIX
+ depth = -1;
+#endif // MSVC_2010_FIX
+ while( i )
+ {
+ i >>= 1;
+ depth++;
+ }
+#if MSVC_2010_FIX
+ depth--;
+#endif // MSVC_2010_FIX
+ offset = ( idx+1 ) - (1<<depth);
+ }
+};
+#endif // BINARY_NODE_INCLUDED
Property changes on: brlcad/trunk/src/libbg/spsr/BinaryNode.h
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/Factor.cpp
===================================================================
--- brlcad/trunk/src/libbg/spsr/Factor.cpp (rev 0)
+++ brlcad/trunk/src/libbg/spsr/Factor.cpp 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,285 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic push
+#endif
+#if defined(__clang__)
+# pragma clang diagnostic push
+#endif
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic ignored "-Wfloat-equal"
+#endif
+#if defined(__clang__)
+# pragma clang diagnostic ignored "-Wfloat-equal"
+#endif
+
+//////////////////////
+// Polynomial Roots //
+//////////////////////
+#include <math.h>
+#include "Factor.h"
+int Factor(double a1,double a0,double roots[1][2],double EPS){
+ if(fabs(a1)<=EPS){return 0;}
+ roots[0][0]=-a0/a1;
+ roots[0][1]=0;
+ return 1;
+}
+int Factor(double a2,double a1,double a0,double roots[2][2],double EPS){
+ double d;
+ if(fabs(a2)<=EPS){return Factor(a1,a0,roots,EPS);}
+
+ d=a1*a1-4*a0*a2;
+ a1/=(2*a2);
+ if(d<0){
+ d=sqrt(-d)/(2*a2);
+ roots[0][0]=roots[1][0]=-a1;
+ roots[0][1]=-d;
+ roots[1][1]= d;
+ }
+ else{
+ d=sqrt(d)/(2*a2);
+ roots[0][1]=roots[1][1]=0;
+ roots[0][0]=-a1-d;
+ roots[1][0]=-a1+d;
+ }
+ return 2;
+}
+// Solution taken from: http://mathworld.wolfram.com/CubicFormula.html
+// and http://www.csit.fsu.edu/~burkardt/f_src/subpak/subpak.f90
+int Factor(double a3,double a2,double a1,double a0,double roots[3][2],double
EPS){
+ double q,r,r2,q3;
+
+ if(fabs(a3)<=EPS){return Factor(a2,a1,a0,roots,EPS);}
+ a2/=a3;
+ a1/=a3;
+ a0/=a3;
+
+ q=-(3*a1-a2*a2)/9;
+ r=-(9*a2*a1-27*a0-2*a2*a2*a2)/54;
+ r2=r*r;
+ q3=q*q*q;
+
+ if(r2<q3){
+ double sqrQ=sqrt(q);
+ double theta = acos ( r / (sqrQ*q) );
+ double cTheta=cos(theta/3)*sqrQ;
+ double sTheta=sin(theta/3)*sqrQ*SQRT_3/2;
+ roots[0][1]=roots[1][1]=roots[2][1]=0;
+ roots[0][0]=-2*cTheta;
+ roots[1][0]=-2*(-cTheta*0.5-sTheta);
+ roots[2][0]=-2*(-cTheta*0.5+sTheta);
+ }
+ else{
+ double s1,s2,sqr=sqrt(r2-q3);
+ double t;
+ t=-r+sqr;
+ if(t<0){s1=-pow(-t,1.0/3);}
+ else{s1=pow(t,1.0/3);}
+ t=-r-sqr;
+ if(t<0){s2=-pow(-t,1.0/3);}
+ else{s2=pow(t,1.0/3);}
+ roots[0][1]=0;
+ roots[0][0]=s1+s2;
+ s1/=2;
+ s2/=2;
+ roots[1][0]= roots[2][0]=-s1-s2;
+ roots[1][1]= SQRT_3*(s1-s2);
+ roots[2][1]=-roots[1][1];
+ }
+ roots[0][0]-=a2/3;
+ roots[1][0]-=a2/3;
+ roots[2][0]-=a2/3;
+ return 3;
+}
+double ArcTan2(double y,double x){
+ /* This first case should never happen */
+ if(y==0 && x==0){return 0;}
+ if(x==0){
+ if(y>0){return PI/2.0;}
+ else{return -PI/2.0;}
+ }
+ if(x>=0){return atan(y/x);}
+ else{
+ if(y>=0){return atan(y/x)+PI;}
+ else{return atan(y/x)-PI;}
+ }
+}
+double Angle(const double in[2]){
+ if((in[0]*in[0]+in[1]*in[1])==0.0){return 0;}
+ else{return ArcTan2(in[1],in[0]);}
+}
+void Sqrt(const double in[2],double out[2]){
+ double r=sqrt(sqrt(in[0]*in[0]+in[1]*in[1]));
+ double a=Angle(in)*0.5;
+ out[0]=r*cos(a);
+ out[1]=r*sin(a);
+}
+void Add(const double in1[2],const double in2[2],double out[2]){
+ out[0]=in1[0]+in2[0];
+ out[1]=in1[1]+in2[1];
+}
+void Subtract(const double in1[2],const double in2[2],double out[2]){
+ out[0]=in1[0]-in2[0];
+ out[1]=in1[1]-in2[1];
+}
+void Multiply(const double in1[2],const double in2[2],double out[2]){
+ out[0]=in1[0]*in2[0]-in1[1]*in2[1];
+ out[1]=in1[0]*in2[1]+in1[1]*in2[0];
+}
+void Divide(const double in1[2],const double in2[2],double out[2]){
+ double temp[2];
+ double l=in2[0]*in2[0]+in2[1]*in2[1];
+ temp[0]= in2[0]/l;
+ temp[1]=-in2[1]/l;
+ Multiply(in1,temp,out);
+}
+// Solution taken from: http://mathworld.wolfram.com/QuarticEquation.html
+// and http://www.csit.fsu.edu/~burkardt/f_src/subpak/subpak.f90
+int Factor(double a4,double a3,double a2,double a1,double a0,double
roots[4][2],double EPS){
+ double R[2],D[2],E[2],R2[2];
+
+ if(fabs(a4)<EPS){return Factor(a3,a2,a1,a0,roots,EPS);}
+ a3/=a4;
+ a2/=a4;
+ a1/=a4;
+ a0/=a4;
+
+ Factor(1.0,-a2,a3*a1-4.0*a0,-a3*a3*a0+4.0*a2*a0-a1*a1,roots,EPS);
+
+ R2[0]=a3*a3/4.0-a2+roots[0][0];
+ R2[1]=0;
+ Sqrt(R2,R);
+ if(fabs(R[0])>10e-8){
+ double temp1[2],temp2[2];
+ double p1[2],p2[2];
+
+ p1[0]=a3*a3*0.75-2.0*a2-R2[0];
+ p1[1]=0;
+
+ temp2[0]=((4.0*a3*a2-8.0*a1-a3*a3*a3)/4.0);
+ temp2[1]=0;
+ Divide(temp2,R,p2);
+
+ Add (p1,p2,temp1);
+ Subtract(p1,p2,temp2);
+
+ Sqrt(temp1,D);
+ Sqrt(temp2,E);
+ }
+ else{
+ R[0]=R[1]=0;
+ double temp1[2],temp2[2];
+ temp1[0]=roots[0][0]*roots[0][0]-4.0*a0;
+ temp1[1]=0;
+ Sqrt(temp1,temp2);
+ temp1[0]=a3*a3*0.75-2.0*a2+2.0*temp2[0];
+ temp1[1]= 2.0*temp2[1];
+ Sqrt(temp1,D);
+ temp1[0]=a3*a3*0.75-2.0*a2-2.0*temp2[0];
+ temp1[1]= -2.0*temp2[1];
+ Sqrt(temp1,E);
+ }
+
+ roots[0][0]=-a3/4.0+R[0]/2.0+D[0]/2.0;
+ roots[0][1]= R[1]/2.0+D[1]/2.0;
+
+ roots[1][0]=-a3/4.0+R[0]/2.0-D[0]/2.0;
+ roots[1][1]= R[1]/2.0-D[1]/2.0;
+
+ roots[2][0]=-a3/4.0-R[0]/2.0+E[0]/2.0;
+ roots[2][1]= -R[1]/2.0+E[1]/2.0;
+
+ roots[3][0]=-a3/4.0-R[0]/2.0-E[0]/2.0;
+ roots[3][1]= -R[1]/2.0-E[1]/2.0;
+ return 4;
+}
+
+int Solve(const double* eqns,const double* values,double* solutions,int dim){
+ int i,j,eIndex;
+ double v,m;
+ int *index=new int[dim];
+ int *set=new int[dim];
+ double* myEqns=new double[dim*dim];
+ double* myValues=new double[dim];
+
+ for(i=0;i<dim*dim;i++){myEqns[i]=eqns[i];}
+ for(i=0;i<dim;i++){
+ myValues[i]=values[i];
+ set[i]=0;
+ }
+ for(i=0;i<dim;i++){
+ // Find the largest equation that has a non-zero entry in the
i-th index
+ m=-1;
+ eIndex=-1;
+ for(j=0;j<dim;j++){
+ if(set[j]){continue;}
+ if(myEqns[j*dim+i]!=0 && fabs(myEqns[j*dim+i])>m){
+ m=fabs(myEqns[j*dim+i]);
+ eIndex=j;
+ }
+ }
+ if(eIndex==-1){
+ delete[] index;
+ delete[] myValues;
+ delete[] myEqns;
+ delete[] set;
+ return 0;
+ }
+ // The position in which the solution for the i-th variable can
be found
+ index[i]=eIndex;
+ set[eIndex]=1;
+
+ // Normalize the equation
+ v=myEqns[eIndex*dim+i];
+ for(j=0;j<dim;j++){myEqns[eIndex*dim+j]/=v;}
+ myValues[eIndex]/=v;
+
+ // Subtract it off from everything else
+ for(j=0;j<dim;j++){
+ if(j==eIndex){continue;}
+ double vv=myEqns[j*dim+i];
+ for(int
k=0;k<dim;k++){myEqns[j*dim+k]-=myEqns[eIndex*dim+k]*vv;}
+ myValues[j]-=myValues[eIndex]*vv;
+ }
+ }
+ for(i=0;i<dim;i++){solutions[i]=myValues[index[i]];}
+ delete[] index;
+ delete[] myValues;
+ delete[] myEqns;
+ delete[] set;
+ return 1;
+}
+
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic pop
+#endif
+#if defined(__clang__)
+# pragma clang diagnostic pop
+#endif
+
Property changes on: brlcad/trunk/src/libbg/spsr/Factor.cpp
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/Factor.h
===================================================================
--- brlcad/trunk/src/libbg/spsr/Factor.h (rev 0)
+++ brlcad/trunk/src/libbg/spsr/Factor.h 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,50 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+#ifndef FACTOR_INCLUDED
+#define FACTOR_INCLUDED
+
+#define PI 3.1415926535897932384
+#define SQRT_3 1.7320508075688772935
+
+double ArcTan2(double y,double x);
+double Angle(const double in[2]);
+void Sqrt(const double in[2],double out[2]);
+void Add(const double in1[2],const double in2[2],double out[2]);
+void Subtract(const double in1[2],const double in2[2],double out[2]);
+void Multiply(const double in1[2],const double in2[2],double out[2]);
+void Divide(const double in1[2],const double in2[2],double out[2]);
+
+int Factor(double a1,double a0,double roots[1][2],double EPS);
+int Factor(double a2,double a1,double a0,double roots[2][2],double EPS);
+int Factor(double a3,double a2,double a1,double a0,double roots[3][2],double
EPS);
+int Factor(double a4,double a3,double a2,double a1,double a0,double
roots[4][2],double EPS);
+
+int Solve(const double* eqns,const double* values,double* solutions, int dim);
+
+#endif // FACTOR_INCLUDED
Property changes on: brlcad/trunk/src/libbg/spsr/Factor.h
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/spsr/FunctionData.h
===================================================================
--- brlcad/trunk/src/libbg/spsr/FunctionData.h (rev 0)
+++ brlcad/trunk/src/libbg/spsr/FunctionData.h 2020-10-01 20:42:02 UTC (rev
77311)
@@ -0,0 +1,109 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
modification,
+are permitted provided that the following conditions are met:
+
+Redistributions of source code must retain the above copyright notice, this
list of
+conditions and the following disclaimer. Redistributions in binary form must
reproduce
+the above copyright notice, this list of conditions and the following
disclaimer
+in the documentation and/or other materials provided with the distribution.
+
+Neither the name of the Johns Hopkins University nor the names of its
contributors
+may be used to endorse or promote products derived from this software without
specific
+prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED
WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT
+SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED
+TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH
+DAMAGE.
+*/
+
+#ifndef FUNCTION_DATA_INCLUDED
+#define FUNCTION_DATA_INCLUDED
+
+#define BOUNDARY_CONDITIONS 1
+
+
+#include "PPolynomial.h"
+
+template<int Degree,class Real>
+class FunctionData{
+ bool useDotRatios;
+ int normalize;
+#if BOUNDARY_CONDITIONS
+ bool reflectBoundary;
+#endif // BOUNDARY_CONDITIONS
+public:
+ const static int DOT_FLAG = 1;
+ const static int D_DOT_FLAG = 2;
+ const static int D2_DOT_FLAG = 4;
+ const static int VALUE_FLAG = 1;
+ const static int D_VALUE_FLAG = 2;
+
+ int depth , res , res2;
+ Real *dotTable , *dDotTable , *d2DotTable;
+ Real *valueTables , *dValueTables;
+#if BOUNDARY_CONDITIONS
+ PPolynomial<Degree> baseFunction , leftBaseFunction , rightBaseFunction;
+ PPolynomial<Degree-1> dBaseFunction , dLeftBaseFunction ,
dRightBaseFunction;
+#else // !BOUNDARY_CONDITIONS
+ PPolynomial<Degree> baseFunction;
+ PPolynomial<Degree-1> dBaseFunction;
+#endif // BOUNDARY_CONDITIONS
+ PPolynomial<Degree+1>* baseFunctions;
+
+ FunctionData(void);
+ ~FunctionData(void);
+
+ virtual void setDotTables(const int& flags);
+ virtual void clearDotTables(const int& flags);
+
+ virtual void setValueTables(const int& flags,const double& smooth=0);
+ virtual void setValueTables(const int& flags,const double&
valueSmooth,const double& normalSmooth);
+ virtual void clearValueTables(void);
+
+ /********************************************************
@@ Diff output truncated at 100000 characters. @@
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
_______________________________________________
BRL-CAD Source Commits mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
