Comp.graphics.algorithms is an unmoderated newsgroup intended as a forum for the discussion of the algorithms used in the process of generating computer graphics. These algorithms may be recently proposed in published journals or papers, old or previously known algorithms, or hacks used incidental to the process of computer graphics. The scope of these algorithms may range from an efficient way to multiply matrices, all the way to a global illumination method incorporating ray tracing, radiosity, infinite spectrum modeling, and perhaps even mirrored balls and lime jello.

It is hoped that this group will serve as a forum for programmers and researchers to exchange ideas and ask questions on recent papers or current research related to computer graphics.

comp.graphics.algorithms is not:

• for requests for gifs, or other pictures
• for requests for image translator software (i.e. gif <--> jpg)

In 2-D, the 2x2 matrix is very simple. If you want to rotate a column vector v by t degrees using matrix M, use

M = {{cos t, -sin t}, {sin t, cos t}} in M*v.

If you have a row vector, use the transpose of M (turn rows into columns and vice versa). If you want to combine rotations, in 2-D you can just add their angles, but in higher dimensions you must multiply their matrices.

Let the point be C (XC,YC) and the line be AB (XA,YA) to (XB,YB). The length of the line segment AB is L:

              ___________________
             |        2         2
        L = \| (XB-XA) + (YB-YA)

and

            (YA-YC)(YA-YB)-(XA-XC)(XB-XA)
        r = -----------------------------
                        L**2

            (YA-YC)(XB-XA)-(XA-XC)(YB-YA)
        s = -----------------------------
                        L**2

Let I be the point of perpendicular projection of C onto AB, the

        XI=XA+r(XB-XA)
        YI=YA+r(YB-YA)

Distance from A to I = r*L
Distance from C to I = s*L

If r < 0 I is on backward extension of AB
If r>1 I is on ahead extension of AB
If 0<=r<=1 I is on AB

If s < 0 C is left of AB (you can just check the numerator)
If s>0 C is right of AB
If s=0 C is on AB

Well, there are two ways, the right way and the fast way. Representing the ellipse in implicit form:

   2              2 
A x  + B x y + C y  + D x + E y + F = 0

The fast way is to calculate

         2              2
      A x  + B x y + C y  + D x + E y + F
-------------------------------------------------
 sqrt((2 A x + B y + D)^2 + (B x + 2 C y + E)^2)

The right way is to rotate and translate the coordinate system so that the ellipse is in canconical form:

   2      2 
a x  + b y  = 1

and solve for lambda in

                2      2
             a x  + b y  = 0
      x + lambda (2 a x) = x0
      y + lambda (2 b y) = y0

This gives a quadric equation in lambda:

    2                2      2                2                     2                2
a x0 (1 + 2 b lambda) + b y0 (1 + 2 a lambda)   =  (1 + 2 a lambda) (1 + 2 b lambda)

You can't. The only case where this is possible is when the bezier can be represented by a straight line. And then the parallel 'bezier' can also be represented by a straight line.

A Bezier curve is a parametric polynomial function from the interval [0..1] to a space, usually 2-D or 3-D. Common Bezier curves use cubic polynomials, so have the form

f(t) = a3 t^3 + a2 t^2 + a1 t + a0,

where the coefficients are points in 3-D. Blossoming converts this polynomial to a more helpful form. Let s = 1-t, and think of tri-linear interpolation:

        F([s0,t0],[s1,t1],[s2,t2]) =
            s0(s1(s2 c30 + t2 c21) + t1(s2 c21 + t2 c12)) +
            t0(s1(s2 c21 + t2 c12) + t1(s2 c12 + t2 c03))
            =
            c30(s0 s1 s2) +
            c21(s0 s1 t2 + s0 t1 s2 + t0 s1 s2) +
            c12(s0 t1 t2 + t0 s1 t2 + t0 t1 s2) +
            c03(t0 t1 t2).

