On Tue, Oct 14, 2008 at 8:20 PM, Simonson, Lucanus J <lucanus.j.simonson@intel.com> wrote:

...

- Why can objects be default-constructible? ... - centroid says it throws an exception if the polygon does not contain any point. But how can a polygon not contain any points!? The constructor of a polygon should force it to have at least one

Mathias wrote: point.

Obviously, a default constructed polygon contains no points. Polygons have to be default constructible so that they can be elements of a std container. Objects that are not default constructible are a pain. I have to deal with a legacy polygon type that is not default constructible, so I have a very clear idea about which way I prefer and why.

STD library containers don't need DefaultConstructible, only CopyConstructible and sometimes Assignable. Could the algorithms be written in such a way that they didn't care whether the polygon type was DefaultConstructible? Is there a different motivating reason to allow empty polygons than dealing with a legacy class? --Michael Fawcett

Simonson, Lucanus J wrote:

...

- Why can objects be default-constructible? ... - centroid says it throws an exception if the polygon does not contain any point. But how can a polygon not contain any points!? The constructor of a polygon should force it to have at least one

Mathias wrote: point.

Obviously, a default constructed polygon contains no points.

Yes, those two comments were indeed linked to show they were maybe inconsistent. Also, if there is a valid reason for polygons to be able to be empty, I think the algorithm shouldn't make the check but assume it as a precondition, possibly with a disactivable assert.

Polygons have to be default constructible so that they can be elements of a std container.

Thankfully, standard containers do not require objects to be default constructible, but only copy constructible. I'd personally be even happier if they required almost nothing at all, which hopefully they will.

Objects that are not default constructible are a pain.

Says who? Two-phase initialization is well known to be a pain. What's the point of having an object constructed that is not really initialized and usable? Constructors are made to enforce invariants. For some objects, it doesn't make sense to make to allow "emptiness" or "default" as a valid state, so we might as well not allow those states to exist. But anyway, that's another debate. I think it was discussed extensively lately in comp.lang.c++.moderated, so there is no need to discuss it here.

I have to deal with a legacy polygon type that is not default constructible, so I have a very clear idea about which way I prefer and why.

Dealing with non default-constructible types is no more difficult than dealing with default-constructible ones. If you have a very clear case where for polygons or something related to geometry it can lead to issues, feel free to expose it.

Just a few additions from my side.

Here are the few answers I can give right now, I let Barend complete or correct me, and answer the questions for which he has a better view of the library than me.

...

- Still about "box", it is really an axis-aligned (whatever that means in non-cartesian geometries) hyperrectangle. Since "box" can not contain an arbitrary box, maybe it should be called minmaxbox or aabox or whatever.

aabox could be good.

...

- envelope calculates the minimum bounding axis-aligned hyperrectangle for an object. It might also be desirable to have a way to calculate the minimum bounding hyperrectangle (non necessary axis-aligned) or the minimum bounding sphere. What names would those algorithms have if envelope is reserved for AABHR?

