Scope
Which entities are considered “geometric”?
What other entity categories are there?
Which geometric entities will be covered?
Data representation and runtime efficiency
Uniqueness.
Caching and redundancies.
Immutability.
Data member protection.
Operator overloading.
Levels of abstraction:
Common base types.
Virtual function interfaces.
Template classes and functions.
Return by value vs. store in reference argument.
Global functions/variables vs. class functions/variables.
Dimensionality
Separate or unified code for 2D and 3D.
Conversions between geometric entities of different dimensionalities.
Extensions to more than 3 dimensions.
Geometric relations and operations
Scope: Given two entity types, T1 and T2,
What subset of relations and operations between T1 and T2 will be implemented?
Which pairs of types will be covered?
Order of operands.
Transformations.
Relations with external packages
cisCommon: RTTI, logging, ... (?)
cisVecs
Serialization and otherwise formatted output.
cisMesh
ITK and VTK
Miscellaneous
Loop operations and iterators.
Representation of errors and exceptions.
A geometric entity is an object with geometric attributes. Geometric attributes are numerical values, such as locations, distances, angles, areas, volumes, etc. That is, any spatial attribute that can be measured. Using this notion, geometric entities include, for example, vectors, lines, rays, line segments, planes, half-planes, polygons, ellipsoids, parametric curves, half-spaces, solids and other three-dimensional objects, etc. The definition we use does not imply any limit on the complexity of the geometric entity.
Other geometry-related entities which are not considered “geometric entities” in this discussion include, for example, transformations (linear on non-liner), quaternions, and operators (e.g. cross- and dot-product, norms, determinants, etc.). To draw the boundary between “geometric” and “non-geometric” entities, we will briefly say that a geometric entity is one that occupies space and stores information about the space it occupies., whereas a non-geometric entity can be anything else.
A geometric entity may be subject to transformations, which change some or all of its geometric attributes. However, a geometric entity usually has some attributes which are invariant to transformations. For example, a triangle subject to an affine transformation remains a triangle. We classify many of the attributes which are invariant to transformations as topological attributes. Examples of topological attributes are the vertex-edge connectivity graph of a polyhedron, or the equal surface continuity of ellipsoids and polyhera.
To help us understand better certain geometric entities, we group them into classes. Classifying a geometric entity allows us to define rules which are applicable to any instance of the class. For example, finding the area of a triangle, or the angle between two lines, involves the notion of the classes “triangle” and “line”, and the procedure for them is identical for any triangle or pair of lines. With this classification, we can use simple geometric entities to describe properties of more complex ones, such as vertices of a triangle, faces of a polyhedron, etc.
Our goal in this document is to design a structured and systematic approach for defining geometric entities. As there are infinitely many possible classes for geometric entities, we cannot cover them all. What we may still be able to do is devise a set of rules for classifying geometric entities, for specifying their essential attributes, and for defining relationships between entities. We can also define a set of elementary quantities and geometric entities, that will allow further definitions of complex entities by applying the rules. When we then seek for a certain entity, attribute, or relation between entities, we should be able to follow the rules again, and find if the desired object exists in our library.
A tentative initial set of geometric primitives could include the following: points and vectors, lines, rays, line-segments, planes, half-planes, and triangles. Note that these are at most 2D objects, and the library should generally be extended to include, for example, more complex polygons and three-dimensional entities. However, this set is rich enough to start with, and to learn about the complexities of defining objects and relations between them.
Considering data representation and runtime efficiency often involves considering tradeoffs. The commonly known tradeoffs between memory, speed, and maintenance, are one set of tradeoffs to be considered here. Another common tradeoff is between abstraction, information preserving, readability, and speed. When coming to resolve a tradeoff, we should bear in mind that speed is never the only concern, and that the cost of bad design increases rapidly with time. On the other hand, while classical, general solutions may be the best asymptotically, we would like to consider some specific solutions for the actual cases that we have to deal with
It is desirable to store a unique representation of the object, so that objects can be compared simply, and operations will yield equal output for equal objects. An example of a non-unique representation for a line is a pair L={Point, Direction}. It is non-unique, as there are infinitely many points on the line, and infinitely many, positive and negative, magnitude values for the direction vector. A unique representation for the line is L={P, D : P is the closest point to the origin on L, D is a unit-length direction vector with last nonzero element positive}. Given any point on the line, and a direction vector of any magnitude, it is easy, though takes a few more CPU cycles, to convert them to the unique representation. However, when we store the line in unique representation, we may lose other properties which we wanted to keep in the first place. For example, if we had two lines intersecting at a point, and we defined them by their point of intersection and their direction vectors, the point of intersection would be preserved in the first representation, and lost in the second. If we had two points, p1 and p2, the line that goes through them will have the same representation, regardless of the order in which they are given, and so we would lose directionality information. The tradeoff for this case can be viewed as specification vs. generalization: we create a general representation of lines, and lose the specific knowledge about the source of information which made us consider the line in the first place. We obtain better abstraction of the class Line, and lose information about the specific Line object.
In general, though, it is a better programming style to use both unique representation and attribute abstraction. The reason is, it generates less coupling between objects, and is therefore easier to maintain. Generally, we would like to assume the least on a Line object, except for those properties that make it a line. We cannot solve the tradeoff for all cases, but we can do well for most. If we still need the information about intersection points, for example, we can use other structures to keep them.
It is often the case that some geometric attributes of an object can be computed immediately as the object is defined, and stored (cached) to be retrieved later. For example, the length of a vector, or the surface normal of a polygon. The same attributes can be evaluated anew each time they are needed. If we choose to cache attributes, we pay with increased memory demands and with increased runtime overhead when the object is created or mutated, as well as with more programming overhead required to maintain the code of the object. We would therefore prefer to store as minimal representation as possible in the object, and, in most cases, store attributes which are used frequently in one particular procedure into a local variable. On the other hand, when the same attribute of the same object is queried frequently in many different locations, and when evaluating it is time consuming, we may wish to cache it in the object itself. Yet, there may be other caching schemes for these cases, such as maintaining parallel data structures – one for the geometric objects and one as a cache – and thus keeping the minimalistic representation of the objects at the cost of having to maintain two separate data structures.
Other aspects of caching have to do with object mutability and with data member protection. When an object is immutable, it is safer and more profitable to cache attributes. When data members are protected, it is easier to measure how frequently each one is accessed, and use that as a criterion to decide whether or not it should be cached.
Immutable objects are safe. They are thread-safe and memory-safe. The price paid for the safety can be increased runtime. Specifically in C++, the copy-constructor mechanism is usually not the most efficient solution we can find, and both creation and destruction of temporary objects are time consuming. But if a variable a represents and immutable object, the only way to assign a new value to a is by creating a new (most likely, temporary) object, then use the assignment operator to store the new value into a. In most cases, there is a fair number of extra CPU operations over what's really necessary, and if we consider runtime efficiency, we would like to reduce them. Even so, there are most likely compiler optimizations for most of these redundant CPU cycles.
As a rule of thumb, I would postulate that mutating complex objects is less frequent than mutating simple objects. The reasoning for that is that we usually work harder to create the complex object, and would like to preserve this work at a minimum cost. This gives us a source of hope.
Another feature of C++ is the ability to block access to the assignment operator, and also for the copy constructor. This enables to declare a completely immutable object, that is, a variable whose value cannot be modified in any way, nor can the variable be re-assigned. In addition to safety, it can also contribute to runtime efficiency, as it would force users of the object to use const references everywhere, instead of copying objects.
The dividing line for mutability should be approximately as follows. Simple quantities, such as Vectors and Scalars, which are subject to many arithmetic operations, may be declared as mutable without any restriction. Objects which are results of operations and returned by value, such as Lines (intersection of planes), Points, etc., should have public copy constructors. It is preferable that the assignment operator be declared private (without implementation) for any class which does not have mutator methods, to keep complete immutability. That means that once an object is created and initialized, it will never be re-assigned. If one has to use temporary variables, one should create and initialize a new object every time a temporary variable is needed. This will make your code that initially looked like this:
cisTriangle tmpTriangle(pt1, pt2, pt3);
/* ... */
tmpTriangle = TransformTriangle(someTransformation, tmpTriangle);
into this:
cisTriangle startTriange(pt1, pt2, pt3);
/* ... */
cisTriangle transformedTriangle( TransformTriangle(someTransformation, startTriange) );
For those concerned about memory overhead, it should be noted that these are all local objects, therefore destroyed at the end of the block. For those concerned about runtime efficiency, it should be pointed out that copying the result of the function call into a vaiable would have taken place anyway, and that an optimized compiler would probably generate no machine code for the copy constructors and destructors involved, so there would not be significant runtime overhead. This scheme is closer in spirit to functional programming languages, and therefore provides cleaner operations. It also enables to backtrack value histories when debugging. One exception to this rule is when a variable is used inside a loop, such as:
for (i = 0; i < n; i++) {
cisTriangle tmpTriangle(a[i], b[i], c[i]);
/* ... */
}
In this case, it is less convenient to track the values of tmpTriangle between iterations with the Visual-C++ debugger, for example. However, note that even in this case, the old values of tmpTriangle are discarded of at the end of the loop, and there is no need to keep the variable beyond it. The only case I can think of where a variable must be declared outside of a loop and re-evaluated through it is with accumulations, such as:
cisTriangle triangle(pt1, pt2, pt3);
for (i = 0; i < n; i++) {
triangle = TransformTriangle(transformation[i], triangle);
}
This can still be solved without using assigment, if we transform the loop into recursion. But very often, we do not want that. Nevertheless, most times we can accumulate the transformation (assuming a Transform object is mutable), and so we get:
cisTriangle triangle(pt1, pt2, pt3);
citTransform transform(); // start with identity transform.
for (i = 0; i < n; i++) {
transform = CompositeTransform(transformation[i], transform);
}
cisTriangle transformedTriangle( TransformTriangle(transform, triangle );
One more thing to note is that overloading the default assigment operator with an identical explicit implementation of memberwise-assignment does not incur cost in runtime. However, writing non-automatic copy constructor might add runtime even when the written constructor is identical to an automatic one, because it may prevent compiler optimizations. Therefore, while rewriting operator=, or a Set() method that replicates a constructor, may be desired in some cases, it is not the same case for copy constructors.
To summarize, we define three levels of mutability: complete immutability (atomic objects), assigment mutability (objects can be re-assigned, but data members cannot be mutated directly), and general mutability (no mutation restrictions). They are prefered in this order. Complex objects are usually not changed once created, and the runtime overhead of creating simple of light objects whenever they are needed, rather than re-assign to existing objects, is usually small enough. General mutability should be reserved to very specific cases, such as accumulators, containers, and iterators. Any mutable object should provide the mutation interface only through mutator methods, and never through public data members. This way, object mutations can be at least potentially monitored, both for debugging and for learned decisions whether the object actually needs to be mutable. Though it may mean some tedious programming for simple objects, such as Vectors, I believe it is worthwhile to follow these principles everywhere, and it pays off in the long run.
As discussed above, data member protection through accessor and mutator methods is a desirable coding pattern. Although it may be initially tedious for programmers to get used to an environment where direct access to data members is disabled, data member protection offers the following advantages:
Consistency: All data members are protected, with no exclusions. There is no question which members should be declared public.
Monitoring: It is easy to find where a member is modified, by placing a breakpoint in the mutator method. This helps detect potential sources for bugs.
Profiling: It is easy to measure the amount of time spent updating each data member, and use that to locate runtime bottlenecks.
Safety: If mutating one member of an object entails other updates immediately, we would like to have the entire mutation as an atomic operation, without worrying about unexpected behavior.
Liner algebra, which is the mathematical basis for our geometric library, defines a set of operators, including various products (dot, cross, matrix), addition, scaling, and maybe others. It is convenient to have a similar set of operators in our library, to reflect our mathematical notion of objects. Using overloaded operators enables to use a more fluent coding syntax, without having to worry about function arguments and method invocation.
However, one should be careful not to over-overload operators. As the syntax of using an overloaded operator is less “formal”, the overloaded operator must be defined in a completely unambiguous way, and have exactly one, well known meaning. Specifically, operator overloading can be useful to represent well-known algebraic structures, such as linear algebra (vectors and matrices), or string operations (concatenation). It should generally not be stretched too far in an attempt to “simplify” unnatural operations.
Another issue concerning operator overloading is discussed under “Global functions, variables”.
Here are some examples for what I consider as misuse of operator overloading, with what I consider as a better solution:
|
Not good |
Maybe better |
|---|---|
cisVec3 v3(x1, y1, z1); |
cisVec3 v3(x1, y1, z1); |
cisTriangle t(p1, p2, p3); |
cisTriangle t(p1, p2, p3); |
cisMatrix m(/* ... */); |
cisMatrix m(/* ... */); |
|
|
|
|
|
|
Sometimes it is desirable to aggregate common features of objects into sets to be used as base classes. Most programming courses or tutorials explain how to use inheritence, and give examples using more or less realistic cases. But I could not find a good source for methodology of deciding when abstraction is necessary, desired, reasonable, or the opposite of any of these. Little abstraction means more work for the programmer when extensions to the library are needed. Over abstraction is when things are too abstract to handle. A good solution is, probably, somewhere between those.
In some cases, some level of abstraction is necessary due to the design philosophy of the library. For example, for most implementations of RTTI, one wants to have a common base-type for all objects (such as a cisObject, or, more concretely, itk::Object, vtkObject, and java.lang.Object in Java). Other examples are generic programming, commonly implemented in C++ through templates, and a large variety of design patterns.
My general approach towards levels of abstraction for our needs is quite loose – use them when it makes sense. It usually makes sense when there is a sufficiently large set of common features (attributes and operations) to justify abstraction. In some cases, one common feature is enough, for example, when we want to have one container to store a set of objects of many different types. In other cases, a few common features may still not make it to an abstract class. For example, individual polygons, which are not related through a mesh structure, can all have GetArea(), GetPerimeter(), and other methods, but since there is no abstract interaction between them, creating the abstraction can inroduce programming obstacles, and may be a waste of time and a source for bugs.
An example of a well-structured and highly abstracted library is ITK. While it is probably useful for many applications, I feel that it is probably too abstract and too large for our needs. Specifically, as of writing this document, ITK does not provide geometric attributes for mesh elements, which are essential for our planned optimization algorithms. Adding these could be troublesome. However, it is a good source of reference for abstraction techniques, and should be studied further to improve our work.
Here are a few issues regarding abstraction, for which I still do not have a definite answer to their necessity in the CIS geometric entities library.
To make each class self-contained, and to increase the portability and modularity of our classes, it is proposed that each class will define its internally meaningful types, in a fashion similar to ITK. That is, all the data and function members of the class will be declared in terms of the class-defined types. Although it requires a more formal programming approach, having class-defined types brings the following advantages:
Readability: it is easier to understand what is the meaning of a parameter
or a return value by looking at a function declaration. For example, a declaration
such as
Line3D::LocalPointType Line3d::EvalAt( const Line3D::ParametrizationType & p );immediately implies that the function takes a parametric representation of a point on the line, and returns a point that is a member of the line (the meaning of LocalPointType is a representation of a point member of the Line object).
Modularity: if we want to redefine a similar set of operations using other underlying types (such as different vector dimensions, integers instead of floats, etc.), most of the changes to be done are to change the class-defined types. Taking this idea further, we can have the geometric class templated by the types of its elements. Taking it even further, we reach the Traits concept used in ITK. But the least we can do is start with the local types, which will simplify these developments.
Portability: following a coding pattern similar to other libraries, such as ITK, will simplify the merging of our code with or into the ITK repository, or with VTK, etc.
The CIS library has a set of standard features that may be applied to any CIS class. These cover the streamed output formatting and level of detail (LOD) control of the cisOSStream package, and partial implementation of RTTI, which is necessitated by the cisOSStream structure. As of March 2003, the implementation of these features is by adding an optional static data member of type cisClass to classes we like to include in this set. That is, there is no full implementation of RTTI by means of deriving everything from a common base class. This decision is reasoned by the wish to keep a subset of lightweight objects, such as Vectors and Matrices, which are used as algebraic entities, and are the core of computational processes that we want to optimize to the maximum. Therefore, we will not assume that there is need for a common base class for all the geometric entities, even though we may conclude it's necessary.
Use of virtual functions offers abstraction and polymorphism, at the cost of adding one implicit member per base class (pointer to the virtual function table), and a slight increase in runtime of virtual methods, which is due to the dereferencing of pointers required for the virtual call. For simple objects, which only take one form, such as Vectors, introducing virtual interfaces would probably introduce more headache than benefit. However, for more complex cases, virtual functions are essential. One such case is use of design patterns, which is not discussed here. Another is the following, which I call “low-level polymorphism”, although it is very close in spirit to the design pattern known as “template method”.
Suppose we want to define a set of geometric entities, whose defining elements may come from different sources. For example, if we take a triangle, the vertices may be given to the constructor as
A set of point values, in which case the triangle is a standalone object.
A set of pointers to memory cells which store the point values. For example, when the triangle represents one facet of a polyhedron, and the vertices are controlled by the polyhedron.
A set of indices into a table that stores the point values. For example, when the triangle is part of a mesh structure, and the mesh stores a table of the vertices.
For each of these cases, the method of computation of area, perimeter, and face normal, for example, work exactly the same, and the difference is in the way they obtain the values of the vertices. We can therefore define our Triangle class as follows:
class Triangle
{
public:
virtual const Point & GetVertex(unsigned int vertexNum) const = 0;
inline Scalar GetPerimeter() const
{
Segment edge0( GetVertex(1), GetVertex(2) );
Segment edge1( GetVertex(2), GetVertex(0) );
Segment edge2( GetVertex(0), GetVertex(1) );
return edge0.GetLength() + edge1.GetLength() + edge2.GetLength();
}
/* ... */
};
class VertexDataTriangle : public Triangle
{
private:
Point Vertices[3];
public:
VertexDataTriangle(const Point & p1, const Point & p2, const Point & p3)
{
Vertices[0] = p1;
Vertices[1] = p2;
Vertices[2] = p3;
}
virtual const Point & GetVertex(unsigned int vertexNum) const
{ return Vertices[vertexNum]; }
};
class VertexPointerTriangle : public Triangle
{
private:
const Point * Vertices[3];
public:
VertexPointerTriangle(const Point * p1, const Point * p2, const Point * p3)
{
Vertices[0] = p1;
Vertices[1] = p2;
Vertices[2] = p3;
}
virtual const Point & GetVertex(unsigned int vertexNum) const
{ return *(Vertices[vertexNum]); }
};
class VertexIndexTriangle : public Triangle
{
private:
const Point * PointTable;
const unsigned int VertexIndices[3];
public:
VertexIndexTriangle(const Point * pointTable, unsigned int p1, unsigned int p2, unsigned int p3)
{
PointTable = pointTable;
VertexIndices[0] = p1;
VertexIndices[1] = p2;
VertexIndices[2] = p3;
}
virtual const Point & GetVertex(unsigned int vertexNum) const
{ return VertexTable[ VertexIndices[vertexNum] ]); }
};
Using this pattern adds some runtime overhead, due to the multiple virtual function invocations in GetPerimeter() and to the creation of the temporary Segment objects. One may choose to optimize the computation of the perimeter to the extreme, by implementing each class independently and getting rid of the Segments. However, the gain in runtime of such an “optimized” implementation would be, in my opinion, very small, especially in the more complex cases, such as VertexIndexTriangle. On the other hand, the unified interface provided by Triangle improves the clarity and readability of the code, and reduces chances of bugs. Furthermore, if a case exists where we need to decide which method to invoke, the virtual function mechanism is superior to any conditional evaluation.
Note that the pattern also implies protected data members, avoidance of caching, and, most likely, immutability. As we concluded before, these are generally desired features. With a little extra work, uniqueness can also be achieved.
C++ templates are means for generic programming. A very useful example of templates is the C++ standard library, known mostly as the Standard Template Library (STL). Use of templates, when appropriate, generates code that can be both highly efficient and general. Generality can be expressed, for example, in handling different data types, different dimensions, and even different definitions of relations (e.g., order relations), which are given as templated functor objects.
A drawback of template use is that it requires a very high level of abstraction. In some cases, solving a problem for the general case is much harder computationally than for the specific cases we have to deal with. In other cases, understanding the general case terms and developing solution methods require much more work than providing an immediate solution to an immediate problem. In yet more cases, learning how to use a templated algorithm can be hard, as it requires the user to acquire all the underlying concepts, in which case it might have been easier for the user just to write the solution by himself.
My feeling is that at the moment, the development of the CIS library should not concern with templates. Firstly, there is already a rich templated library, ITK, which will always be ahead of us in developing templated solutions. Secondly, our goals involve practical real-time programming and developing answers to specific problems, rather than solving general cases. Thirdly, while we do wish to own a reliable, consistent, and extensible infrastructure library, we also want it to be simple enough for new students to learn. I believe that as of now, ITK and templated programming are too cumbersome for these needs. While solving general-case problems and developing algorithms using templates are interesting research topics, they are not in our focus of research. Nevertheless, use of existing templated code, such as STL, is often a good idea.
The question, “How to store the value(s) computed in a function?”, typically has two answers: return value and store in a reference argument. As one main theme in this document is object immutability, we would naturally prefer the first answer. That is, when an operation evaluates a new outcome, we will return this outcome by value. However, for efficiency concerns, and for cases where the computed value is not represented by a single object, we would like to outline exceptions to the rule, and ways to handle them.
“Store in place” operations. For operations involving accumulators, it is often more efficient to store the computed value directly into the target accumulator rather than create a new object to represent the return value, then copy the new object into the target. For example, the operations v += u and v *= M are probably more efficient than v = v + u and v = v * M. Note, however, that most geometric entities (lines, planes, shapes, ...) are not accumulators. Rather, their components may be. Although we could still represent in-place operations on them, it seems to me reasonably preferable to accumulate operation results into true accumulators (such as Vectors), and store the final Vector values into fresh geometrical objects. Though one may wish to “simplify” this work by allowing in-place operations on geometric entities, one may end up having to implement many more functions than necessary to achieve this goal. In my opinion, it is better to limit the set of permitted operations to a logical minimum, and provide necessary, and possibly efficient conversions between geometric entities and sets of vectors.
More than one return value. Consider an operation such as GetClosestPointPair(const Line & l1, const Line & l2, Point & p1, Point & p2), which evaluates p1 and p2 as the closest point on each of l1 and l2 to the other line respectively. For this case, there are two possible solutions which do not involve two return values:
Segment Line::GetShortestSegmentTo(const Line & l2) const; and
Point Line::GetClosestPointTo(const Line & l2)
const;
Point Line::GetPointProjection(const Point & p) const);
While it remains to be questioned why one would use the pair of closest points without considering the segment between them, there are still cases where a function returns a tuple completely unrelared objects, and further yet, we do not wish to create an artificial class just to accommodate this tuple.
We could, for example, have the following definition:
class EmptyClass {};
template<class X, class Y = EmptyClass>
class Tuple{
public:
typedef X FirstType;
typedef Y SecondType;
FirstType & GetFirst() {
return First;
}
SecondType & GetRest() {
return Rest;
}
/* ... */
};
The second type of a Tuple can be another Tuple, much like Lisp programming style. Then a function can return a Tuple of any object types. This imposes cumbersome syntactic structures on the programmer. In particular, returning a Tuple by value may incur significant runtime overhead in copying all members of the Tuple. For those cases where there is no logical structure to the return values of a function, we may be forced to use referenced result variables. We should adopt a convention that places all the return values last (or, alternatively, first, but consistently), and possibly names them resultXXX to identify them as result values.
Containers. Returning a container by value is generally extremely inefficient. When writing a function that modifies a container argument, and the container is the only modified object, we can adopt the following convention:
ContainerType & ModifyingFunction(ArgTypes args, ContainerType & result)
{ /* do something */
return result;
}
Runtime Exceptions. These are to be discussed elsewhere.
Efficiency considerations. As long as the returned value is a small object, that is, fixed-size with not too many fields, we should not worry much about copying overhead. For start, storing values into a referenced variable takes its own time. Then, when using automatic copy operations, they are likely to be optimized by the compiler. Large, complex object would usually fall under a container category, which is handled explicitly. For everything else, we may have to make exception to the rule.
In general, global functions should be avoided. If a function has one main operand, it can be made an object method. If it has symmetric operands of equal type, it can be made class static method. For any case where a function cannot be associated with any object class, it is better to have it as a static method in an artificial class, made to aggregate related functions which are defined in one C++ file. This way, it is at least easy to track down the location of the function code by looking at its class name.
One case that cannot yet be resolved with static functions it overloaded operators which take the form of global function. For example, ostream & operator <<(ostream & out, const MyClass & obj);. Global operators defined this way should be coded in the same C++ file as their operand that does not belong to the standard C++ family. When both operands are new types, the operator should be defined in the same C++ file as the type of larger size out of the two. For example, Matrix * Vector is to be defined in the Matrix file, and so is Vector * Matrix, as long as the latter is declared as a global function.
Generally, I believe that except for in-place operators (such as +=), operators which return new objects are better defined as global functions. This has to do with symmetry and consistency issues. For example, Scalar * Vector is a very useful notation, and must be defined as global. Symmetrically, Vector * Scalar should also be global. The globality of other operators follows a similar logic.
Like global functions, global variables should be avoided. Global constants of specific class types, such as XVector, should be defined as static members of their class: const cisVec3cisVec3::XVector(1, 0, 0);. Global constants which relate to attributes of classes (such as the number 3 representing the dimension of a Vector) should also be declared as class constants. I believe that global variables can be completely eliminated this way. That is because any variable denoted a value of some attribute, and therefore can be associated with a class. There are other substitutes to global variables, such as Singleton objects, etc., to take this approach even further.
Different space dimensions offer different varieties of geometric entities. m dimensional entities defined in an n dimensional space (m <= n) are still n dimensional entities, and are given by n dimensional coordinates. To deal with spaces of different dimensions, it is necessary to define consistent projection rules between them.
As the underlying vector types of spaces of different dimensions are different, we cannot assume, for example, that the two-dimensional XY plane is a subset of the three-dimensional XYZ space. Rather, we say that the two-dimensional plane can be mapped into the 3D space, and that 3D objects can be projected on a plane. It is important to remember that there is no algebra that involves directly pairs of vectors of different dimensions without a projection rule (e.g., in the form of m*n matrices). Furthermore, the algorithms to evaluate certain relations, such as finding orthogonal vectors, may be different in spaces of different dimensions.
That said, we conclude that it's important to completely separate geometric entities of different dimensions, and restrict our operations to objects of equal dimensionality.
It should be asked if we should develop one package to accommodate all geometric objects of all dimensions, or a separate package for each dimension. The answer has to do more with the development process than with abstract definitions. As we gave up the dimensionality abstraction level (which exists in ITK, for example), it seems to me that we are bound to provide specific solutions for the problems of each dimension. That means we can separate the development process, and build classes incrementally and independently for each dimension. Of course, it's a good idea to have as similar interfaces as possible, so that if ever we wish to abstract dimensionality, we can do it easily. For my convenience, I would develop the 3D geometric objects without regarding the 2D case. Those who wish to extend the library may do it according to their needs.
Extension to more than three dimensions is currently out of the list of our goals. High dimensions are diffcult to handle. ITK is templated over the dimension of the geometry, but on the other hand, it does not contain geometric entities.
The set of operations between geometric entities (in the same space) can be classified as follows:
Projections
Parametric evaluations
Closest point
Intersections
Angles
Membership tests
Projections map objects of high dimensions into objects of lower dimensions. The returned value can be given either in high-dimension global coordinates, or in low-dimension local coordinates. For example, the global projection of a 3D point on a line in 3D space is a 3D point on the line. The local projection of the same point on the same line is a scalar that parametrizes the projection point with respect to the unique representation of the line. Similar relations can be defined for other cases.
Parametric evaluations take a local coordinate and return the global coordinate evaluated with respect to the target object. Certain entities may have more than one parametrization. For example, simplexes can use spatial or barycentric coordinates. A plane in 3D must have its origin and local axes well defined for parametric evaluation.
Closest point is usually identical to projection when the target object is an unbounded entity. When the target is bounded, the closest point may be a boundary point, for which special tests need to be performed. When the target is not convex, the closest point may be non-unique. As our list of entities currently contains only convex objects, we will not deal with this problem now.
Intersections of objects generally return objects more complext than a point. If the return value can be classified in our library, we can implement the interaction operation. Otherwise, it is likely to turn more and more complex.
Angles must be defined uniquely. For two unit vectors, the angle is defined as acos( v1 * v2 ), which returns a value between zero and pi. As lines are defined uniquely with restrictions on the last element of the direction vector, the range of angles between lines may be affected. For example, we will not be able to represent lines that have the same point but opposite directions. We will be able to do it with rays. But how should we choose the angle between a line and a ray? A similar problem may exist for lines and planes. This issue should be discussed further.
Membership tests are subject to numerical tolerances. For example, testing if a point is on a given line involves measuring the distance between the point and the line, and deciding on a threshold under which we return “true”. Generally, we cannot assume infinite precision, or zero distance, even if the point was obtained by a parametric evaluation of the line.
For a well defined operation, we will follow this convention. A projection or closest-point operation always returns a point on the target object. That is, we can have definitions such as
Point Line::ProjectPointGlobal(const Point & p) const; // return a point on this line
Scalar Line::ProjectPointLocal(const Point & p) const; // return a parameter for evaluating a point
Line3D Plane3D::ProjectLineGlobal(const Line3D & l) const;
Line2D Plane3D::ProjectLineLocal(const Line3D & l) const;
Point Segment::ClosestPointToLineGlobal(const Line & l) const; // return the closest point on this segment to the line, in global coordinates.
Scalar Segment::ClosestPointToLineLocal(const Line & l) const; // return a value in [0,1]
Parametric evaluation takes as input the result type of the local projection operation:
Line l;
Point p;
l.EvalAt( l.ProjectPointLocal(p) ) == l.ProjectPointGlobal(p);
Plane pl;
pl.EvalAt( pl.ProjectLineLocal(l) ) == pl.ProjectLineGlobal(l);
For intersections, the operations will be defined for the object of higher dimension, with an argument of lower of equal dimension. For objects of equal dimension, the operation will be on the less bounded object, and the argument will have less or the same amount of boundaries. If both objects are unbounded, such as a line and a curve, the operation will be on the one that is defined by a larger set of parameters, i.e., the curve. The general idea is that the operation is on the “bigger” object.
Segment::Intersect(Segment);
Ray::Intersect(Segment);
Ray::Intersect(Ray);
Line::Intersect(Segment);
Line::Intersect(Ray);
Line::Intersect(Line);
Triangle::Intersect(Line);
HalfPlane::Intersect(Triangle);
Plane::Intersect(HalfPlane);
Ellipsoid:Intersect(Plane); // ellipsoid is defined by three orthogonal axes and three magnitude values
Tetrahedron::Intersect(Ellipsoid); // tetrahedron is defined by four vertices
...
Transformations are not geometric entities. Our main goal here is to define a standard way for applying transformations to geometric entities. For this discussion, we will assume the existence of a method Transformation::Apply(Vector) for a single vector. We need to question transformation invariants. That is, if a transformation does not preserve the class-defining properties of an object, is the application of the transformation to that object defined at all?
We can see three syntactic structure for applying a transformation to a geometric object:
ObjType2 Transformation::ApplyTo(ObjType1)
ObjType2 ObjType1::ApplyTransformation(Transformation)
ObjType2 o2( transformation.ApplyTo( o1.ExtractFeaturesToTransform() )
In the general case, we do not know if the output of the transformation has the same type as the input. The third structure uses an intermediate data type to store, for example, a set of vectors, to which the transformation is applied in a standard fashion. The output is another set of vectors, which are used to initialize o2. The underlying assumption is that the objects can be defined completely by sets of vectors, or in another standard way. That may not always be the case.
We also need to question which transformations will be applicable to which objects. We cannot hope to cover the full matrix of transformation types and object types. Therefore, we will have to design the application syntax in a way that is simple and easilty extendable. This is yet another good reason to get rid of operator overloading when applying transformations to other objects than the simplest ones.
My feeling is that we should use the second structure, that is, have the target object be the subject of operation, rather than the operation be parametrized with the object. I think this approach is cleaner and more extensible than the first. Just as we want a minimal representation of the geometric object, we also want a minimal representation of the operation. As the amount of target object types may increase beyond the amount of available operations, and we do not expect to fill the matrix, I'd say it's safer that each object specify which operations are applicable to it, than each operation specify to which objects it can be applied.
The concepts to be discussed in this section are more general than the ones before, as they relate to larger-scale development philosophy. While an attempt was made to outline such philosophy for the geometric entities library, defining one for the CIS project requires a more joint effort. I will try to outline here a few subjects for further discussion.
The cisCommon package is supposed to provide common functionalities and features for all other CIS packages. These can include standard message logging, type-information structures (RTTI), basic type and constant definitions, and more. The scope of the cisCommon package needs to be determined, and its rules applied to any other package of CIS.
cisVecs provides the basic computational elements for the geometric package. Therefore, the geometric package depends directly on cisVecs. Many of the guidelines defined here can be applied to the cisVecs package. But the design of cisVecs is not withing the scope of this document.
I would strongly recommend the design of a global serialization scheme for the CIS library objects. Existence of serialization interface should be a necessary condition for a class to be included in CIS. Serialization is needed for distributed systems, and also for planning systems (such as surgical planning), when the planner wants to save the plan for later editing or use. This document did not specify a serialization scheme, as the scope of such scheme is beyond the document's goals.
In addition, it is very useful to include a standard scheme for output formatting. The overloading of operator<<(ostream &, ...) by itself is insufficient, first because operator<< is not a formatted output (while C printf(), for example, is), and second, because we might want to format the output regardless of a stream object (e.g., to display on a pop-up window). A standard method, such as ToString() could help much more than any other solution. My belief is it should exist for any CIS class.
As of the creation of this document, the borderline between cisMesh and the geometric package is not clear. cisMesh contains several geometric objects, such as triangles and tetrahedra, which in my opinion should be in the geometric package. On the other hand, there is lack of functionality in the existing cisGeomObj package, which may be needed for mesh operations.
A mesh object should be built on top of a powerful geometric library. A mesh is, in the general sense, a container of geometric objects, and can provide higher-level operations that involve groups of mesh elements (usually called “cells”). A good insightful example to is the itk::Mesh class, and its related classes. Theirs is very generic, and we could probably do better with a simpler structure. Yet we can use many of their ideas.
The cisMesh package is yet to be defined in the future.
These two projects compete, in a way, with our internal CIS code. They provide reliable and rich infrastructures for many types of applications. One of VTK's strengths is the built-in visualization tools, which eliminate the need to create our own interface to OpenGL. Another is its very structured nature. ITK's main strengths are its genericity and the built-in image processing tools.
Formerly, Jianhua Yao has written his own interfaces for both OpenGL and image processing. I would like to replace them with external, reliable libraries. For example, Yao's implementation of DICOM readers was very limited, and did not handle well all the cases. Rather than fix it, I would replace it with ITK's (or VTK's) DICOM reader, and get better compliance with the DICOM standard guaranteed. Likewise, to visualize complex structures, I would rather follow the well-defined rules of VTK, than try to devise my own rules, or to re-implement VTK's.
With all this, I feel that ITK and VTK are, at least at the current stage, an overkill for our projects. In addition to efficiency, I believe we strive also for simplicity and for ease of use and learning. These do not always exist in ITK and VTK. In addition, I could not find some features I needed for my research in either of them, and I would not be so happy about trying to write my ideas in the ITK language from the beginning.
My perception is that CIS can offer a simple testbed for ideas which are not part of ITK and VTK. After testing them, we can launch projects for implementing them in either, or let the ITK/VTK developer groups incorporate our paper results in their libraries. Gradually, I believe, CIS will form a structure that will be comparable with ITK or VTK. When that happens, we will be able to consider porting CIS into one of them, or merging the development process of CIS with them. At this moment, I do not think we are in that stage.