Month: March 2019

Barycentric Coordinates

thequietlearnerLeave a comment

An interesting concept I came across recently in an introductory geometry book is barycentric coordinates. The fundamental idea is defining the location of a point by its relationship to the points of a triangle. (Well, technically you can use any polygon, but we will keep things simple.)

So, each of the vertices of the triangle (A, B, and C) are given a sort of weight or pull factor. These three numbers uniquely represent the position of P. The theorem goes

If A, B, C are three non-collinear points, any vector P may be expressed in terms of the vectors A, B, C thus:

P = xA + yB +zC, where x + y + z = 1

Geometry: A Comprehensive Course, Dan Pedoe

If P is pulled all the way over to vertex A, then the coordinates become (1, 0, 0). Likewise at B it is (0, 1, 0) and for C it is (0, 0, 1). Another coordinate inside the triangle would be (⅓, ⅓, ⅓). If the point moves outside the triangle, some of those coordinates become negative.

Of course, instead of boring old vectors, you can use… complex numbers!

From the formula, it is obvious how to start with some barycentric coordinates, add a triangle, and find the point. Here is some Racket code:

(define (barycentric->complex x y z A B C)
   (+ (* x A) (* y B) (* z C)))

Note, that you can use any triangle you want with the same barycentric coordinates. This could allow you to stretch the graph of some set of coordinates in various ways.

But, how do you convert in the other direction, having a point and a triangle, and wanting the barycentric coordinates? This is a little more complicated. Fortunately, it happens that the factor for each of the vertices is equal to the proportional area of the triangle opposite that vertex. So, looking at this diagram again

x = [BCP] / [ABC], y = [CAP] / [ABC], z = [ABP] / [ABC]
where [BCP] is the area of triangle BCP, etc.

Now, a tricky part here is that the areas must be signed areas, or the formula doesn’t necessarily work correctly. Signed area means that you have to be consistent in assigning a negative or positive value to the area of the triangle, depending on the orientation of the vertices. So, if you happen to plug in the vertices in a counter-clockwise direction, you get a negative area, and in a clockwise direction, you get a positive area. Thankfully, there is an elegant formula for determining this, which I will represent by this Racket function

(define (orientation Z1 Z2 Z3)
  (let* ([rise1 (- (imag-part Z2) (imag-part Z1))]
         [rise2 (- (imag-part Z3) (imag-part Z2))]
         [run1 (- (real-part Z2) (real-part Z1))]
         [run2 (- (real-part Z3) (real-part Z2))]
         [diff (- (* rise1 run2) (* rise2 run1))])
    (if (< diff 0) 'ccw
        (if (> diff 0) 'cw 'col))))

The function returns ‘ccw for counter-clockwise, ‘cw for clockwise, and ‘col for collinear (i.e., the points are on a straight line, which is not relevant in this application).

This allows us to have a signed-area function, using the common formula for determining the area of a triangle from just the length of the sides. We get the length of the sides by subtracting the complex numbers from each other and pulling the magnitudes of the resulting complex numbers.