The names of geometric objects are, as much as possible, taken over from OGC (see also http://geometrylibrary.geodan.nl/ogc.html). Those names (point, linestring, linear ring, polygon) are also used in Google's KML and in most databases with geometry support. Indeed the "circle" and the "box" were not there, we've added them, and they are open for discussion. Also the names of the algorithms are borrowed from OGC, and are used in many GIS applications and spatially enabled databases. It seems therefore useful to have an algorithm "envelope" there. See for example http://dev.mysql.com/doc/refman/5.1/en/geometry-property-functions.html#func... (which, as OGC describes, returns a polygon - we thought a box was more appropriate for C++ applications, but it is open for discussion).

...

- Modeling sequences of polygons as iterator pairs seems an annoyance to me, ranges are much prettier to handle since one variable is thus sufficient to represent one object (while two are required for iterators).

Maybe we can double each algorithm to have the choice between passing a pair and a range. AFAIK, this will be the approach of the algorithms in C++0x, am I right?

None of the algorithms use iterator pairs for polygons. Polygons might have holes, in our library, so having just a sequence of points is not sufficient to define a polygon. On the other hand, algorithms take linestrings by iterator pairs, which was in fact a positive result of the discussions of this list. Return of linestrings / polygons is done using output iterators.

...

- Algorithms should be overloaded to handle all geometric objects. I assume the reason it wasn't done is simply because of time to do so.

Yes I suppose that we'll be able to add those overloads once the general design will be validated.

why isn't envelope returning in result rather than taking a reference? The reason is, like with centroid, that the output box / point might have a different type than the input geometry. If it was a result, the

- Why are intersection_segment arguments template parameters instead of segment, why is it not simply an intersection overload and why does it return something different than intersection? Maybe the design of the generic intersection should be extended so that it can work with segments rather than segments be made a special case. Right. The case of the segment, in fact, is actually an internal implementation function. It should have been in namespace "impl", will be changed. Agree, all intersection algorithms should return the same

Agree, where it makes sense. The "within" algorithm is not symmetric, a point can be within a polygon but a polygon cannot be within a point. library user has, when calling the function, to specify the type as template parameter. But it might be changed, it is open. thing and should have the same name. Internally it is useful that it exactly returns how it intersects.

- An overlap algorithm, that would be equivalent to testing not_empty(intersection) || within, would probably be interesting, especially since it is very useful for collision detection or spatial indexing. It can obviously be much faster than performing intersection and within calls. Agree. - centroid says it throws an exception if the polygon does not contain any point. But how can a polygon not contain any points!? The constructor of a polygon should force it to have at least one point. A polygon with only one point has also not much meaning. See also discussion about construction. As soon as the polygons follow the concept Polygon, the construction is open to the implementator / library user. Then the centroid algorithm can still get an empty polygon.

Barend

Barend Gehrels wrote:

None of the algorithms use iterator pairs for polygons. Polygons might have holes, in our library, so having just a sequence of points is not sufficient to define a polygon.

I didn't say for polygons, I said for sequences of polygons (what used to be called multipolygons I think). It's true however no algorithm take these at the moment, so let's consider sequences of points, aka linestrings.

On the other hand, algorithms take linestrings by iterator pairs, which was in fact a positive result of the discussions of this list.

It is good, but ranges would be more practical and less verbose. Instead of passing two arguments (the begin iterator and the end iterator), you only pass one, a range (and by const-reference, not by value). A range is an object from which you can extract a begin and an end iterator (using boost::begin and boost::end). A std::pair of iterators is a valid range, but so is also any container. Consider std::for_each, here how you make it range-aware: template<typename Range, typename F> F for_each(const Range& r, F f) { return std::for_each(boost::begin(r), boost::end(r), f); } That way, you can directly pass a linestring, a vector of points, a deque of points, a lazily computed adaptor or whatever to the algorithms with concise syntax: ie. just 'foo' instead of 'foo.begin(), foo.end()'.

Agree, where it makes sense. The "within" algorithm is not symmetric, a point can be within a polygon but a polygon cannot be within a point.

If the case can never happen, why not make an overload that simply returns false? If I have generic code and I want to know whether an instance of T1 is within an instance of T2, I don't really want to make special the case where T2 is a point before calling the generic algorithm.

Barend wrote:

Besides this, in the past discussions concerning the previous previews,

and the other submissions, discussions concentrated on the points. How it was modelled and how it should look like. Until now I didn't saw any

comments on this (besides construction of generic geometry objects). Can we conclude that the list is now satisfied with the provided point concept?

template <> struct access<example_point_1> { template <int I, typename PP> static typename forward_const<PP, typename support::coordinate<example_point_1>::type>::type& get(PP& p) { return I == 0 ? p.x : p.y; } }; No, not really. First of all, what about constructing and setting values of the point? It appears to be missing from the generic interface. Also, your get function returns a reference to coordinate. It should return by value. Many legacy point types may have a public API that prohibits direct access to their members because that is what good little C++ programmers were taught to do in school. Those point classes could not conform to your concept. A polar coordinate point could give a Cartesian coordinate view of itself, but not if it has to return is Cartesian coordinate values by reference. It could not conform to your concept either. A point is not a tuple. I don't think it is valuable that tuple conforms to your point concept. An interval has a low and high value, could be implemented as a tuple of size two, and would incorrectly conform to your point concept. That's not a good thing. If you had an interval concept in your geometry you would have to ensure that it was different from the point concept in some way. Since you are defining traits for your geometry objects, why not define their conceptual type as a trait? template<> struct geometry_concept<example_point_1> { typdef point_concept type; } This takes the guess work out of things and allows you to do tag dispatching based on conceptual type. Luke

On Tue, Oct 21, 2008 at 9:24 AM, Bruno Lalande <bruno.lalande@gmail.com> wrote:

Luke, I come back to you about this point from you:

...

Also, your get function returns a reference to coordinate. It should return by value. Many legacy point types may have a public API that prohibits direct access to their members because that is what good little C++ programmers were taught to do in school. Those point classes could not conform to your concept.

[snip]

If you are talking about using the getter for write access (that is: get<0>(p) = value) it's more problematics. We've discussed about that with Barend and it appears that as you say, a legacy point without direct access to its coordinates will never be adaptable to such traits. I hadn't thought about that before.

Could get<0>(p) return an object (convertible to the coordinate type, so it can also be used for read access) with an overloaded assignment operator that does the appropriate thing for write access in such a point class? It seems like that would solve the problem, but I'm not sure about any drawbacks. Kind regards, Stjepan P.S. kudos for all the geometry library development and discussions. I'm looking forward to a well-designed geometry library in boost.

Luke wrote:

...

Also, your get function returns a reference to coordinate. It should return by value. Many legacy point types may have a public API that prohibits direct access to their members because that is what good little C++ programmers were taught to do in school. Those point classes could not conform to your concept.

Bruno wrote:

If you are talking about using the getter for read access, this way of writing traits specialization is not an obligation and the user can always return by value in his own one.

Except that wouldn't work with library code that used it as an lvalue. If you are talking about using the getter for write access (that is: get<0>(p) = value) it's more problematics. We've discussed about that with Barend and it appears that as you say, a legacy point without direct access to its coordinates will never be adaptable to such traits. I hadn't thought about that before.

The problem is that this approach (value = get<0>(p) for read access and get<0>(p) = value for write access) has been chosen after having been advocated by some people on this list. The argument was that it is a very common approach in Boost, for Tuple and Fusion. And this argument was quite convincing.

Must I deduce that the Boost.Fusion approach is not totally generic and that a structure with unreferencable private members couldn't be adapted to it?

And what would finally be the best way to handle point write access if the Boost.Fusion one is wrong?

I don't think it's a matter of right and wrong. It's a matter of goals and objectives. My goal is to work with all geometry types and not impose unnecessary restrictions on the way geometry types need to be implemented. Fusion is a way to work with plain old data. Because point objects have no invariant to enforce, they should be plain old data (though often aren't defined that way), but that is not true of geometry objects in general. 1D intervals are almost as simple, but do have an invariant, and therefore should not provide direct access to their data members. I was just looking a Brandon's library, and it seems he has no write access to point, but does have a construct function that accepts the coordinates and returns a point by value. This could be used to overwrite a point type that provides an assignment operator, but would be cumbersome to use to write to a single coordinate of a point. My traits looks like the following: template <typename T> struct point_traits { typedef typename T::coordinate_type coordinate_type; static inline coordinate_type get(const T& point, orientation_2d orient) { return point.get(orient); } static inline void set(T& point, orientation_2d orient, coordinate_type value) { point.set(orient, value); } static inline T construct(coordinate_type x_value, coordinate_type y_value) { return T(x_value, y_value); } }; This is specifically the Cartesian 2D point traits, the 3D point traits are a different structure. It needs to be specialized for point types that don't conform to my expected default behavior. Like Brandon, I provide a construct function that returns a point by value. I also provide the set function that takes a point by reference, an orientation class encapsulated enum and the coordinate value. I prefer runtime indexing, as I've mentioned before. There is a certain degree of redundancy in the way I define the traits, since construct and set are both ways of getting data into an object. I think that Brandon may have the right idea with a getter function and a construct function. It is rare to write to only part of an unknown object type. I conclude that the point accessor should return by value or const reference, and that the traits should also provide a construct function to construct a point from coordinate values. A given point class should additionally provide an assignment operator and copy constructor to conform to the point concept. If we are able to get away without requiring default constructor that would be a good thing. For an ND point, the construct function should take an iterator pair over coordinate values and the accessor function must be runtime indexing. In looking more at Brandon's code, I have to repeat my observation that he has done an impressive job of defining a generic geometry type system and using boost libraries to do it. I think that if we make his accessor dimension agnostic (either runtime or compile time indexing) we have a winner. Regards, Luke

