[Visualpython-users] fix for twist in thick curve

[<Awa...@ao...>]

included in this email is a modifed routine to plot thick curves without 
twisting ...

(This code can replace the same routine in curve.cpp)

The only modification is to the code which calculates the orthogonal vectors 
which span the plane at the intersection of line segments.  The modification 
constructs orthogonal vectors which span the first plane (perpendicular to the 
first line segment) and then parallel transports these vectors along the 
curve.

regards

Andrew Wadsley

### start of code

void
curve::thickline( rView& view)
{
    // double* v, double* color, double radius, int count
    const int MAX_SEGMENTS = 1024;
    const int curve_around = 4;
    static vertex projected[MAX_SEGMENTS * curve_around ];
    static float light[MAX_SEGMENTS * curve_around * 3];
    int step = (count + MAX_SEGMENTS-1) / MAX_SEGMENTS;

    vertex *pr = projected;
    float *lt = light;

    float *cost = curve_sc;
    float *sint = cost + curve_around;

    //Vector lastx(0,1,0), lasty(0,0,1);
    vector x, y, lastx, lasty;

        vector lastB, B;
        float  tmp;


    int npoints = 0;
    
    // pos and color iterators
    const double* v_i = (double*)( data(pos));
    const double* c_i = (double*)( data(color));

    for(size_t corner=0; corner < count/step; ++corner, v_i += step*3, c_i += 
step*3) {
        // The vector to which v_i currently points towards.
        vector current( v_i[0], v_i[1], v_i[2] );
        view.ext_point( current);

        /* To construct the circle of vertices around a point on the curve,
         * I need an orthonormal basis of the plane of the circle.  A and B,
         * normalized vectors pointing toward the next and from the previous
         * points on the curve, are vectors in the plane of the *curve* at
         * this point.  The cross product and the difference of these vectors
         * are orthogonal and lie in the appropriate plane.
             * Then parallel transport the last frame to this plane
         */

        if (corner && (corner+1 < count/step)) {
            // Any point except the first and last.
            vector next( v_i[3], v_i[4], v_i[5]);    // The next vector in 
pos
            vector prev( v_i[-3], v_i[-2], v_i[-1]); // The previous vector 
in pos
            vector A = next - current;
            B = current - prev;
                  if ( !B ) B = lastB;

            x = A-B;
            y = A.cross(B);

            if (!x || !y) {
                x  = lastx;
                y  = lasty;
                      B  = lastB;
                  }
                  else {
                        vector pz = x.cross(y);
                        pz  = pz.norm();
                        tmp = pz.dot(lastx)/pz.dot(B);
                        x   = lastx - B * tmp;
                        tmp = pz.dot(lasty)/pz.dot(B);
                        y   = lasty - B * tmp;
                  }
        }
        else if (corner) {
            // The last point
            vector prev( v_i[-3], v_i[-2], v_i[-1]);
            B = current - prev;
                  if ( !B ) B = lastB;
            x = B.cross(vector(0,1,0));
            if (!x) x = B.cross(vector(1,0,0));
            y = B.cross(x);

            if (!x || !y) {
                x  = lastx;
                y  = lasty;
                      B  = lastB;
                  }
                  else {
                        vector pz = x.cross(y);
                        pz   = pz.norm();
                        tmp  = pz.dot(lastx)/pz.dot(B);
                        x    = lastx - B * tmp;
                        tmp  = pz.dot(lasty)/pz.dot(B);
                        y    = lasty - B * tmp;
                  }
        }
        else {
            // The first point
            vector next( v_i[3], v_i[4], v_i[5]);
            B = next - current;
            x = B.cross(vector(0,1,0));
            if (!x) x = B.cross(vector(1,0,0));
                        y  = B.cross(x);
                        if (!x || !y ) {
                            x = vector(1,0,0);
                            y = vector(0,0,1);
                            B = vector(0,1,0);
                        }
        }

            x      = x.norm();
            y      = y.norm();
        lastx  = x;
            lasty  = y;
            lastB  = B;

        float lx[2], ly[2];
        for(int li=0; li < view.lights.n_lights; li++) {
            lx[li] = view.lights.L[li].dot( x );
            ly[li] = view.lights.L[li].dot( y );
        }

        x = x * radius;
        y = y * radius;
        
        // Manipulate colors when we are in stereo mode.
        rgb rendered_color( static_cast<float>(c_i[0]), 
                            static_cast<float>(c_i[1]),
                            static_cast<float>(c_i[2]));
        if (view.anaglyph) {
            if (view.coloranaglyph) {
                rendered_color = rendered_color.unsaturate();
            }
            else {
                rendered_color.r = 
                rendered_color.g = 
                rendered_color.b = rendered_color.grayscale();
            }
        }
        
        npoints++;
        for (int a=0; a < curve_around; a++) {
            float c = cost[a];
            float s = sint[a];

            float illum = view.lights.ambient;
            for(int li=0; li < view.lights.n_lights; li++) {
                float dot = c*lx[li] + s*ly[li];
                if (dot > 0.0)
                    illum += dot;
            }

            view.wct.project(current + x*c + y*s, *pr++);
            *(lt++) = illum * rendered_color.r;
            *(lt++) = illum * rendered_color.g;
            *(lt++) = illum * rendered_color.b;
        }
    }

    glEnableClientState(GL_VERTEX_ARRAY);
    glEnableClientState(GL_COLOR_ARRAY);
    glShadeModel(GL_SMOOTH);

    int *ind = curve_slice;
    for(int a=0; a<curve_around; a++) {
        int ai = a;
        if (a == curve_around-1) {
            ind += 256;
            ai = 0;
        }

        for(int i=0; i<npoints; i+=127) {
            glVertexPointer(4, GL_DOUBLE, sizeof(vertex), &
projected[i*curve_around + ai].x);
            glColorPointer (3, GL_FLOAT, sizeof(float)*3, &
light[(i*curve_around + ai)*3] );

            if (npoints-i < 128)
                glDrawElements(GL_TRIANGLE_STRIP, 2*(npoints-i), 
GL_UNSIGNED_INT, ind);
            else
                glDrawElements(GL_TRIANGLE_STRIP, 256, GL_UNSIGNED_INT, ind);
        }
    }
}


### end of code