Most Useful Math

[jemidiah]

I'm interested in writing another tutorial on some math topic. I'm curious, is there anything people would be particularly interested in? I'm especially interested in something that doesn't already exist--say a programmer's approach to some math problem. There are a zillion topics that would fit the bill. I suppose I could list a few:

* Counting things, and coding it
* Calculating sin(x) from scratch
* A survey of calculus
* Numerical integration
* Basic probability
* Polynomials and You
* Basic real analysis
* Trig functions in code
* Modular arithmetic
* Theoretical CS Logic

Thanks for any responses .

[wossname]

Personally, I'd be very interested in some of the following, in descending order of interest (some of these may be more logical or algorithmic than mathematical so I'm just brainstorming really...)

* 4 dimensional ballistic trajectory plotting (x,y,z and time). Basically this would help me out with shooting a cannonball at a moving target that may or may not be at the same altitide as the cannon.

* Plotting (projecting) 3D points onto a 2D plane (a computer screen). The mathematics of this I have always found difficult to understand.

* Probability - the example that you NEVER see on the internet... I have n independent events (where n>3), each of which has a known probability. What is the probability of ANY of those events occurring? What is the probability of ALL of those events occurring? The only example I saw of this was very badly explained and was filled with dreadful errors in the calculations.



I am pretty thick-skulled though so go easy on the jargon eh? I almost failed A-Level maths at college, I was too busy writing programs on my TI-82 at the back of the class to pay any attention

[jemidiah]

The ballistic trajectory one sounds good. Mapping 3D points to a 2D plane would unfortunately get into linear algebra, and use a lot of pictures I'm not willing to draw outside my head :P. Probability could work too, though the independent probability you asked for was a question in this forum a while back.

Any other thoughts?

[wossname]

* Monte-Carlo methods or the math behind fractalline patterns are always popular subjects.

* How about computing the surface-areas of cross sections of irregular solids!

* Fast ways to compute the area of irregular (2d) polygons.

[jemidiah]

The 2D polygon area came up a while back. So long as you have an ordered sequence of points, it turns out there's an easy formula for it. A "real" discussion of Monte-Carlo methods would require some somewhat in-depth stats. Computing the surface areas of cross sections of irregular solids might be interesting--what format would the solid take? That is, how would one figure out the boundary? (I can think of many ways; I just want to know which you had in mind.)

[wossname]

While we're on the subject, the actual volume of a 3d solid would be useful too.

As for how you'd represent the solid - I guess the normal way... a list of all vertices and a list of faces that group the vertices together.

I think this would be a bit heavy on the algorithmic side, possibly less mathematical.

[jemidiah]

Analytic volume of a 3D solid given nothing but a graph of vertex connections... I've never seen a way to find a volume given that information, but I'm sure one exists. Interesting; I might try to find a way to do that.

[wossname]

Hehe I never said it was easy One exists, in fact many exist, some more efficient than others.

[krtxmrtz]

I'm interested in polygon interpolation.

Assume you have the 2 cross sections of a 3D object corresponding to the intercepts by 2 parallel planes at different heights. These cross sections could be closed curves but assume you have a representation of them consisting in 2 sets of points, in other words, you have 2 polygons (that may have a different number of vertices).

The problem is how to find the interpolated cross section (polygon) at any intermediate height.

I have been working on this for a while and have come up with a solution, but my approach is not dealing with degenerate cases (i.e., cross sections having multiple closed contours, etc). Also, if the cross sections are overlapping each other I don't know how to proceed.

[jemidiah]

Can we assume the points in each cross section make possibly several convex hulls? While such an assumption would make things easier, in general I don't think the problem is determined enough. Take the points in the following configuration:

    *   *
*           *
    *   *

This configuration could be generated by two triangular shafts placed right next to each other, or it could be made by a single hexagonal shaft, but how are we to know which is the correct interpretation?

I'll have to think about it more, though. It's an interesting question, and the above ambiguity might not make a whole lot of difference given that you only want interpolated cross sections, presumably given by sets of points.

[jemidiah]

Yup, I believe this problem has too little information to be exactly or approximately solvable in its current state. Take the following two slices:

Top slice:

*  *  *  *
*  *  *  *

Bottom slice:

*  *  *  *
*  *  *  *

