• Projecting and un-projecting 3D point onto sphere

    From Rick C. Hodgin@21:1/5 to All on Wed Nov 11 20:20:22 2020
    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.

    http://www.libsf.org/images/sphere_wrap.png

    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.

    --
    Rick C. Hodgin

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From =?UTF-8?Q?Hans-Bernhard_Br=c3=b6ker@21:1/5 to All on Thu Nov 12 21:47:17 2020
    Am 12.11.2020 um 02:20 schrieb Rick C. Hodgin:

    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);
    x = t1 * radius;
    z = t2 * radius;

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Rick C. Hodgin@21:1/5 to All on Sun Nov 15 19:57:53 2020
    On 11/12/20 3:47 PM, Hans-Bernhard Bröker wrote:
    Am 12.11.2020 um 02:20 schrieb Rick C. Hodgin:
         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.

    --
    Rick C. Hodgin

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Scott Hemphill@21:1/5 to HBBroeker@t-online.de on Wed Nov 25 13:38:59 2020
    Hans-Bernhard Bröker <HBBroeker@t-online.de> writes:

    Am 12.11.2020 um 02:20 schrieb Rick C. Hodgin:

    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
    --
    Scott Hemphill hemphill@alumni.caltech.edu
    "This isn't flying. This is falling, with style." -- Buzz Lightyear

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)