Simonson, Lucanus J wrote:

Also, your get function returns a reference to coordinate. It should return by value. Many legacy point types may have a public API that prohibits direct access to their members because that is what good little C++ programmers were taught to do in school. Those point classes could not conform to your concept.

This is not a problem, you just have to write a proxy. That kind of thing is done very often to achieve operator[] overloading with different semantics for reads and writes. It's also done a lot for iterators, which also return references.

Simonson, Lucanus wrote:

Here I disagree. First of all, it is an absurd amount of work.

Yes, you have to make an overload per combination for optimality. However, since most concepts are refinements of polygon, handling polygon in polygon should get you support for mostly everything but spheres and interpolated stuff, but not in an optimal way.

Fully generic functions cannot be overloaded directly, so you would have to use tag dispatching.

While that solution should work fine, I suppose using SFINAE (which might be considered "direct overloading") would be best.

You can't specialize it either because that would mean rewriting the entire library for every new type. I want to emphasize that we are not talking about making a boost::point and boost::polygon objects and writing a bunch of functions that accept those types. We are talking about writing fully generic functions that work with anyone's point and polygon objects.

Overloading over concepts is not much more difficult than overloading over types using SFINAE. There are no best matches with SFINAE however. (there may be ways to emulate that, any one knows of any trick?)