Suppose that these slices are generated by two square shafts. Suppose the leftmost 4 points of the top slice are due to shaft 1, leaving the rightmost 4 points of the top slice to be due to shaft 2. Which shaft generates each square of 4 points in the bottom slice?

Case 1: If shaft 1 generates the leftmost 4 points of the bottom slice, shaft 2 must generate the rightmost 4 points. This would cause the shaft axes to be parallel to each other as well as parallel to the normal of the planes--the shaft axis would be straight up and down with nothing strange going on.

Case 2: However, if instead shaft 2 generates the leftmost 4 points of the bottom slice, shaft 1 must generate the rightmost 4 points. This means that the shafts must cross through each other. Their shaft axes would no longer be parallel. Viewed from another angle, the shafts would form an "X".

This is important because it shows that the exact same data may result in different intermediate cross sections. In Case 1, the cross section midway between the top and bottom slices would be simple:

Case 1 middle slice:

*  *  *  *
*  *  *  *

However, in Case 2, the cross section would be even simpler, cutting through the middle of the "X":

    *  *
    *  *

So, no matter what interpolation algorithm you use given the information above, you'll get this particular interpolation wrong sometimes. Even adding requirements that individual polygons form convex hulls and that we somehow know which line segments correspond to edges on the 3D surface, the same result follows.

Is there something I'm missing?

[krtxmrtz]

All right, I realize that more clarification is necessary.

First of all, let me explain a possible practical side of this problem, actually the one I'm interested in.

In medical imaging, you may have, say, 2 axial CT slices with anatomical information at 2 different levels (in the head-feet direction) -if you're not familiar with CT scan images I suggest you make a search with Google images. Then you draw on each slice to outline a specific anatomical structure, for example the liver. Usually the 2 levels are very close to each other (a few mm) but in some cases the separation is higher and, for some purposes I need not go into now, it may be necessary to figure out the interpolated liver contour in between.

What all this means is:

1. The contours may be concave.
2. There may be multiple contours belonging to the same structure (for example, some pelvic bones).
3. In actual practice the polygons from different levels may not be too different.
4. But even if they are very different, the nature of the problem doesn't change. If we forget about this medical application and consider the general problem it's obvious that there are infinite possible solutions. However, the algorithm must be such as to yield a "reasonable" and possibly the simplest solution, i.e., assuming that the 3D structure is as "orthodox" so to speak as possible, not subject to strange kinks in between the given CT slices.

[krtxmrtz]

I have attached my polygon interpolation project. Of course I don't pretend you work on it, even with comments is messy so I'd rather explain the algorithm I've used, which I will do eventually as I haven't been working on it for quite some time and have forgotten most of the details.

The idea is you can play around with the interface and get a feeling for the problem.

[krtxmrtz]

...Fast ways to compute the area of irregular (2d) polygons.

Yes, there's an easy formula. Based on the attached image, all you have to do is compute A₁ and A₂, and each one of these can be split in the sum of the areas of the various trapezoids separated by the dashed lines. And each trapezoid in its turn can be split in a rectangle and a triangle.

Then, for example, for the trapezoid bound by segment 1-2 the area would be:

A_1-2 = (x₂ - x₁)y₁ + 1/2 (x₂ - x₁)(y₂ - y₁) = 1/2 (x₂ - x₁)(y₂ + y₁)

So, in general, the total area A will be the sum of terms like the above for all trapezoids:

A = SUM(for i=1 to N) [1/2 (x_i+1 - x_i)(y_i+1 + y_i)]

Of course, for this to work you have to add one extra point N+1 duplicating the first one.

Notice that when x_i > x_i+1 the addend may be negative thus contributing to substract the A₂ part.

[zaza]

Incidentally, just to add 2c, there's a link in my sig to some coordinates and vectors stuff that I did a while ago. There's some 3D solids generation, converting 3D points into differnet representations and plotting to screen etc.
VB.Net, I'm afraid, but still.....

[jemidiah]

I looked at your code, which seems pretty reasonable. Of course I haven't checked every piece of reasoning, but only got a general idea of how everything works. To outline the algorithm...

0. Start with two slices. Each slice ("Polygon" in the code) contains some number of non-self-intersecting polygons ("Contours" in the code). I'll use the code's wording instead of math's wording when they conflict. None of the contours intersect each other inside a given polygon, either.