The data points c30, c21, c12, and c03 have been used in such a way as to make the resulting function give the same value if any two arguments, say [s0,t0] and [s2,t2], are swapped. When [1-t,t] is used for all three arguments, the result is the cubic Bezier curve with those data points controlling it:

          f(t) = F([1-t,t],[1-t,t],[1-t,t])
               =  (1-t)^3 c30 +
                 3(1-t)^2 t c21 +
                 3(1-t) t^2 c12 +
                 t^3 c03.

Notice that

           F([1,0],[1,0],[1,0]) = c30,
           F([1,0],[1,0],[0,1]) = c21,
           F([1,0],[0,1],[0,1]) = c12, _
           F([0,1],[0,1],[0,1]) = c03.

In other words, cij is obtained by giving F argument t's i of which are 0 and j of which are 1. Since F is a linear polynomial in each argument, we can find f(t) using a series of simple steps. Begin with

        f000 = c30, f001 = c21, f011 = c12, f111 = c03.

Then compute in succession:

        f00t = s f000 + t f001, f01t = s f001 + t f011,
        f11t = s f011 + t f111,
        f0tt = s f00t + t f01t, f1tt = s f01t + t f11t,
        fttt = s f0tt + t f1tt.

This is the de Casteljau algorithm for computing f(t) = fttt.

It also has split the curve for the intervals [0..t] and [t..1]. The control points for the first interval are f000, f00t, f0tt, fttt, while those for the second interval are fttt, f1tt, f11t, f111.

If you evaluate 3 F([1-t,t],[1-t,t],[-1,1]) you will get the derivate of f at t. Similarly, 2*3 F([1-t,t],[-1,1],[-1,1]) gives the second derivative of f at t, and finally using 1*2*3 F([-1,1],[-1,1],[-1,1]) gives the third derivative.

This procedure is easily generalized to different degrees, triangular patches, and B-spline curves.

Interestingly enough, bezier curves can approximate a circle but not perfectly fit a circle.