Second of all, if you write code that is doing something nonsensical it is more likely that you are doing something wrong conceptually and would prefer a compiler error to remind you that a polygon cannot be within a point.

I don't see how checking if an object is within another, be the combination impossible to be true, nonsense. The type of your object may not be known by your (meta-)program and it may need to know whether that information is true. Would you rather have a trait can_be_within<T1, T2>::value than just overloading within to return false? While it could be useful, requiring to check that before using within is quite a big hassle IMO.

The purpose of templates is not to make every line of code the developer could write compile by breaking down the type safety of C++.

You'll have to show how it is type unsound.

It is desirable to provide conceptual type safety for geometry types such that you cannot pass point and polygon to the within() function in the wrong order and get a bug that you have to figure out during testing.

You are assuming you're never using 'within' from a generic context. It seemed obvious to me a geometry within boost should be generic, and enable people to abstract the type (or meta-type/concept) of objects they're dealing with. If you have to special case everything, then you can't make a generic program. As for 'within', note also that any object may be in a point if you consider bounds and if that object is reducible to that same point (a polygon with all points being the same or having only one point, a sphere with zero radius, etc.). Which is quite funny if you think about it: the point is actually the most refined concept, since it's a polygon, a sphere, a line and everything at the same time, even though it's the basic building block.

Box has a dimension-agnostic flavor to me (if not to all geometers out there). I believe the orthodox term would be (axis-aligned) parallelotope, some people also use hyperrectangle. Box is shorter and implies the pairwise orthogonality of the sides, if not the axis- alignment. Personally I'd prefer to keep box for most common use, with a clear definition in the docs, and polytope (most general) and parallelopiped (2^n sides not necessarily orthogonal) for the more exotic uses... My two cents... Sent from my iPhone On Oct 15, 2008, at 11:23 AM, Jeff Flinn <TriumphSprint2000@hotmail.com> wrote:

...

...

- Also, while most names are 2D-centric, "box" is 3D-centric. This seems fairly inconsistent. Or maybe it's just me and box is actually dimension agnostic. box doesn't sound less dimension agnostic than anything else to me, but this can be discussed...

Convex Hull?

Jeff

_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

Back to the issue of default vs. non-default constructible polygons, I believe that while the boost geometry library may choose to define a polygon data type that may not be default constructible, it cannot make that choice for legacy polygon types. We have to assume that the library's generic interfaces will work with polygons that are default constructible and that the library's algorithms need to be implemented in a way that is robust to all kinds of things that might normally be considered an invariant of a well designed polygon class.

