Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I found this picture on mathoverflow, which I find very intriguing and so I like to know how to draw such an image with a simple computer program.

To calculate the rational point, I can draw a line from P_0(0,0,1) and P_1(u,v,0) and calculate the intersection with the sphere as follows:

\begin{equation} x=\frac{2u}{u^2+v^2+1};y=\frac{2v}{u^2+v^2+1};z=\frac{u^2+v^2-1}{u^2+v^2+1} \end{equation}

So the intersection coordinates are rational numbers where

\begin{equation} a=2u;b=2v;c=u^2+v^2-1;d=u^2+v^2+1 \end{equation}

with

\begin{equation} a^2+b^2+c^2=d^2 \end{equation}

If I understand correctly the guy who posted the picture he would mark the points depending on the value of d. So a value of d below a certain threshold would override the pixel color to white.

Now inside a program we would loop through all pixels along the x and y axis using the following algorithm:

for(int y = 0; y < height; ++y)
{
    for(int x = 0; x < width; ++x)
    {
        // does pixel ray intersect with sphere?
        // using orthgraphic projection ( all rays are parallel ) of 
        // sphere with radius of image height
        if(intersect(sphere,x,y))
        {
            // calculate the distance of the intersection of 
            // pixel ray and sphere
            double z = calc_distance(sphere,x,y);

            // the further away the darker the color
            rgb_color color = make_color(z);

            // calculate d
            // the gcd of a,b,c,d will be used to make d as small as possible
            int d = calculate_denominator(x,y)

            // 100 is made up
            if(d < 100)
            {
                color = white;
            }

            set_pixel(img,x,y,color);
        }
    }
}

Could someone help in correcting my algorithm?

share|improve this question

1 Answer 1

up vote 2 down vote accepted

I'd iterate over $u$ and $v$ instead, since the point which actually has easy rational coordinates is not exactly at a pixel position.

double d;
for (int u = 0; u <= 100; ++u) {
  for (int v = 0; (d = u*u+v*v+1) <= 100; ++v) {
    markBlack(2*u/d, 2*v/d); // arguments are double in [0,1]
  }
}

or, if you want to stay in integer arithmetic:

int d;
for (int u = 0; u <= 100; ++u) {
  for (int v = 0; (d = u*u+v*v+1) <= 100; ++v) {
    int dh = d/2; // rounded down
    int x = (size*2*u + dh)/d;
    int y = (size*2*v + dh)/d;
    set_pixel(img,x,y,black);
  }
}

In both cases, I'm ommitting the third coordinate since we are projecting onto a coordinate plane.

The picture still doesn't look like the one you indicated, though. That's likely because you'd have to consider fractional $u,v$ as well. You might be better off using Pythagorean quadruples instead. Then you can use nested loops there.

int size = 2047, limit = 100;
int dm, dn, dp, d;
for (int m = 0; (dm = m*m) <= limit; ++m)
  for (int n = 0; (dn = dm + n*n) <= limit; ++n)
    for (int p = 0; (dp = dn + p*p) <= limit; ++p)
      for (int q = 0; (d = (dp + q*q) <= limit; ++q) {
        if (d == 0) continue;
        int a = abs(m*m + n*n - p*p - q*q);
        int b =     2*(m*q + n*p) ;
        int c = abs(2*(n*q - m*p));
        int dh = d/2;
        int x = (size*a + dh)/d;
        int y = (size*b + dh)/d;
        int z = (size*c + dh)/d;
        set_pixel(img,x,y,black);
        set_pixel(img,x,z,black);
        set_pixel(img,y,x,black);
        set_pixel(img,y,z,black);
        set_pixel(img,z,x,black);
        set_pixel(img,z,y,black);
      }

However, instead of marking pixels either black or white, you might also consider marking them with a suitable gray level, with the gray value depending on the number of rational points within a given pixel. If you do that, then you want to make sure that you mark every rational point exactly once. I did so in the code used for the animation below, building on this paper.

Animation

share|improve this answer
    
Thanks for your solution. Once I understand it I'll mark your answer. Could you post the whole code so I can test it out? I tried the the second loop but don't know what size means. –  chhenning Aug 31 at 16:24
    
@chhenning: size is the size of the image, in pixels, Actually one less, since coordinates range from 0 to size inclusive. I never ran this C code, but I ran its Python analog in sage, and posted that here. –  MvG Aug 31 at 19:23

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.