Michael Goldapp, "Approximation of circular arcs by cubic polynomials" Computer Aided Geometric Design (#8 1991 pp.227-238)

Tor Dokken and Morten Daehlen, "Good Approximations of circles by curvature-continuous Bezier curves" Computer Aided Geometric Design (#7 1990 pp. 33-41).

This is the Sutherland-Hodgman algorithm, an old favorite that should be covered in any textbook. Any of the references listed in the FAQ should have it. According to Vatti (q.v.) "This [algorithm] produces degenerate edges in certain concave / self intersecting polygons that need to be removed as a special extension to the main algorithm" (though this is not a problem if all you are doing with the results is scan converting them.)

A simple method is to put the bitmap through the filter:

        -1    -1    -1
        -1     8    -1
        -1    -1    -1

This will highlight changes in contrast. Then any part of the picture where the absolute filtered value is higher than some threshold is an "edge".

We describe the shape and position of objects using numbers, usually (x,y,z) coordinates. For example, the corners of a cube might be {(0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1)}. A deep understanding of what we are saying with these numbers requires mathematical study, but I will try to keep this simple. At the least, we must understand that we have designated some point in space as the origin--coordinates (0,0,0)--and marked off lines in 3 mutually perpendicular directions using equally spaced units to give us (x,y,z) values. It might be helpful to know if we are using 1 to mean 1 foot, 1 meter, or 1 parsec; the numbers alone do not tell us.

A picture on a screen is two steps removed from the 3D world it depicts. First, it is a 2D projection; and second, it is a finite resolution approximation to the infinitely precise projection. I will ignore the approximation (sampling) for now. To know what the projection looks like, we need to know where our viewpoint is, and where the plane of the projection is, both in the 3D world. Think of it as looking out a window into a scene. As artists discovered some 500 years ago, each point in the world appears to be at a point on the window. If you move your head or look out a different window, everything changes. When the mathematicians understood what the artists were doing, they invented perspective geometry.

If your viewpoint is at the origin--(0,0,0)--and the window sits parallel to the x-y plane but at z=1, projection is no more than (x,y,z) in the world appears at (x/z,y/z,1) on the plane. Distant points will have large z values, causing them to shrink in the picture. That's perspective.

The trick is to take an arbitrary viewpoint and plane, and transform the world so we have the simple viewing situation. There are two steps: move the viewpoint to the origin, then move the viewplane to the z=1 plane. If the viewpoint is at (vx,vy,vz), transform every point by the translation (x,y,z) --> (x-vx,y-vy,z-vz). This includes the viewpoint and the viewplane. Now we need to rotate so that the z axis points straight at the viewplane, then scale so it is 1 unit away.

After all this, we may find ourselves looking out upside- down. It is traditional to specify some direction in the world or viewplane as "up", and rotate so the positive y axis points that way (as nearly as possible if it's a world vector). Finally, we have acted so far as if the window was the entire plane instead of a limited portal. A final shift and scale transforms coordinates in the plane to coordinates on the screen, so that a rectangular region of interest (our "window") in the plane fills a rectangular region of the screen (our "canvas" if you like).

I have left out details of how you define and perform the rotation of the viewplane, but I'm sure someone else will be happy to supply those if there is demand. It requires knowing how to describe a plane, and how to rotate vectors to line up. Neither is difficult, but this is already using a lot of net space. One further practical difficulty is the need to clip away parts of the world behind us, so -(x,y,z) doesn't pop up at (x/z,y/z,1). (Notice the mathematics of projection alone would allow that!) But all the viewing transformations can be done using translation, rotation, scale, and a final perspective divide. If a 4x4 homogeneous matrix is used, it can represent everything needed, which saves a lot of work.

Assuming you want to rotate vectors around the origin of your coordinate system. (If you want to rotate around some other point, subtract its coordinates from the point you are rotating, do the rotation, and then add back what you subtracted.) In 3-D, you need not only an angle, but also an axis. (In higher dimensions it gets much worse, very quickly.) Actually, you need 3 independent numbers, and these come in a variety of flavors. The flavor I recommend is unit quaternions.

Paul S. Heckbert, "Survey of Texture Mapping", IEEE Computer Graphics and Applications V6, #11, Nov. 1986, pp 56-67 revised from Graphics Interface '86 version

Eric A. Bier and Kenneth R. Sloan, Jr., "Two-Part Texture Mappings", IEEE Computer Graphics and Applications V6 #9, Sept. 1986, pp 40-53 (projection parameterizations)

Complete source code appears in Graphics Gems IV pp. 175-192. A fairly complete sketch of the code appeared in the original article, in Graphics Interface 92 Proceedings, available from Morgan Kaufmann.

If the plane is defined as:

a*x + b*y + c*z + d = 0

and the line is defined as:

x = x1 + (x2 - x1)*t = x1 + i*t
y = y1 + (y2 - y1)*t = y1 + j*t
z = z1 + (z2 - z1)*t = z1 + k*t

Then just substitute these into the plane equation. You end up with:

t = - (a*x1 + b*y1 + c*z1 + d)/(a*i + b*j + c*k)

If the denominator is zero, then the vector (a,b,c) and the vector (i,j,k) are perpendicular. Note that (a,b,c) is the normal to the plane and (i,j,k) is the direction of the line. It follows that the line is either parallel to the plane or contained in the plane. In either case there is no unique intersection point.

Represent the ray as:

x = a + t b

where a, b and x are 3-vectors

Represent the quadric as:

x' A x + x' B + C = 0

where A is a symmetric 3x3 matrix, B a 3x1 vector, and C a scalar. Note that x' A x is the dot product of x and A x.

Substituting in gives the quadratic equation:

              2
        l2 * t  + l1 * t + l0 = 0

where:

	l2 = b' A b

        l1 = 2 b' A a + b' B

        l0 = a' A a + b' B + C

In this notation, If the determinant of this equation is less than 0, the line does not intersect the quadric. Solving the equation and substituting the values of t into the ray equation will give you the points.

The marching cubes algorithm is used in volume rendering to construct an isosurface from a 3D field of values.

The 2D analog would be to take an image, and for each pixel, set it to black if the value is below some threshold, and set it to white if it's above the threshold. Then smooth the jagged black outlines by skinning them with lines.

The marching cubes algorithm tests the corner of each cube (or voxel) in the scalar field as being either above or below a given threshold. This yields a collection of boxes with classified corners. Since there are eight corners with one of two states, there are 256 different possible combinations for each cube. Then, for each cube, you replace the cube with a surface that meets the classification of the cube. For example, the following are some 2D examples showing the cubes and their associated surface.

    - ----- +       - ----- -       - ----- +       - ----- +
    |:::'   |       |:::::::|       |::::   |       |   ':::|
    |:'     |       |:::::::|       |::::   |       |.    ':|
    |       |       |       |       |::::   |       |::.    |
    + ----- +       + ----- +       - ----- +       + ----- -

The result of the marching cubes algorithm is a smooth surface that approximates the isosurface that is constant along a given threshold. This is useful for displaying a volume of oil in a geological volume, for example.

Just define all points of all polygons in clockwise order.

            c = (x3*((z1*y2)-(y1*z2))+
                (y3*((x1*z2)-(z1*x2))+
                (z3*((y1*x2)-(x1*y2))

x1,y1,z1, x2,y2,z2, x3,y3,z3 = the first three points of the polygon.

If c is positive the polygon is visible. If c is negative the polygon is invisible (or the other way).

c = (x1-x2)*(y3-y2)-(y1-y2)*(x3-x2)

x1,y1, x2,y2, x3,y3 = the first three points of the polygon.

If c is positive the polygon is visible. If c is negative the polygon is invisible (or the other way).

This routine is `off the net somewhere'

unsigned short psi_sqrt(unsigned long v)
/*
// Calculates the square root of a 32-bit number.
*/
{
    register long t = 1L << 30, r = 0, s;

#define STEP(k) \
    s = t + r; \
    r >>= 1; \
    if (s <= v) { \
        v -= s; \
        r |= t; \
    }

    STEP(15); t >>= 2;
    STEP(14); t >>= 2;
    STEP(13); t >>= 2;
    STEP(12); t >>= 2;
    STEP(11); t >>= 2;
    STEP(10); t >>= 2;
    STEP(9); t >>= 2;
    STEP(8); t >>= 2;
    STEP(7); t >>= 2;
    STEP(6); t >>= 2;
    STEP(5); t >>= 2;
    STEP(4); t >>= 2;
    STEP(3); t >>= 2;
    STEP(2); t >>= 2;
    STEP(1); t >>= 2;
    STEP(0);

    return (unsigned short) r;
}

Sometimes you have a 3D vector and you want some vector perpendicular to it. Of course there are an infinite number of possible ones, but you just need any old perpendicular. The thing to do is to find some other vector that is not parallel to the given one and to calculate the cross product, which gives your answer. The tricky thing is guaranteeing that the "other" vector is not parallel to the original one and ideally that it is maximally non-parallel to ensure an accurate cross product. What we do is to find which of the principal axes is most perpendicular to the vector and then cross-product with that. Any other way will lead you to grief in the end. The following routine is not the fastest way of doing this but it is clear.

Vector3 get_perpendicular(const Vector3& A) 
// returns a line perpendicular to A
{
  // find smallest projection onto a principal axis
  // (ie the nearest to 90 degrees)
  // cross this with the original and return
  
  Vector3 X(1,0,0);
  Vector3 Y(0,1,0);
  Vector3 Z(0,0,1);

  Vector3& axis = Z;
  if (A[0] < A[1]) {
    // Not Y 
    if (A[0] < A[2])
      axis = X;
  }
  else
    // Not X
    if (A[1] < A[2]) {
      axis = Y;
    }
  
  return (A ^ axis).normalize();
}