Indeed, my preference towards non-default constructors is available for a class but becomes much less ideal when it comes to generalize it to a concept. In that precise context, runtime checks inside the algorithms are probably preferable. Bruno

Simonson, Lucanus J wrote:

How can you model a rectangle that is not axis aligned and is still a rectangle? The only way I can thing of is to model it is width, height, center point and rotation. If you snap its corners to the coordinate grid it will, in general, no longer be a rectangle. Similar to the circle, I like rectangle, prism and hyperprism for 2D, 3D and ND rectangles. While non-axis-aligned rectangles could be useful, you have the option of modeling them as polygons with 4 vertices that are as nearly rectangular as you can make them within your coordinate system.

In my domain we call this an oriented bounding box. Typically in 2d its represented by a center point, x and y axis vectors, width and height. Something like: template< typename T > struct obb { point<T,2> center; vector<T,2> x_axis, y_axis; T width, height; }; You could also represent it at 3 points I think. Although I suppose in your domain where you use mostly integers these definitions might not be workable. -- Michael Marcin

Simonson, Lucanus J wrote:

...

My preference would be to have a circle data type, a sphere data type and a hypersphere data type. Only the hypersphere would be parameterized on dimension. The user would pick the object that is named for what they want. If a user wants a circle they are going to be frustrated and confounded that we provide them nsphere<2>. There is something of a conceptual disconnect between circle and sphere.

Steven Watanabe:

What is wrong with typedef nsphere<2> circle;

Nothing. I'm not against saving typing with templates, I just don't want to do it at the expense of readability and intuitiveness of the library API. Some compilers will report the typedef name circle in errors and warnings, while others will report nsphere<2>, which could hurt readability of compiler errors. Also, there is more at issue here than what to name the geometry objects the library provides. The more important issue is whether there is a separate circle concept, sphere concept and nsphere concept. Circle and sphere could potentially be related to nsphere in some manner, both circle and sphere might satisfy nsphere, for instance, but not necessarily the other way around. We might use inheritance or specialization to produce a certain effect with concepts. Something that satisfies circle would also satisfy nsphere, because circle could be derived from nsphere, but nsphere would not satisfy circle unless it was 2D. I was thinking something like the following might be a good idea: template<int n> struct hypersphere_concept { ... }; struct sphere : hypershpere_concept<3> { //optionally add stuff specific to sphere concept not applicable to nsphere ... }; struct circle : hypersphere_concept<2> { //optionally add stuff specific to circle not applicable to nsphere ... }; This general pattern should be applicable everywhere, and is similar to what I'm doing with the generic geometry type system in my library. Luke

Simonson, Lucanus J wrote:

Bruno wrote:

...

Indeed it's something that should be done, we had a quick discussion about that with Barend and decided to postpone that for now. This kind of subject is never easy. For example, I'd like to rename circle into hypersphere or nsphere, but then we'll also have to clarify the meaning of some algorithms. For example what we currently call the area of a circle is actually the volume of an hypershere<2>...

My preference would be to have a circle data type, a sphere data type and a hypersphere data type. Only the hypersphere would be parameterized on dimension. The user would pick the object that is named for what they want. If a user wants a circle they are going to be frustrated and confounded that we provide them nsphere<2>.

circle would just be a typedef for nsphere<2>. Still, that has issues, since you will only be able to use that in 2D space. So if you want to use a 'circle' in 3D, it won't work, since it needs to know in which plan the circle is. You might as well use CSG to intersect a sphere with a plan and construct a circle in 3D. By the way, I think the geometry library will need to consider CSG at a point or another so that intersection can always return something, at least combinations with spheres. Unless it is chosen to approximate spheres to polyhedrons to avoid the issue.

How can you model a rectangle that is not axis aligned and is still a rectangle?

I think the minimal representation just needs three points (while an aligned one requires only two), whatever the dimension, but I'm not sure about this. Of course, you can choose to represent it as a quadrilateral, i.e. a polygon with 4 points and no holes. I fail to see how it not being axis-aligned makes its existence as a rectangle more difficult. Considering 2D, take for example that non axis-aligned rectangle (actually square because I'm lazy) defined by those points, in direct order: (0, 1), (1, 0), (0, -1), (-1, 0) Yes, that can model a rectangle that's not axis aligned, even if the coordinate systems is made with integers.