1. Invert the orientation of each contour in each polygon, if needed, so that everything has a consistent clockwise orientation.

2. Renumber the vertices in each contour in each polygon in such a way that first vertices of different contours are in approximately the "same" location on their contour, relatively speaking. The algorithm used takes the x-midpoint of the contour (the point on the x axis halfway between the farthest left and farthest right points on that contour). It intersects the vertical line passing through that midpoint with the edges of the contour and adds the topmost point of intersection to that contour. Vertices are renumbered so that this intersection is point number 1 and the clockwise ordering is kept.

3. The contours are resampled in such a way that the number of vertices on every contour is the same. This is done by sampling at a high enough rate to guarantee that the shortest edge will still have at least one point. (Actually, I think there's a bug that's causing the sampling rate to be even higher than it should be, though aside from speed this isn't really a problem. The length calculation was too high when I looked at a couple of triangles. I also think there's another bug causing some points to be taken several times, which sometimes causes distorted and eventually self-intersecting interpolations.)

4. Every contour in the first polygon is interpolated with every contour in the second polygon. That is, a new contour made out of somehow combining the two contours is generated. For n contours in the first polygon and m in the second polygon, this generates n*m new contours--which is of course very wrong. The actual interpolation used is quite simple. Since the contours all have the same number of points and each point is at approximately the same "place" on every contour from renumbering, a simple linear interpolation is used on pairs of same-numbered vertices from each contour in the interpolation. This generates a very smooth transition between contours.


First, aside from bugs, the interpolation part of this algorithm seems pretty good to me. There's some bias involved in choosing the number 1 vertex as described above, but not a whole lot, and something had to be chosen, somehow, regardless. Perhaps there are some annoying edge cases where this would become a problem, but improving the number 1 choice algorithm doesn't seem like it would net a large gain. So, steps 1-3 seem fine.

Step 4 is the real culprit, especially in your difficulty with handling multiple contours. I would suggest a different approach in determining which pairs of contours should be interpolated.

First, calculate the center of mass of every contour in both polygons. Formulas for discrete center of mass calculations can be found online. You would assume that the polygon bounds a mass of uniform density and find the center of that mass. If you wanted to approximate this instead, you could more easily take the average of the vectors describing the edges of the contour. Physically, this would assume that all of your mass is uniformly distributed at only the vertices of a contour.

Second, determine which contours in the first slice should morph into which contours in the second slice. The most reasonable thing I can think of here is for each contour in the first slice, find the closest center of mass for a contour in the second slice. If Contour A in the first slice has its center of mass located 10, 20, and 30 units away from Contours 1, 2, and 3 in the second slice, then only interpolate between Contour A and Contour 1. Don't include the interpolations between Contour A and Contours 2 or 3.

A slight problem crops up when contours aren't uniquely paired off. For instance, you might have Contours A and B both morphing into Contour 1 if Contour B is close to Contour 1 and far from Contours 2 and 3. Only a slight modification of this algorithm should be needed here. We generate Contours A1 and B1 using the existing linear interpolation scheme. We should then simply polygon union A1 and B1. If A1 and B1 don't intersect each other, this will not affect them. If A1 and B1 do intersect each other, we would want to be sure we cut out the edges that accidentally stay internal to the polygon.

Another slight problem happens when a contour in one slice has no corresponding contour in another slice. It seems reasonable to me to simply shrink the given contour down to nothing in such a case. Shrinking would be done by scaling using the center of mass as the center of scaling. However, if slices are taken sufficiently close to each other, this shouldn't happen.


What did you mean when you said, "Also, if the cross sections are overlapping each other I don't know how to proceed."?

[krtxmrtz]

I looked at your code, which seems pretty reasonable. Of course I haven't checked every piece of reasoning, but only got a general idea of how everything works. To outline the algorithm...

Congratulations! Deciphering someone else's code is one of the greatest achievements of mankind, the more so because I consider myself a sloppy coder.

Steps 1-3 look right but I don't quite like the algorithm after all. That's because I've found some examples of misbehavior, for example you can see in the attached image that the resulting contour is improperly clipping both input contours, so something is wrong or maybe the sampling rate should be higher.

