[Kstars-devel] XYZ Rotation matrix in 3.5.x

[James Bowlin]

I implemented the X-Y-Z rotation matrix for converting from celestial
coordinates to screen coordinates in the 3.5.x branch.

The good news is that it was very easy to implement.  The bad news is that
it did not significantly speed things up.  I used the QTime code that Jason
put in the 4.0 trunk in order to time the drawing of the Milky Way.  I chose
to time this because I was familiar with the code and I already had two
different versions of it in the program.  I didn't cull out any points so
the new version had to deal with 250 more points, roughly an 8% increase.
The new version also did line clipping but I arranged the display so this
only occurred in a few places.  

Here are typical results for drawing the outline of the Milky Way:

toXY() Milky Way took 7.23 ms  (old way)
 X-Y-Z Milky Way took 6.00 ms  (new way)

I put in a loop in each branch to run the code 100 times which gave me two
extra digits of timing resolution.

I did a simple timing experiment to compare the speed of a multiplication
to the speed of doing sin() and in this experiment the multiplication was
faster but not by much, maybe 20% or so.

So it appears that modern processors can do simple sin() and cos()
calculations almost as fast as multiplication.  This implies that there
is probably not much speed boost from caching sin() and cos() values,
at least for modern processors.

For completeness, I've included the new code below.  Only the first two
routines are new.  98% of slowToScreenXYZ() came from the original toXY().
Maybe there is a more efficient way of doing the matrix calculation, making
use of the SSE registers or something.  I tried removing the sqrt() from
both branches but this didn't significantly alter the timing results.

I'll leave the code in my 3.5.x branch so we can make further experiments
but for now, I don't think it is a worthwhile change.  It could possibly
be useful for people with processors that take a longer time doing trig
functions but I would like to verify that situation before proceeding.

I am certain that the newer code is being used because until I got the
unit vectors in calcRotMatrix() correct, things that were drawn with
the new code ended up in the wrong places on the screen.  It was neat to
see the grid moving along with the cursor and the constellations moving
in the opposite direction.


QPoint SkyMap::toScreen( SkyPoint *skyP, bool *clipped) {
    QPoint o;
    double Width = ( width() * screen_scale );
    double Height = ( height() * screen_scale );

    double cx = skyP->cx;
    double cy = skyP->cy;
    double cz = skyP->cz;
    
    double x = cx * rot_xx + cy * rot_yx + cz * rot_zx;
    double y = cx * rot_xy + cy * rot_yy + cz * rot_zy;
    double z = cx * rot_xz + cy * rot_yz + cz * rot_zz;
    double k = sqrt( 2.0/( 1 + z ) );

    o.setX( int( 0.5*Width  + k*x ) );
    o.setY( int( 0.5*Height + k*y ) );
    *clipped = ( z < 0.0 );
    return o;
}

void SkyMap::calcRotMatrix( bool Horiz, bool doRefraction, double scale) {
    SkyPoint px = SkyPoint( 0.0,  0.0 );
    SkyPoint py = SkyPoint( 6.0,  0.0 );
    SkyPoint pz = SkyPoint( 0.0, 90.0 );
    screen_scale = scale;
    slowToScreenXYZ( &px, &rot_xx, &rot_xy, &rot_xz);
    slowToScreenXYZ( &py, &rot_yx, &rot_yy, &rot_yz);
    slowToScreenXYZ( &pz, &rot_zx, &rot_zy, &rot_zz);
}

void SkyMap::slowToScreenXYZ( SkyPoint *o, double *sx, double *sy, double *sz) {

    bool Horiz = false;  //Options::useAltAz();
    bool doRefraction = false;

    double Y, dX;
    double sindX, cosdX, sinY, cosY, sinY0, cosY0;

    double pscale = Options::zoomFactor() * screen_scale;

    if ( Horiz ) {
        if ( doRefraction ) Y = refract( o->alt(), true ).radians(); //account for atmospheric refraction
        else Y = o->alt()->radians();

        if ( focus()->az()->Degrees() > 270.0 && o->az()->Degrees() < 90.0 ) {
            dX = 2*dms::PI + focus()->az()->radians() - o->az()->radians();
        } else {
            dX = focus()->az()->radians() - o->az()->radians();
        }

        focus()->alt()->SinCos( sinY0, cosY0 );

    } else {
        if (focus()->ra()->Hours() > 18.0 && o->ra()->Hours() < 6.0) {
	    dX = 2*dms::PI + o->ra()->radians() - focus()->ra()->radians();
        } else {
            dX = o->ra()->radians() - focus()->ra()->radians();
        }
        Y = o->dec()->radians();
        focus()->dec()->SinCos( sinY0, cosY0 );
    }

    //Convert dX, Y coords to screen pixel coords.
    #if ( __GLIBC__ >= 2 && __GLIBC_MINOR__ >=1 ) && !defined(__UCLIBC__)
    //GNU version
    sincos( dX, &sindX, &cosdX );
    sincos( Y, &sinY, &cosY );
    #else
    //ANSI version
    sindX = sin(dX);
    cosdX = cos(dX);
    sinY  = sin(Y);
    cosY  = cos(Y);
    #endif

    *sz = sinY0*sinY + cosY0*cosY*cosdX;
    
    //double k = sqrt( 2.0/( 1 + z ) );
    double k = 1.0;
    
    *sx = - pscale*k*cosY*sindX ;
    *sy = - pscale*k*( cosY0*sinY - sinY0*cosY*cosdX ) ;
}

-- 
Peace, James