The only way I can thing of is to model it is width, height, center point and rotation. If you snap its corners to the coordinate grid it will, in general, no longer be a rectangle.

I don't really see how it stops being a rectangle once rotated. Are you thinking in terms of pixels maybe, because straight lines that are not axis-aligned aren't straight when rendered as pixels?

Similar to the circle, I like rectangle, prism and hyperprism for 2D, 3D and ND rectangles. While non-axis-aligned rectangles could be useful, you have the option of modeling them as polygons with 4 vertices that are as nearly rectangular as you can make them within your coordinate system.

Rectangles have some nice proprieties regular quadrilateral do not have, so you might be able to use more efficient algorithms by making use of a rectangle (meta-)type (or concept) instead of a quadrilateral one.

Michael Marcin wrote:

Do you mean typedef in the literal sense?

Yes and no. The circle concept could be equal to the 2-nsphere concept, if you prefer. And the circle class (or template class) that implements this concept could be a real typedef of the nsphere template class.

Template typedefs aren't available until C++0x.

Yet you can do it in C++03, it's just slightly more verbose (::type). I don't see how template typedefs are related to concepts though.

Is this library still targeting the current standard?

I suppose so. It's true using C++0x would make it cleaner, but you can have the same features in C++03 for not much more verbosity.

Phil wrote:

Personally I prefer to call a spade a spade, a circle a circle and a rectangle a rectangle. "Two dimensional hypershepere"??? "Two dimensional paralleloflange"??? "Quadropoint cultivational utensil"!!!

Whats wrong with typedef digging_object<long_handled, wood_shaft, foot_assist, pointed_end> spade; ? Sorry Steven, I couldn't resist ;) Luke

Mathias Gaunard wrote:

Apart from that, allowing multiple coordinate systems within a given space would be fairly useless, and in Euclidian space the Cartesian system is the best choice since it was made for algebra and computations. For example, points expressed in polar coordinates should automatically perform conversions from/to cartesian when using get.

I disagree with this. The best choice of coordinate system for internal representation is problem dependent. Imposing transformations to and from Cartesian internally can cause both numeric stability problems and performance problems. Every coordinate system was created for algebra, calculus and computations. Making the best choice of coordinates for a representation can speed calculations by orders of magnitude, or even in some cases make problems solvable that weren't in Cartesian coordinates (because of a lack of an appropriate representation for the solution). This is a huge deal in analytic calculations (See, for example the books by Morse and Feshback, or the early editions of Arfken where chapters are dedicated to finding and operating in the best coordinate system for the job. Also see facilities for this work in Maple or Mathematica.), but it is a sometimes under appreciated aspect of numeric calculation, as well. This is a conceptually similar problem to the discussions about unit systems from a couple of years ago. High accuracy and high performance calculations can't afford to have hidden conversion costs in computation and precision imposed on them. A well designed library should make such hidden conversions possible for those who want them, while also not making them mandatory for those who can't afford to pay for them. John Phillips

Mathias Gaunard wrote:

John Phillips wrote:

...

Every coordinate system was created for algebra, calculus and computations.

Yet when dealing with vectors and matrices, I'm pretty sure the cartesian coordinate system is what everyone uses.

It is usually the first thing to try, but it is not always the best idea. I won't claim that a large fraction of the computing community is conscious of the other possibilities, or of when they are useful. However, much of the scientific computing part of the community is familiar with multiple systems and reasons to use them.

...

Making the best choice of coordinates for a representation can speed calculations by orders of magnitude, or even in some cases make problems solvable that weren't in Cartesian coordinates (because of a lack of an appropriate representation for the solution). This is a huge deal in analytic calculations (See, for example the books by Morse and Feshback, or the early editions of Arfken where chapters are dedicated to finding and operating in the best coordinate system for the job. Also see facilities for this work in Maple or Mathematica.), but it is a sometimes under appreciated aspect of numeric calculation, as well.

I've heard that various analytical problems are much more easily represented in other coordinates. But is really the kind of problems you want a boost geometry library to deal with?

Apart from that, I've seen people use homogeneous coordinates (which are very similar to cartesian) to avoid doing some divisions that might lead to numerical instability. That should be an easy addition.