As for the 4th step, I must admit it's difficult as the solution is somewhat ambiguous. What you suggest about using the center of mass when there are multiple contours in both input sections seems reasonable. Anyway, I must think about it.

As for your question about overlapping: if the first cross section has one closed contour, say a circle, and the second has 2 circles then as you interpolate ever closer to the first cross section, the output will be 2 circles ever closer to one another until eventually the will overlap. At this stage, a polygon clipping algorithm is required to figure out the union. For this reason I mentioned I didn't know how to proceed (general polygon clipping is no joke). But I've since remembered I once downloaded a dll from this site and I performed some successful tests with it, though then I left it aside in my code bank and it didn't strike me I could use it in this problem. But then, the polygon structures used in this dll may not be compatible with mine so that complicates the whole thing even further.

[jemidiah]

You'll almost certainly have to use polygon clipping eventually, I'd say. Imagining the bones of the forearm, a slice taken when they're separate and a slice taken towards the base when they're not (they aren't still separate, right? I never really bothered with biology/anatomy) would require a union operation to interpolate properly.

I thought a bit more about my center of mass idea. My original idea was as follows: for each contour in slice 1, find out which contour in slice 2 it should evolve into. Evolve it. However, this leads to asymmetry in that if we switch the first and second slices we may not get the same interpolation. (Take two contours in slice 1 going to one contour in slice 2, and then reverse the situation to see the asymmetry in a simple case.) This doesn't seem that bad to resolve, though. The following should work:

1. For every contour in slice 1, calculate the contour in slice 2 that it would "like" to evolve into, using the minimum center of mass distance (or a random number generator, or whatever you feel is appropriate).
2. For every contour in slice 2, calculate the contour in slice 1 that it would "like" to evolve into.
3. Generate an interpolated contour for every pair of contours where one would like to evolve into the other, even if the feeling isn't mutual, i.e. even if one contour is closer to some other contour in the other slice but another contour in the other slice is closest to it, as in the 1-2 situation above. As an optimization, if two contours want to morph into each other, generate only one interpolation.
4. Union everything, always. Who knows when the contours will overlap, so at each step of the interpolation just cut out internal edges.

If the trouble in the image you attached is that the edges of the interpolation are sometimes not sharp points, I noticed the same thing. Something even worse happened if I tried to have one triangle morph into another, as the attached picture shows. The second picture shows the interpolation lying exactly on top of the second slice, which shows that most of the points are accidentally "doubled up"--several number labels are overlapping each other. As expected this screws up the interpolation routine.

[krtxmrtz]

Given the attached figure, there are indeed many possibilities. For example, contours A and B could evolve into B' and at the same time A could also evolve into A' as if this were some sort of N-shaped 3D structure. Or A could evolve into B' and B into A', or... or...

But a reasonable guess -though not necessarily corresponding to the real case- is that A evolves into A' and B into B'. In this event you have to refine the center of mass rule which would pair A with B'. Offhand, you could adopt this extra rule (distances between CMs are meant):

If Abs(AA' - BB') < Abs(AB' - A'B) then pair AA' and BB'
If Abs(AA' - BB') > Abs(AB' - A'B) then pair BA' and AB'

Anyway, I don't know if this fails for some other example and how to extend it to cases with more contours complicating it further.
------
I'm not sure I see what you mean when you mention the asymmetry of the problem when switching slices. Using my (otherwise faulty) algorithm I think the solution should be the same. Maybe you can provide a clearer example.

In the meantime I think I shall try to take time to review the algorithm in its basic part dealing with single contours in each slice..

[jemidiah]

The algorithm I listed about just taking the nearest center of mass definitely fails when the slices are far apart in the z direction, but I think that's a fault of the setup and not the algorithm. You should just take closer slices if that's the problem. Eventually the center of masses will get close to their original position.

I think that if you try to cover all possible cases where the slices can be taken to be relatively far away, you could come up with all sorts of unintuitive edge cases, or edge cases where a particular algorithm won't "guess" in the same way that you might "guess". At some point, you'd probably just go for machine learning to see what is and is not a good evolution strategy, but that seems... silly.

However, if you just wanted to extend things so that, given equal numbers of contours in each slice, each contour is matched with exactly one contour in the other slice, you could perhaps make that work. This is a very special case and seems a bit hacky to me considering the above, but ah well. In this case, you could minimize the total (variance in?) "travel distance" between all pairs of contours, that is, the sum of all distances of the form XY' where contour X is paired with contour Y'. There are many edge cases where this won't work well, though, so... meh.


I'll try to come up with a clearer example of asymmetry. The basic trouble is that if you implement my "for contour in polygon1: interpolate(contour, contour.nearestNeighbor(polygon2), currentZ)" algorithm, you can often get two different results if you swap polygon1 and polygon2. This is despite the fact that the interpolation should physically be invariant when you swap the two polygons and adjust the currentZ appropriately.

[krtxmrtz]

The problem with my algorithm is, to begin with, it doesn't account well for concave cases. Maybe I'll have to start a new one from scratch.

In the mean time, about polygon boolean operations, I coded my own algorithm some time ago but it didn't cope with degenerate cases, for example, a vertex from a polygon lying on the side of another. Then I found this paper which solves the problem in the general case, but after I read a few pages finally I got lost and didn't understand the rest.

[jemidiah]

Implementing that paper, at a glance, looks like a mess to code and debug. Using a library is probably the way to go, even if it's pretty inconvenient in other places.

The general algorithm I outlined above should deal pretty well with concave cases. Is there a specific part of your implementation that doesn't deal with them well? (I know there are some bugs as pointed out above; perhaps that's the problem.)

