Projecting and un-projecting 3D point onto sphere

[Rick C. Hodgin]

I have an Open GL project I'm working on and it's almost working. I
can't figure out what I'm doing wrong on this last part.

I have a quad with 4 points oriented in this manner:

1 _ 2
|_|
3 4

I have a sphere that's oriented like a globe with the north pole up. I
use two angles I call rho and theta. Rho goes from 0 at the north pole
to 180 degrees at the south pole. Theta goes from 0 and the meridian
around to 360 degrees back to itself.

I have the quad oriented above the north pole and project it around the
sphere. The link below shows an image which has an array of these quads
in an early and incorrect projection that has since been corrected, but
the image illustrates generally what I'm attempting.

Loading Image...

To project the quad I use these trig sequences. The iRotateX() and
iRotateZ() algorithms work correctly, and rotate the image around those
axes. The axes are oriented left/right = x, up/down = y, to/fro = z.

struct SVert { double x, y, z; };

void iProjectVert(SVert* v, double radius)
{
double x, y, z, t1, t2, t3;
t1 = v->x / radius;
t2 = v->z / radius;
t3 = t1 / cos(t2); // Adjust great circle to lesser circle

// Position starting point at the north pole
x = 0.0;
y = radius;
z = 0.0;

// Rotate on X axis to move it to/fro based on z
iRotateX(x, y, z, t2);

// Rotate on the lesser circle based for x and y
iRotateZ(x, y, z, t3);

// Store the result
v->x = x;
v->y = y;
v->z = z;
}

The above seems to work. I am able to project the test pattern used in
the attached image onto the sphere and I believe it's correct, but it
may not be. I'm not sure how to validate it other than what I think it
should look like. :-) I suppose a calculation of surface area on each
projected quad onto the sphere. What's the formula for computing the
area of a general 4-point polygon on a sphere?

When I go to un-project those points back to their flat plane above the
north pole, I think this is where my math is wrong.

void iUnprojectVert(SVert* v,
double rho, double theta, double radius)
{
double t1, t2, x, y, z;

// Compute thetas
t1 = atan2(v->z, v->y);
t2 = atan2(v->x, v->y);

// Un-project
v->x = radius * t2 * cos(t1);
v->y = radius;
v->z = radius * t1 * cos(t2);
}

It's very close. Quads that are very small (like by proportion the size
of Australia compared to the size of the Earth) are nearly perfect. But
the further it wraps around the sphere the more distorted it gets.

I have no doubt my math is wrong somewhere, but I don't know how to fix
it. I need some help here please. Any help is appreciated.

Thank you in advance.

[Hans-Bernhard Bröker]

The above seems to work.  I am able to project the test pattern used in
the attached image onto the sphere and I believe it's correct, but it
may not be.

To be able to tell, you would have to have a separate way of expressing
what it is you're trying to achieve. "It should look look like ..."
doesn't exactly qualify for that.

I suppose a calculation of surface area on each
projected quad onto the sphere.

That might qualify as a criterion if the unstated goal you set out to
achieve was an area-preserving projection. But be warned that goal is
impossible for non-infinitesimal quadrangles. That's one way of saying
that wall-papering a sphere is simply not possible if you want to avoid
all wrinkles and tears.

What's the formula for computing the
area of a general 4-point polygon on a sphere?

I suspect there isn't one. There is one for triangles, but IIRC that
only works if they don't include either pole. Quadrangles have to be
split into triangles, but the decision how to do that correctly in the
general case is ... complicated.

    void iUnprojectVert(SVert* v,
                           double rho, double theta, double radius)
    {
        double t1, t2, x, y, z;
        // Compute thetas
        t1 = atan2(v->z, v->y);
        t2 = atan2(v->x, v->y);
        // Un-project
        v->x = radius * t2 * cos(t1);
        v->y = radius;
        v->z = radius * t1 * cos(t2);
    }

If you want to revert your original projection, why don't you do just
that, step by step?

t3 = atan2(v->x, v->z); // recover t3 from rotation around Y
rho = sqrt(v->x*v->x + v->z*v->z); // projected radius
t2 = atan2(rho,v->y) // recover t3 from rotation around z
t1 = t3 * cos(t2);
v->x = t1 * radius;
v->z = t2 * radius;

[Rick C. Hodgin]

Post by Hans-Bernhard Bröker

     void iUnprojectVert(SVert* v,
                            double rho, double theta, double radius)
     {
         double t1, t2, x, y, z;
         // Compute thetas
         t1 = atan2(v->z, v->y);
         t2 = atan2(v->x, v->y);
         // Un-project
         v->x = radius * t2 * cos(t1);
         v->y = radius;
         v->z = radius * t1 * cos(t2);
     }

If you want to revert your original projection, why don't you do just
that, step by step?
    t3 = atan2(v->x, v->z); // recover t3 from rotation around Y
    rho = sqrt(v->x*v->x + v->z*v->z);  // projected radius
    t2 = atan2(rho,v->y)  // recover t3 from rotation around z
    t1 = t3 * cos(t2);
    v->x = t1 * radius;
    v->z = t2 * radius;

I tried this, but it didn't work. I wasn't sure why. After much
debugging, I discovered that I had the projection wrong for my display,
which meant the rotation I had for the X,Y,Z axes was also backwards.

I was able to step through it using a grid for the spherical texture
showing degrees (radians) until I got it sorted through all the steps.

Now that I've done that, my original algorithms are all working. They
just needed theta to be theta and not 2PI-theta. :-)

I always appreciate your help, Hans-Bernhard Bröker. Thank you for your
reply. And thank you for trying to make my life a little better by
incorporating some of your own knowledge, skill, and expertise into my
project.

If there's ever anything I can do for you, please ask.

[Scott Hemphill]

Post by Hans-Bernhard Bröker

The above seems to work.  I am able to project the test pattern used
in the attached image onto the sphere and I believe it's correct,
but it may not be.

To be able to tell, you would have to have a separate way of
expressing what it is you're trying to achieve. "It should look look
like ..." doesn't exactly qualify for that.

I suppose a calculation of surface area on each projected quad onto
the sphere.

That might qualify as a criterion if the unstated goal you set out to
achieve was an area-preserving projection. But be warned that goal is
impossible for non-infinitesimal quadrangles. That's one way of
saying that wall-papering a sphere is simply not possible if you want
to avoid all wrinkles and tears.

What's the formula for computing the area of a general 4-point
polygon on a sphere?

I suspect there isn't one. There is one for triangles, but IIRC that
only works if they don't include either pole. Quadrangles have to be
split into triangles, but the decision how to do that correctly in the
general case is ... complicated.

[snip]

There is a simple formula for the area of a polygon on a sphere, if you
have the angles at each of the corners. It doesn't depend on a
coordinate system, and therefore there are no issues with poles.

Area of triangle = r^2 * (sum of angles - pi)

For a general polygon with n sides, area = r^2 * (sum of angles -
pi*(n-2)).

Note that the polygon doesn't have to be convex.

Scott