I would like a boost library that serves as large a segment of the problem space as is reasonably possible. I haven't tried to design one, so I don't know if the extension somehow collides with the abstractions and implementation details needed for a clearly designed Cartesian library. I don't think so, with the small amount of thought I've put into it so far, but that isn't the same as a real study of the problem. Some problems are simplified incredibly by a good choice of coordinate system. As a trivial example, consider a circle in two dimensional coordinates. There is no properly defined function in Cartesian coordinates that produces a circle. (There is a relation, but a circle is not single valued in Cartesian coordinates, and so can't be written as a function.) The function in plane polar coordinates is r = c (The radius equals a constant), when talking about a circle of radius c centered at the origin. This is a function, and a particularly simple one. Physics problems in orbital mechanics and the quantum mechanics of electron orbitals, and mathematical problems of interpolation using radial basis functions are less trivial examples. In each case, the selection of a good coordinate system for the problem allows for a clear separation of concerns in the calculation. Analytically this can make it possible to apply techniques to the isolated components of the problem. Computationally, this reduces the multi-dimensional problem into a set of independent or nearly independent one dimensional problems. Since the computational effort of many common numerical techniques scales exponentially in the number of dimensions, this separation moves us out of exponential scaling and into linear scaling in the number of dimensions. I've never met a computer programmer who wouldn't approve of that trade.

...

This is a conceptually similar problem to the discussions about unit systems from a couple of years ago. High accuracy and high performance calculations can't afford to have hidden conversion costs in computation and precision imposed on them. A well designed library should make such hidden conversions possible for those who want them, while also not making them mandatory for those who can't afford to pay for them.

I thought it was better for simplicity, dealing with different geometries seems hard enough, but a prerequisite of the library since it aims at dealing with GIS; and hopefully, there won't be too many geometries to handle. On the other hand, different representations of points in the same space are all equivalent, and one representation is enough to describe everything we want to deal with, even if not in the most precise way.

The problem is that it seems hard for algorithms to consider all the possibilities in existence for the points they're dealing with to be represented and select the right way to compute something accordingly. Latitude/Longitude seems close enough to cartesian so I supposed that a lot of code can be shared.

All of the typical algebraic, differential and integral structures have well known translations into this world of different translations. For ortho-normal systems the key is the Jacobian. Even systems that are not ortho-normal can be approached with tensor methods. For some very common structures, like the dot product or the cross product, the algorithm doesn't even change. For differential and integral structures, the change is typically covered by throwing in a multiplicative factor. Details of this are available in most vector calculus books, or in the above mentioned books on mathematical physics from my previous note.

Maybe only doing vector operations and generalizing to any orthogonal coordinate system is possible without much hassle by working with normalized basis vectors, but wouldn't that be quite similar to converting to cartesian everytime?

No, it wouldn't. The Cartesian unit vectors are constants, and the Cartesian area and volume elements do not depend on current coordinate values. These facts are generally not true for non-Cartesian systems. (In plane polar coordinates, the same change of angle and the same difference in inner radius compared to outer radius sweeps out different areas depending on the actual value of the inner radius. As the inner radius gets bigger, the area swept out increases.) For problems where the natural description of the paths and shapes looks like rectangles, Cartesian coordinates provide a simple and accurate basis for calculation. However, if the natural paths and shapes don't look like parts of rectangles, Cartesian coordinates provide a bad description of the system and there is a need for substantial extra calculation to process the system in those coordinates. The answer is to look for a different system that matches the natural paths and shapes. (Better put, it shares the same symmetries as the problem.) These symmetries are automatically preserved just by writing the problem in the new coordinate system. There is no need to put computational time and effort into calculating them, or into sanity check functions that check for them and correct the results if they stray away from the required symmetries. In a best choice coordinate system, one characteristic change in the system can be well represented as one small change in a single system variable. (If it is a geometric system we are describing, the variables are coordinates.) This is good for the exact same reason that isolating effects in different modules in a computer program is good.

I honestly have zero knowledge of non-cartesian coordinates, so I have no knowledge of the vector theory behind it (and I have no knowledge of how vectors apply to alternative geometries either). Maybe that could somehow also be applied to help dealing with other geometries that are locally euclidian?