[krtxmrtz]

Implementing that paper, at a glance, looks like a mess to code and debug. Using a library is probably the way to go, even if it's pretty inconvenient in other places.

The thing is I have a full implementation of a clipping algorithm that I worked out on my own after I had read a different paper dealing with some basic algorithm. But the part yet to be implemented is precisely the extension to degenerate cases explained in this paper. I'd love to understand it so I could take the challenge of finishing my clipping algorithm.

The general algorithm I outlined above should deal pretty well with concave cases. Is there a specific part of your implementation that doesn't deal with them well? (I know there are some bugs as pointed out above; perhaps that's the problem.)

I'm not sure, sometimes it deals well with concave cases, sometimes it doesn't. As you point out it may be that some of the bugs that crept in are making the algorithm fail in some special examples. As I said, I must find time to calmly sit down and try take care of the bugs.

[tiguy]

I would be interested in Aspheric surfaces. maybe Sag tables to go with it.

[jemidiah]

I'm sorry, I know almost nothing about such surfaces. The conic section lenses and perturbations are geometric and algebraic and therefore general, and the basic idea is pretty simple. I'll give a very brief explanation (which you can probably find with an internet search).

Physically creating a lens can be difficult. It's relatively easy to create a sphere and to cut that sphere by "intersecting" it with a plane. This creates a standard lens, which focuses idealized light sources. However, in the real world, light sources are non-ideal and other shapes may work better, depending on the specifics of your application. A lens which isn't cut from a sphere is called an aspheric lens or surface. One important property that aspheric lenses enjoy is "rotational symmetry": rotating the lens about its optic axis results in no apparent change. For instance, the function e^-(x^2 + y^2) is rotationally symmetric. This property corresponds to a human tilting their head without (appreciably) distoring the image their eyes see.

The "sag" of an aspheric surface is a way to completely describe the lens' shape. This article has a decent picture. The "Z" is the "sag" (the arrow below the text "'sag' height" seems useless). The equation for the sag Z as a distance of "r" as defined in the image in the linked article can be quite arbitrary, though it must be an even function so the top half of the lens and the bottom half are mirror images. Since the curve Z(r) should really be infinitely differentiable (actually, I meant analytic), one can describe it using a Maclaurin series. This leads to a very general definition of the sag Z(r) as a high-order polynomial with some fixed constants determining exactly how it curves. Evenness forces us to drop each term with odd degree in the Maclaurin series, and a constant term is pointless so it's generally set to 0.

For convenience, you can also add another non-power-series term, though if you wanted, this could be eliminated by changing infinitely many of the polynomial constants. This term gives the lens the character of a conic section--a circle, which gives the basic spherical lens, assuming all polynomial correction terms are 0; an ellipse; a parabola; or a hyperbola. The use of the different correction terms and/or different conic sections is entirely application-dependent. For application-specific details, I have to say I know nothing.

For further reference, the Wikipedia article lists a number of articles.