(define (signed-area A B C)
  (let* ([a (magnitude (- C B))]
         [b (magnitude (- A C))]
         [c (magnitude (- A B))]
         [s (* (+ a b c) 0.5)]
         [area (sqrt (* s (- s a) (- s b) (- s c)))]
         [orient (orientation A B C)])
    (if (eq? orient 'ccw)
        (* -1 area)
        area)))

And finally we determine x, y, and z:

(define (complex->barycentric A B C P)
  (let ([ABC (signed-area A B C)])
    (list (/ (signed-area B C P) ABC)
          (/ (signed-area C A P) ABC)
          (/ (signed-area A B P) ABC))))

See, simple as Pi! Heh, heh.

racket@barycentric.rkt> (complex->barycentric 0 6 2+2i 2+i)
 '(0.33333333333333376 0.16666666666666657 0.5000000000000009)
 racket@barycentric.rkt> (barycentric->complex 0.33333333333333376 0.16666666666666657 0.5000000000000009 0 6 2+2i)
 2.0000000000000013+1.0000000000000018i

So, you can see it works, though there is some very small inaccuracy introduced by the approximations that occur during the conversion operations.

I hope to explore, in future posts, some interesting applications of barycentric coordinates.

Advertisement

Trigonometry Problem without Trigonometry: Clarification/Correction

thequietlearner1 Comment

In the previous post I was exulting about the coolness of solving a trigonometry problem without trigonometry, using complex numbers. Well, that isn’t quite the full truth. You can certainly do the problem described easily and elegantly with a complex number capable calculator. However, there is some hidden trigonometry in the initial rotation operation:

The complex number ei has to be raised to a power. If you raise to an integer power, you could simply multiply out the rectangular values algebraically. But how do you raise a complex number to a fractional power, like 2.5? The standard approach, near as I can tell, is to convert the complex number to polar format, with a magnitude r and an angle ϕ. Assuming complex number a+bi, pulling out the magnitude requires only the sqrt function (pythagorean theorem) but pulling out the angle requires the arctangent of b/a. Then, use the formula Zn = rn (cos nϕ + i sin nϕ). I really don’t know if this is how all calculators do it, but here is the code behind the operations in the racket interpreter:

  Scheme_Complex *cb = (Scheme_Complex *)base;
  Scheme_Complex *ce = (Scheme_Complex *)exponent;
  double a, b, c, d, bm, ba, nm, na, r1, r2;
  int d_is_zero;

  if ((ce->i == zero) && !SCHEME_FLOATP(ce->r)) {
    if (SCHEME_INTP(ce->r) || SCHEME_BIGNUMP(ce->r))
      return scheme_generic_integer_power(base, ce->r);
  }

  a = scheme_get_val_as_double(cb->r);
  b = scheme_get_val_as_double(cb->i);
  c = scheme_get_val_as_double(ce->r);
  d = scheme_get_val_as_double(ce->i);
  d_is_zero = (ce->i == zero);

  bm = sqrt(a * a + b * b);
  ba = atan2(b, a);

  /* New mag & angle */
  nm = scheme_double_expt(bm, c) * exp(-(ba * d));
  if (d_is_zero) /* precision here can avoid NaNs */
    na = ba * c;
  else
    na = log(bm) * d + ba * c;

  r1 = nm * cos(na);
  r2 = nm * sin(na);

The code there is just slightly more complicated, as it has separate angle finding conditional branches for real number and complex number exponents.

Of course, you don’t need to deal with any of that if you are just manipulating complex numbers in the interpreter or calculator.

I feel this would be a great place to dive into the interesting subject of complex roots, but that I should hold off until I have a more thorough grasp of the subject.

Trigonometry Problem without Trigonometry

thequietlearnerLeave a comment

Here is a classic trigonometry problem:

With the satellite dish pointed up 13 degrees above the ground, how far up will the Line of Sight (LOS) be after passing over 100 meters of ground? (To save a few sentences in our post, we will ignore the fact that the satellite dish will usually start a few meters above the ground.)

This is simple using any calculator with the trigonometry function tangent:

tan(ϕ) = opp/adj, so height = 100 * tan 13° (after 13° is converted to radians).

If you didn’t want to use a trigonometric function, however, what could you do? Another approach is complex numbers, and it is pretty simple. You simply need to rotate ϕ degrees up the unit circle (a basic operation with complex numbers), get the real and imaginary parts of that complex number, and use that proportion to figure out the opposite side of your problem triangle.

So, after raising e^i by 13 degrees (converted to radians) multiply 100 x (imag Z / real Z). In programming code, it looks like this:

#lang racket

(define tau (* 2 pi))

(define rfact (/ tau 360))

(define (rise gnd deg)
  (let* ([rad (* deg rfact)]
         [Z (exp (* 0+1i rad))])
    (/ (* gnd (imag-part Z)) (real-part Z))))

> (rise 100 13)
23.08681911255631

It gives geometrically sane results also if you use negative degrees or negative distances:

> (rise 100 -13)
-23.08681911255631
> (rise -100 -13)
23.08681911255631
> (rise -100 13)
-23.08681911255631

What if you don’t have or want to use the exp function (to raise e to some value). You instead start with an approximation of e^i, which happens to be

0.5403023058681398+0.8414709848078965i

and then just raise that to ϕ.

Now, why would you want to use the complex number approach rather than the trig function? Well, to be honest I couldn’t think of a really compelling reason off the top of my head. An easy one would be if you were already using complex numbers for other reasons. Another would be if you wanted to use purely algebraic operations in your programming, say, on a computer board that did not have trigonometric functions built in. I suspect that the complex number approach is ultimately simpler and more efficient on the computational level, but since trigonometric functions are usually approximated with tables and with an algebraic exponential function, I couldn’t really say for sure off the top of my head.

In any case, it seems cool that complex numbers allow you to do an angle and distance trigonometry problem without trigonometry.

Harmongrams on Arduino Via Serial Commands

thequietlearner1 Comment

Previously I was experimenting with trying to write advanced programs in uLisp, embedded on an Arduino with a TFT display/touchscreen. Ultimately this proves not the best approach due to limitations in the memory and also the uLisp interpreter. A different approach is to keep the uLisp interpreter on the Arduino, but limit the on-board programming to a set of core interface commands, e.g., (drawline) and (fillscr). Then have the desktop computer or server sending uLisp commands over the serial connection. Of course, the uLisp commands can be constructed in Racket Scheme, allowing the power of the Racket syntax and libraries, and the great CPU and memory access, in controlling the arduino. I.e., remote control by the Mother Brain.

The demo I’ve programmed here is drawing the harmonograms slow enough to watch, so as to simulate the harmonograph experience. I do not have a camera setup suitable to show this in live video, but here are some snapshots:

The Racket code for the demo is available here:

ftp://lavender.qlfiles.net/Racket/arduino-harmon-demo.7z (md5sum 7e2bc9bc2ea25d048577b28c95c4fca5)

The project to adapt uLisp with TFT control primitives is here:

http://savannah.nongnu.org/projects/ulisp-tft

Though control over serial cables is doable for some applications, the longer-term goal would be to have them transmitted over Wi-Fi.

Complex Numbers in Emacs Calc Tutorial Video (Part I)

thequietlearner1 Comment

After learning the basics of how to use complex numbers, I find myself suddenly more interested in buttons on scientific calculators, and the built-in functions of calculator and programming software. Something which I thought might be of interest to my small (but distinguished!) audience would be a tutorial on using complex numbers in the Gnu/Emacs Calc program.

Please forgive the quiet audio — you’ll need to turn up the volume a lot. I’m dealing with a not-ideal microphone set-up at the moment, but I hope to have access to a better microphone set-up next week.

I have been advised by well-meaning relatives that, by continuing to upload boring tutorials and technical articles, it is likely that my readership will never expand beyond two dozen people. There is hope, they say, if I were instead to upload cute bunny rabbit photos. Consequently, I am going to start the practice of randomly attaching bunny photos taken from my mother’s Web site about rabbit farming.

I am hoping that, since we are related, my Mom will overlook this blatant act of copyright violation. This image is reproduced for scholarship purposes only under the “fair use” doctrine.