If there are ways to write algorithms generically in a coordinate-system-independent way, and maybe even geometry-independent way, then it would be great, for sure. But is that really possible?

This is the real question. My first guess is that it is possible for many algorithms. However, this is only an educated guess, and not the result of time spent producing real and practical implementations, so it could be wrong. In the cases where it isn't possible, we should try to do the same thing all generic programming does. Define a minimal set of concepts that must be modeled by the coordinate system for the algorithm to work correctly and embed such concepts and concept checks in the library. John

Mathias wrote:

I don't even know how one could define a polygon as a function in any coordinate system, by the way.

You define it as a piecewise linear logic inequality. A planar polygon can be thought of as a function on x and y such that for any x and y that lies inside the polygon the function returns 1 and for any x and y that lies outside the polygon the function returns 0. The actual mathematic expression for doing that isn't really interesting, but thinking of a polygon as a mathematical function is. I take the partial derivative of a polygon in x and y to reduce it from something that is conceptually two dimensional to a series of zero dimensional quantities. These quantities are impulse functions with either positive or negative sign at a given point with a given direction. This means I can model them as positive or negative unit vectors. When I scan those quantities, I am taking the integral with respect to x then y and I get back my original polygon. When I scan many polygons this way I can do boolean logic operations on the 2D logic functions that are the input polygons. Because of the way I model polygons the algorithm employs arithmetic rather than flow control to do the core computational effort. Obviously, this approach extends easily into higher dimensional space at the conceptual level, but in practice, it is very hard to write an efficient n-1 dimensional sweep-plane algorithm due to a lack of suitable data structures, but it is easy to write a scan-line that handles the 2D case in near linear time using, for example, the stl red-black tree. Figures that have curves are not piecewise linear and harder to take the derivative/integrate because there is no order reduction taking place. In a different coordinate system the curves might be analogous to straight lines, allowing for order reduction, but the more common approach in computational geometry is to simply model the curves as their piecewise linear approximation and solve the problem in the Cartesian coordinate system. While non-Euclidean geometries are undeniably valuable, I think that it is not worthwhile to try to generalize the geometry away from the geometry library. Refactoring algorithms that are similar in different coordinate systems to be shared in their implementation is probably harder than duplicating that effort in separate libraries that target separate geometries. I think it is preferable to design the library to interoperate well with other geometries, rather than try to encompass them. Luke

Phil wrote:

- I note that you're supporting homogeneous coordinates in arbitrary dimensions. My current requirements are for at least partially heterogeneous coordinates in 2 or 3 dimensions. Of course you can't support all permutations and combinations of requirements, but you should perhaps qualify some names to indicate their limited scope.

Can your coordinates cast to and from one coordinate type at the interface between your data structure and the algorithms provided by a generic library?

- I have always used a very simple point type: struct point { T x; T y; };. It seems that I could continue to use "user" code in this style with your library; however, as a potential library author, I need to replace "p.x" with something like "access<P>::get<0>(p)". I'm not wildly enthusiastic about this verbosity. I accept that I'm outvoted on this, but I'm going to continue to remind you (all) that I'm unhappy about it.

You only need to use the accessor function when you don't know the data type of the object. If you are writing code within the library you can use std::pair<coordinate_type> for point and throw it in a vector and sort it and just get on with your life. The accessors are just for the interfaces.

- I have previously noted that you can advertise your distance function as returning a "type convertible to T" in order to avoid unnecessary sqrt calculations while, for example, scanning a list of points to find the one furthest from p.

I have a euclidean_distance_squared function on geometries that does what you say. Distance calls distance squared and returns its square root.

- I look forward to seeing concepts for shapes. - Is there a strategy for "within" that can be used when the shapes being tested store their bounding boxes?

I have a winding direction accessor function on my generic interface to unknown polygon type T. If the polygon knows its winding direction in constant time it can return it, otherwise it returns unknown and I compute the winding direction. A similar thing could be done for bounding box, where the cached bounding box can be returned in constant time. Because these things make the interface bigger, I'd like to come up with a way to make them syntactically optional. Any competent implementation of within on polygon would do bounding box comparison first to rule out simple cases. Getting the bounding box would go to the bounding box accessor first and could get your cached bounding box to speed up the check. Luke