#help_index "Graphics/Math"

public I64  gr_x_offsets[8]  = {-1,  0,  1, -1, 1, -1, 0, 1}, 
            gr_y_offsets[8]  = {-1, -1, -1,  0, 0,  1, 1, 1}, 
            gr_x_offsets2[4] = { 0, -1,  1,  0}, 
            gr_y_offsets2[4] = {-1,  0,  0,  1};

public Bool Line(U8 *aux_data, I64 x1, I64 y1, I64 z1, I64 x2, I64 y2, I64 z2, 
                Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z),
                I64 step=1, I64 start=0)
{//Step through line segment calling callback.
//Uses fixed-point.
    I64  i, j, d, dx = x2 - x1, dy = y2 - y1, dz = z2 - z1, _x, _y, _z, adx = AbsI64(dx), ady = AbsI64(dy), adz = AbsI64(dz);
    Bool first = TRUE;

    if (adx >= ady)
    {
        if (adx >= adz)
        {
            if (d = adx)
            {
                if (dx >= 0)
                    dx = 0x100000000;
                else
                    dx = -0x100000000;
                dy = dy << 32 / d;
                dz = dz << 32 / d;
            }
        }
        else
        {
            if (d = adz)
            {
                dx = dx << 32 / d;
                dy = dy << 32 / d;
                if (dz >= 0)
                    dz = 0x100000000;
                else
                    dz = -0x100000000;
            }
        }
    }
    else
    {
        if (ady >= adz)
        {
            if (d = ady)
            {
                dx = dx << 32 / d;
                if (dy >= 0)
                    dy = 0x100000000;
                else
                    dy = -0x100000000;
                dz = dz << 32 / d;
            }
        }
        else
        {
            if (d = adz)
            {
                dx = dx << 32 / d;
                dy = dy << 32 / d;
                if (dz >= 0)
                    dz = 0x100000000;
                else
                    dz = -0x100000000;
            }
        }
    }
    x1 <<= 32;
    y1 <<= 32;
    z1 <<= 32;
    for (j =0; j < start; j++)
    {
        x1 += dx;
        y1 += dy;
        z1 += dz;
    }
    if (step != 1 && step != 0)
    {
        dx *= step;
        dy *= step;
        dz *= step;
        d  /= step;
    }
    for (i = start; i <= d; i++)
    {
        if ((_x != x1.i32[1] || _y != y1.i32[1] || _z != z1.i32[1] || first) &&
                !(*fp_plot)(aux_data, x1.i32[1], y1.i32[1], z1.i32[1]))
            return FALSE;
        first = FALSE;
        _x = x1.i32[1];
        _y = y1.i32[1];
        _z = z1.i32[1];
        x1 += dx; y1+=dy; z1+=dz;
    }
    if (step == 1 && (_x != x2 || _y != y2 || _z != z2) && !(*fp_plot)(aux_data, x2, y2, z2))
        return FALSE;

    return TRUE;
}

#help_index "Graphics/Math/3D Transformation"
public I64 *Mat4x4MulMat4x4Equ(I64 *dst, I64 *m1, I64 *m2)
{//Multiply 4x4 matrices and store in dst. Uses fixed-point.
//Conceptually,  the transform m1 is applied after m2
    I64 i, j, k;
    F64 sum;

    for (i = 0; i < 4; i++)
    {
        for (j = 0; j < 4; j++)
        {
            sum = 0;
            for (k = 0; k < 4; k++)
                sum += ToF64(m1[k + 4 * j]) * ToF64(m2[i + 4 * k]);
            dst[i + 4 * j] = sum / GR_SCALE;
        }
    }

    return dst;
}

public I64 *Mat4x4MulMat4x4New(I64 *m1, I64 *m2, CTask *mem_task=NULL)
{//Multiply 4x4 matrices. Return MAlloced matrix. Uses fixed-point.
//Conceptually, the transform m1 is applied after m2
    return Mat4x4MulMat4x4Equ(MAlloc(sizeof(I64) * 16, mem_task), m1, m2);
}

public I64 *Mat4x4Equ(I64 *dst, I64 *src)
{//Copy 4x4 Rot matrix.
    MemCopy(dst, src, sizeof(I64) * 16);

    return dst;
}

public I64 *Mat4x4New(I64 *src, CTask *mem_task=NULL)
{//Return MAlloced copy of 4x4 Rot matrix.
    return Mat4x4Equ(MAlloc(sizeof(I64) * 16, mem_task), src);
}

public I64 *Mat4x4RotX(I64 *m, F64 phi)
{//Rot matrix about X axis. Uses fixed-point.
    F64 my_cos = Cos(phi) * GR_SCALE, my_sin = Sin(phi) * GR_SCALE;
    I64 r[16], r2[16];

    MemSet(r, 0, sizeof(r));
    r[5]  = my_cos;
    r[10] = my_cos;
    r[9]  = my_sin;
    r[6]  = -my_sin;
    r[0]  = GR_SCALE;
    r[15] = GR_SCALE;

    return Mat4x4Equ(m, Mat4x4MulMat4x4Equ(r2, r, m));
}

public I64 *Mat4x4RotY(I64 *m, F64 omega)
{//Rot matrix about Y axis. Uses fixed-point.
    F64 my_cos = Cos(omega) * GR_SCALE, my_sin = Sin(omega) * GR_SCALE;
    I64 r[16], r2[16];

    MemSet(r, 0, sizeof(r));
    r[0]  = my_cos;
    r[10] = my_cos;
    r[8]  = -my_sin;
    r[2]  = my_sin;
    r[5]  = GR_SCALE;
    r[15] = GR_SCALE;

    return Mat4x4Equ(m, Mat4x4MulMat4x4Equ(r2, r, m));
}

public I64 *Mat4x4RotZ(I64 *m, F64 theta)
{//Rot matrix about Z axis. Uses fixed-point.
    F64 my_cos=Cos(theta)*GR_SCALE, my_sin=Sin(theta)*GR_SCALE;
    I64 r[16], r2[16];

    MemSet(r, 0, sizeof(r));
    r[0]  = my_cos;
    r[5]  = my_cos;
    r[4]  = my_sin;
    r[1]  = -my_sin;
    r[10] = GR_SCALE;
    r[15] = GR_SCALE;

    return Mat4x4Equ(m, Mat4x4MulMat4x4Equ(r2, r, m));
}

public I64 *Mat4x4Scale(I64 *m, F64 s)
{//Scale 4x4 matrix by value.
    I64 i;

    for (i = 0; i < 16; i++)
        m[i] *= s;

    return m;
}

public U0 DCThickScale(CDC *dc=gr.dc)
{//Scale device context's thick by norm of transformation.
    I64 d;

    if (dc->flags & DCF_TRANSFORMATION)
    {
        if (dc->thick)
        {
            d = dc->thick * dc->r_norm + 0x80000000; //Round
            dc->thick = d.i32[1];
            if (dc->thick < 1)
                dc->thick = 1;
        }
    }
}

public I64 *Mat4x4TranslationEqu(I64 *r, I64 x, I64 y, I64 z)
{//Set translation values in 4x4 matrix. Uses fixed-point.
    r[0 * 4 + 3] = x << 32;
    r[1 * 4 + 3] = y << 32;
    r[2 * 4 + 3] = z << 32;
    r[3 * 4 + 3] = GR_SCALE;

    return r;
}

public I64 *Mat4x4TranslationAdd(I64 *r, I64 x, I64 y, I64 z)
{//Add translation to 4x4 matrix. Uses fixed-point.
    r[0 * 4 + 3] += x << 32;
    r[1 * 4 + 3] += y << 32;
    r[2 * 4 + 3] += z << 32;
    r[3 * 4 + 3] = GR_SCALE;

    return r;
}

public Bool DCSymmetrySet(CDC *dc=gr.dc, I64 x1, I64 y1, I64 x2, I64 y2)
{//2D. Set device context's symmetry.
    F64 d;

    if (y1 == y2 && x1 == x2)
        return FALSE;
    dc->sym.snx = y2 - y1;
    dc->sym.sny = x1 - x2;
    dc->sym.snz = 0;
    if (d = Sqrt(SqrI64(dc->sym.snx) + SqrI64(dc->sym.sny) + SqrI64(dc->sym.snz)))
    {
        d = GR_SCALE / d;
        dc->sym.snx *= d;
        dc->sym.sny *= d;
        dc->sym.snz *= d;
    }
    dc->sym.sx = x1;
    dc->sym.sy = y1;
    dc->sym.sz = 0;

    return TRUE;
}

public Bool DCSymmetry3Set(CDC *dc=gr.dc, I64 x1, I64 y1, I64 z1, I64 x2, I64 y2, I64 z2, I64 x3, I64 y3, I64 z3)
{//3D. Set device context's symmetry.
    F64 d, x, y, z, xx, yy, zz;
    I64 xx1, yy1, zz1, xx2, yy2, zz2, *r;

    xx1 = x1 - x2;
    yy1 = y1 - y2;
    zz1 = z1 - z2;
    xx2 = x3 - x2;
    yy2 = y3 - y2;
    zz2 = z3 - z2;

    if (!xx1 && !yy1 && !zz1 || !xx2 && !yy2 && !zz2 || xx1 == xx2 && yy1 == yy2 && zz1 == zz2)
        return FALSE;

    x = yy1  * zz2 - zz1 * yy2;
    y = -xx1 * zz2 + zz1 * xx2;
    z = xx1  * yy2 - yy1 * xx2;
    if (dc->flags & DCF_TRANSFORMATION)
    {
        r = dc->r;
        xx = x * r[0] + y * r[1] + z * r[2];
        yy = x * r[4] + y * r[5] + z * r[6];
        zz = x * r[8] + y * r[9] + z * r[10];
        x = xx;
        y = yy;
        z = zz;
    }
    if (d = Sqrt(Sqr(x) + Sqr(y) + Sqr(z)))
    {
        d = GR_SCALE / d;
        dc->sym.snx = d * x;
        dc->sym.sny = d * y;
        dc->sym.snz = d * z;
    }
    if (dc->flags & DCF_TRANSFORMATION)
        (*dc->transform)(dc, &x1, &y1, &z1);
    dc->sym.sx = x1;
    dc->sym.sy = y1;
    dc->sym.sz = z1;

    return TRUE;
}

public U0 DCReflect(CDC *dc, I64 *_x, I64 *_y, I64 *_z)
{//Reflect 3D point about device context's symmetry. Uses fixed-point.
    I64 x = *_x << 32,
        y = *_y << 32,
        z = *_z << 32, 
        xx = *_x - dc->sym.sx,
        yy = *_y - dc->sym.sy,
        zz = *_z - dc->sym.sz, 
        d = (xx * dc->sym.snx + yy * dc->sym.sny + zz * dc->sym.snz) >> 16, 
        xn, yn, zn, xx2, yy2, zz2;

    xn = d * dc->sym.snx >> 15;
    yn = d * dc->sym.sny >> 15;
    zn = d * dc->sym.snz >> 15;
    xx = x - xn;
    yy = y - yn;
    zz = z - zn;
    xx2 = x + xn;
    yy2 = y + yn;
    zz2 = z + zn;
    if (SqrI64(xx >> 16 - dc->sym.sx << 16) + SqrI64(yy >> 16 - dc->sym.sy << 16) + SqrI64(zz >> 16 - dc->sym.sz << 16) <
        SqrI64(xx2 >> 16 - dc->sym.sx << 16) + SqrI64(yy2 >> 16 - dc->sym.sy << 16) + SqrI64(zz2 >> 16 - dc->sym.sz << 16))
    {
        *_x = xx.i32[1];
        *_y = yy.i32[1];
        *_z = zz.i32[1];
    }
    else
    {
        *_x = xx2.i32[1];
        *_y = yy2.i32[1];
        *_z = zz2.i32[1];
    }
}

#help_index "Graphics/Math"
#define GR_SCALE1_BITS  24
#define GR_SCALE2_BITS  8
public Bool Circle(U8 *aux_data, I64 cx, I64 cy, I64 cz, I64 radius, Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z), 
                   I64 step=1, F64 start_radians=0, F64 len_radians=2*pi)
{//Step through circle arc calling callback.
    I64 i, j, len = Ceil(len_radians * radius), x, y, x1, y1, s1, s2, c;
    F64 t;

    if (radius <= 0 || !step)
        return TRUE;
    t = 1.0 / radius;
    c = 1 << GR_SCALE1_BITS * Cos(t);
    if (step < 0)
    {
        step = -step;
        s2 = 1 << GR_SCALE1_BITS * Sin(t);
        s1 = -s2;
    }
    else
    {
        s1 = 1 << GR_SCALE1_BITS * Sin(t);
        s2 = -s1;
    }
    if (start_radians)
    {
        x = radius * Cos(start_radians);
        y = -radius * Sin(start_radians);
    }
    else
    {
        x = radius;
        y = 0;
    }
    x <<= GR_SCALE2_BITS;
    y <<= GR_SCALE2_BITS;
    for (i = 0; i <= len; i += step)
    {
        if (!(*fp_plot)(aux_data, cx + x >> GR_SCALE2_BITS, cy + y >> GR_SCALE2_BITS, cz))
            return FALSE;
        for (j = 0; j < step; j++)
        {
            x1 =(c * x + s1 * y) >> GR_SCALE1_BITS;
            y1 =(s2 * x + c * y) >> GR_SCALE1_BITS;
            x = x1;
            y = y1;
        }
    }

    return TRUE;
}

public Bool Ellipse(U8 *aux_data, I64 cx, I64 cy, I64 cz, I64 x_radius, I64 y_radius,
                    Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z), F64 rot_angle=0,
                    I64 step=1, F64 start_radians=0, F64 len_radians=2*pi)
{//Step through ellipse arc calling callback.
    I64  i, j, len, x, y, _x, _y, x1, y1, x2, y2,  s1, s2, c,  s12, s22, c2;
    F64  t;
    Bool first = TRUE;

    if (x_radius <= 0 || y_radius <= 0 || !step)
        return TRUE;
    if (x_radius >= y_radius)
    {
        t = 1.0 / x_radius;
        len = Ceil(len_radians * x_radius);
    }
    else
    {
        t = 1.0 / y_radius;
        len = Ceil(len_radians * y_radius);
    }

    c = 1 << GR_SCALE1_BITS * Cos(t);
    if (step < 0)
    {
        step = -step;
        s2 =1 << GR_SCALE1_BITS * Sin(t);
        s1 =-s2;
    }
    else
    {
        s1 = 1 << GR_SCALE1_BITS * Sin(t);
        s2 = -s1;
    }

    c2 = 1 << GR_SCALE1_BITS * Cos(rot_angle);
    s12 = 1 << GR_SCALE1_BITS * Sin(rot_angle);
    s22 = -s12;

    if (start_radians)
    {
        x = x_radius * Cos(start_radians);
        y = -x_radius * Sin(start_radians);
    }
    else
    {
        x = x_radius;
        y = 0;
    }
    x <<= GR_SCALE2_BITS;
    y <<= GR_SCALE2_BITS;
    x2 = x;
    y2 = y;

    y1 = y2 * y_radius / x_radius;
    x = (c2 * x2 + s12 * y1) >> GR_SCALE1_BITS;
    y = (s22 * x2 + c2 * y1) >> GR_SCALE1_BITS;

    for (i = 0; i <= len; i += step)
    {
        if ((x >> GR_SCALE2_BITS != _x || y >> GR_SCALE2_BITS != _y || first) &&
                !(*fp_plot)(aux_data, cx + x >> GR_SCALE2_BITS, cy + y >> GR_SCALE2_BITS, cz))
            return FALSE;

        _x = x >> GR_SCALE2_BITS;
        _y = y >> GR_SCALE2_BITS;
        first = FALSE;
        for (j = 0; j < step; j++)
        {
            x1 = (c * x2 + s1 * y2) >> GR_SCALE1_BITS;
            y1 = (s2 * x2 + c * y2) >> GR_SCALE1_BITS;
            x2 = x1;
            y2 = y1;
            y1 = y1 * y_radius/x_radius;
            x = (c2 * x1+ s12 * y1) >> GR_SCALE1_BITS;
            y = (s22 * x1 + c2 * y1) >> GR_SCALE1_BITS;
        }
    }

    return TRUE;
}

public Bool RegPoly(U8 *aux_data, I64 cx, I64 cy, I64 cz, I64 x_radius, I64 y_radius, I64 sides, 
                    Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z), 
                    F64 rot_angle=0, I64 step=1, F64 start_radians=0, F64 len_radians=2*pi)
{//Step through regular polygon calling callback.
    I64 i, n, x, y, x1, y1, x2, y2, xx1, yy1, xx2, yy2, s1, s2, c,  s12, s22, c2;
    F64 angle_step;

    if (sides <= 0 || x_radius <= 0 || y_radius <= 0)
        return TRUE;

    angle_step = 2 * pi / sides;
    n = len_radians * sides / (2 * pi);

    s1 = 1 << GR_SCALE1_BITS * Sin(angle_step);
    s2 = -s1;
    c = 1 << GR_SCALE1_BITS * Cos(angle_step);

    s12 = 1 << GR_SCALE1_BITS * Sin(rot_angle);
    s22 = -s12;
    c2 = 1 << GR_SCALE1_BITS * Cos(rot_angle);

    if (start_radians)
    {
        x = x_radius * Cos(start_radians);
        y = -x_radius * Sin(start_radians);
    }
    else
    {
        x = x_radius;
        y = 0;
    }
    x <<= GR_SCALE2_BITS;
    y <<= GR_SCALE2_BITS;
    x2 = x;
    y2 = y;

    y1 = y2 * y_radius / x_radius;
    x = (c2 * x2 + s12 * y1) >> GR_SCALE1_BITS;
    y = (s22 * x2 + c2 * y1) >> GR_SCALE1_BITS;

    xx1 = cx + x >> GR_SCALE2_BITS;
    yy1 = cy + y >> GR_SCALE2_BITS;
    for (i = 0; i < n; i++)
    {
        x1 = (c * x2 + s1 * y2) >> GR_SCALE1_BITS;
        y1 = (s2 * x2 + c * y2) >> GR_SCALE1_BITS;
        x2 = x1;
        y2 = y1;
        y1 = y1 * y_radius / x_radius;
        x = (c2 * x1 + s12 * y1) >> GR_SCALE1_BITS;
        y = (s22 * x1 + c2 * y1) >> GR_SCALE1_BITS;
        xx2 = cx + x >> GR_SCALE2_BITS;
        yy2 = cy + y >> GR_SCALE2_BITS;
        if (!Line(aux_data, xx1, yy1, cz, xx2, yy2, cz, fp_plot, step))
            return FALSE;
        xx1 = xx2;
        yy1 = yy2;
    }

    return TRUE;
}

#help_index "Graphics/Data Types/D3I32;Math/Data Types/D3I32;Data Types/D3I32"
public F64 D3I32Dist(CD3I32 *p1, CD3I32 *p2)
{//Distance
    return Sqrt(SqrI64(p1->x - p2->x) + SqrI64(p1->y - p2->y) + SqrI64(p1->z - p2->z));
}

public I64 D3I32DistSqr(CD3I32 *p1, CD3I32 *p2)
{//Distance Squared
    return SqrI64(p1->x - p2->x) + SqrI64(p1->y - p2->y) + SqrI64(p1->z - p2->z);
}

public F64 D3I32Norm(CD3I32 *p)
{//Norm
    return Sqrt(SqrI64(p->x) + SqrI64(p->y) + SqrI64(p->z));
}

public I64 D3I32NormSqr(CD3I32 *p)
{//Norm Squared
    return SqrI64(p->x) + SqrI64(p->y) + SqrI64(p->z);
}

#help_index "Graphics/Math"
public Bool Bezier2(U8 *aux_data, CD3I32 *ctrl, Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z), Bool first=TRUE)
{//Go in 2nd order bezier calling callback.
    I64 x, y, z, xx, yy, zz, dx, dy, dz, d_max;
    F64 x0 = ctrl[0].x,
        y0 = ctrl[0].y,
        z0 = ctrl[0].z, 
        x1 = ctrl[1].x - x0,
        y1 = ctrl[1].y - y0,
        z1 = ctrl[1].z - z0, 
        x2 = ctrl[2].x - x0,
        y2 = ctrl[2].y - y0,
        z2 = ctrl[2].z - z0, 
        t, d = D3I32Dist(&ctrl[0], &ctrl[1]) + D3I32Dist(&ctrl[1], &ctrl[2]) + D3I32Dist(&ctrl[2], &ctrl[0]), 
        s = 0.5 / d, t1, t2;

    xx = x0;
    yy = y0;
    zz = z0;
    if (first && !(*fp_plot)(aux_data, xx, yy, zz))
        return FALSE;
    for (t = 0.0; t <= 1.0; t += s)
    {
        t1 = t * (1.0 - t);
        t2 = t * t;
        x = x0 + x1 * t1 + x2 * t2;
        y = y0 + y1 * t1 + y2 * t2;
        z = z0 + z1 * t1 + z2 * t2;
        dx = AbsI64(x - xx);
        dy = AbsI64(y - yy);
        dz = AbsI64(z - zz);
        if (dx > dy)
            d_max = dx;
        else
            d_max = dy;
        if (dz > d_max)
            d_max = dz;
        if (!d_max)
            s *= 1.1;
        else
        {
            s *= 0.9;
            if (!(*fp_plot)(aux_data, x, y, z))
                return FALSE;
            xx = x;
            yy = y;
            zz = z;
        }
    }
    x = ctrl[2].x;
    y = ctrl[2].y;
    z = ctrl[2].z;
    if ((xx != x || yy != y || zz != z) && !(*fp_plot)(aux_data, x, y, z))
        return FALSE;

    return TRUE;
}

public Bool Bezier3(U8 *aux_data, CD3I32 *ctrl, Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z), Bool first=TRUE)
{//Go in 3rd order bezier calling callback.
    I64 x, y, z, xx, yy, zz, dx, dy, dz, d_max;
    F64 x0 = ctrl[0].x,
        y0 = ctrl[0].y,
        z0 = ctrl[0].z, 
        x1 = ctrl[1].x - x0,
        y1 = ctrl[1].y - y0,
        z1 = ctrl[1].z - z0, 
        x2 = ctrl[2].x - x0,
        y2 = ctrl[2].y - y0,
        z2 = ctrl[2].z - z0, 
        x3 = ctrl[3].x - x0,
        y3 = ctrl[3].y - y0,
        z3 = ctrl[3].z - z0, 
        t, d =  D3I32Dist(&ctrl[0], &ctrl[1]) +
                D3I32Dist(&ctrl[1], &ctrl[2]) +
                D3I32Dist(&ctrl[2], &ctrl[3]) +
                D3I32Dist(&ctrl[3], &ctrl[0]),
        s = 0.5 / d, nt, t1, t2, t3;

    xx = x0;
    yy = y0;
    zz = z0;
    if (first && !(*fp_plot)(aux_data, xx, yy, zz))
        return FALSE;
    for (t = 0.0; t <= 1.0; t += s)
    {
        nt = 1.0 - t;
        t1 = t * nt * nt;
        t2 = t * t * nt;
        t3 = t * t * t;
        x = x0 + x1 * t1 + x2 * t2 + x3 * t3;
        y = y0 + y1 * t1 + y2 * t2 + y3 * t3;
        z = z0 + z1 * t1 + z2 * t2 + z3 * t3;
        dx = AbsI64(x - xx);
        dy = AbsI64(y - yy);
        dz = AbsI64(z - zz);
        if (dx > dy)
            d_max = dx;
        else
            d_max = dy;
        if (dz > d_max)
            d_max = dz;
        if (!d_max)
            s *= 1.1;
        else
        {
            s *= 0.9;
            if (!(*fp_plot)(aux_data, x, y, z))
                return FALSE;
            xx = x;
            yy = y;
            zz = z;
        }
    }
    x = ctrl[3].x;
    y = ctrl[3].y;
    z = ctrl[3].z;
    if ((xx != x || yy != y || zz != z) &&!(*fp_plot)(aux_data, x, y, z))
        return FALSE;

    return TRUE;
}

public Bool BSpline2(U8 *aux_data, CD3I32 *ctrl, I64 count, Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z), Bool closed=FALSE)
{//Go in 2nd order spline calling callback.
    I64      i, j;
    CD3I32  *c;
    Bool     first;

    if (count < 3)
        return FALSE;
    first = TRUE;
    if (closed)
    {
        count++;
        c = MAlloc(sizeof(CD3I32) * (count * 2 - 1));
        j = 1;
        for (i = 0; i < count - 2; i++)
        {
            c[j].x = (ctrl[i].x + ctrl[i + 1].x) / 2.0;
            c[j].y = (ctrl[i].y + ctrl[i + 1].y) / 2.0;
            c[j].z = (ctrl[i].z + ctrl[i + 1].z) / 2.0;
            j += 2;
        }
        c[j].x = (ctrl[0].x + ctrl[count - 2].x) / 2.0;
        c[j].y = (ctrl[0].y + ctrl[count - 2].y) / 2.0;
        c[j].z = (ctrl[0].z + ctrl[count - 2].z) / 2.0;

        c[0].x = (c[1].x + c[j].x) / 2.0;
        c[0].y = (c[1].y + c[j].y) / 2.0;
        c[0].z = (c[1].z + c[j].z) / 2.0;
        j = 2;
        for (i = 0; i < count - 2; i++)
        {
            c[j].x = (c[j - 1].x + c[j + 1].x) / 2.0;
            c[j].y = (c[j - 1].y + c[j + 1].y) / 2.0;
            c[j].z = (c[j - 1].z + c[j + 1].z) / 2.0;
            j += 2;
        }
        c[j].x = c[0].x;
        c[j].y = c[0].y;
        c[j].z = c[0].z;
    }
    else
    {
        c = MAlloc(sizeof(CD3I32) * (count * 2 - 1));
        c[0].x = ctrl[0].x;
        c[0].y = ctrl[0].y;
        c[0].z = ctrl[0].z;
        c[count * 2 - 2].x = ctrl[count - 1].x;
        c[count * 2 - 2].y = ctrl[count - 1].y;
        c[count * 2 - 2].z = ctrl[count - 1].z;
        j = 1;
        for (i = 0; i < count - 1; i++)
        {
            c[j].x = (ctrl[i].x + ctrl[i + 1].x) / 2.0;
            c[j].y = (ctrl[i].y + ctrl[i + 1].y) / 2.0;
            c[j].z = (ctrl[i].z + ctrl[i + 1].z) / 2.0;
            j += 2;
        }
        j = 2;
        for (i = 0; i < count - 2; i++)
        {
            c[j].x = (c[j - 1].x + c[j + 1].x) / 2.0;
            c[j].y = (c[j - 1].y + c[j + 1].y) / 2.0;
            c[j].z = (c[j - 1].z + c[j + 1].z) / 2.0;
            j += 2;
        }
    }
    for (i = 0; i < count * 2 - 2; i += 2)
    {
        if (!Bezier2(aux_data, &c[i], fp_plot, first))
            return FALSE;
        first = FALSE;
    }
    Free(c);

    return TRUE;
}

public Bool BSpline3(U8 *aux_data, CD3I32 *ctrl, I64 count, Bool (*fp_plot)(U8 *aux, I64 x, I64 y, I64 z), Bool closed=FALSE)
{//Go in 3rd order spline calling callback.
    I64      i, j;
    F64      x, y, z;
    CD3I32  *c;
    Bool     first;

    if (count < 3)
        return FALSE;
    first = TRUE;
    if (closed)
    {
        count++;
        c = MAlloc(sizeof(CD3I32) * (count * 3 - 2));
        j = 1;
        for (i = 0; i < count - 2; i++)
        {
            x = ctrl[i].x;
            y = ctrl[i].y;
            z = ctrl[i].z;
            c[j].x = (ctrl[i + 1].x - x) / 3.0 + x;
            c[j].y = (ctrl[i + 1].y - y) / 3.0 + y;
            c[j].z = (ctrl[i + 1].z - z) / 3.0 + z;
            j++;
            c[j].x = 2.0 * (ctrl[i + 1].x - x) / 3.0 + x;
            c[j].y = 2.0 * (ctrl[i + 1].y - y) / 3.0 + y;
            c[j].z = 2.0 * (ctrl[i + 1].z - z) / 3.0 + z;
            j += 2;
        }
        x = ctrl[count - 2].x;
        y = ctrl[count - 2].y;
        z = ctrl[count - 2].z;
        c[j].x = (ctrl[0].x - x) / 3.0 + x;
        c[j].y = (ctrl[0].y - y) / 3.0 + y;
        c[j].z = (ctrl[0].z - z) / 3.0 + z;
        j++;
        c[j].x = 2.0 * (ctrl[0].x - x) / 3.0 + x;
        c[j].y = 2.0 * (ctrl[0].y - y) / 3.0 + y;
        c[j].z = 2.0 * (ctrl[0].z - z) / 3.0 + z;

        c[0].x = (c[1].x + c[j].x) / 2.0;
        c[0].y = (c[1].y + c[j].y) / 2.0;
        c[0].z = (c[1].z + c[j].z) / 2.0;

        j = 3;
        for (i = 0; i < count - 2; i++)
        {
            c[j].x = (c[j - 1].x + c[j + 1].x) / 2.0;
            c[j].y = (c[j - 1].y + c[j + 1].y) / 2.0;
            c[j].z = (c[j - 1].z + c[j + 1].z) / 2.0;
            j += 3;
        }
        c[j].x = c[0].x;
        c[j].y = c[0].y;
        c[j].z = c[0].z;
    }
    else
    {
        c = MAlloc(sizeof(CD3I32) * (count * 3 - 2));
        c[0].x = ctrl[0].x;
        c[0].y = ctrl[0].y;
        c[0].z = ctrl[0].z;
        c[count * 3 - 3].x = ctrl[count - 1].x;
        c[count * 3 - 3].y = ctrl[count - 1].y;
        c[count * 3 - 3].z = ctrl[count - 1].z;
        j = 1;
        for (i = 0; i < count - 1; i++)
        {
            x = ctrl[i].x;
            y = ctrl[i].y;
            z = ctrl[i].z;
            c[j].x = (ctrl[i + 1].x - x) / 3.0 + x;
            c[j].y = (ctrl[i + 1].y - y) / 3.0 + y;
            c[j].z = (ctrl[i + 1].z - z) / 3.0 + z;
            j++;
            c[j].x = 2.0 * (ctrl[i + 1].x - x) / 3.0 + x;
            c[j].y = 2.0 * (ctrl[i + 1].y - y) / 3.0 + y;
            c[j].z = 2.0 * (ctrl[i + 1].z - z) / 3.0 + z;
            j += 2;
        }
        j = 3;
        for (i = 0; i < count - 2; i++)
        {
            c[j].x = (c[j - 1].x + c[j + 1].x) / 2.0;
            c[j].y = (c[j - 1].y + c[j + 1].y) / 2.0;
            c[j].z = (c[j - 1].z + c[j + 1].z) / 2.0;
            j += 3;
        }
    }
    for (i = 0; i < count * 3 - 3; i += 3)
    {
        if (!Bezier3(aux_data, &c[i], fp_plot, first))
            return FALSE;
        first = FALSE;
    }
    Free(c);

    return TRUE;
}

#define CC_LEFT         1
#define CC_RIGHT        2
#define CC_TOP          4
#define CC_BOTTOM       8

public Bool ClipLine(I64 *_x1, I64 *_y1, I64 *_x2, I64 *_y2, I64 left, I64 top, I64 right, I64 bottom)
{//Clip x1,y1 x2,y2 with left,top,right,bottom.
    I64 x, y, x1 = *_x1, y1 = *_y1, x2 = *_x2, y2 = *_y2, cc, cc1, cc2;

    if (y1 > bottom)
        cc1 = CC_BOTTOM;
    else if (y1 < top)
        cc1 = CC_TOP;
    else
        cc1 = 0;
    if (x1 > right)
        cc1 |= CC_RIGHT;
    else if (x1 < left)
        cc1 |= CC_LEFT;

    if (y2 > bottom)
        cc2 = CC_BOTTOM;
    else if (y2 < top)
        cc2 = CC_TOP;
    else
        cc2 = 0;
    if (x2 > right)
        cc2 |= CC_RIGHT;
    else if (x2 < left)
        cc2 |= CC_LEFT;

    while (TRUE)
    {
        if (!(cc1 | cc2))
            return TRUE;
        if (cc1 & cc2)
            return FALSE;

        if (cc1)
            cc = cc1;
        else
            cc = cc2;

        if (cc & CC_BOTTOM)
        {
            x = x1 + (x2 - x1) * (bottom - y1) / (y2 - y1);
            y = bottom;
        }
        else if (cc & CC_TOP)
        {
            x = x1 + (x2 - x1) * (top - y1) / (y2 - y1);
            y = top;
        }
        else if (cc & CC_RIGHT)
        {
            y = y1 + (y2 - y1) * (right - x1) / (x2 - x1);
            x = right;
        }
        else
        {
            y = y1 + (y2 - y1) * (left - x1) / (x2 - x1);
            x = left;
        }

        if (cc == cc1)
        {
            *_x1 = x1 = x;
            *_y1 = y1 = y;
            if (y1 > bottom)
                cc1 = CC_BOTTOM;
            else if (y1 < top)
                cc1 = CC_TOP;
            else
                cc1 = 0;
            if (x1 > right)
                cc1 |= CC_RIGHT;
            else if (x1 < left)
                cc1 |= CC_LEFT;
        }
        else
        {
            *_x2 = x2 = x;
            *_y2 = y2 = y;
            if (y2 > bottom)
                cc2 = CC_BOTTOM;
            else if (y2 < top)
                cc2 = CC_TOP;
            else
                cc2 = 0;
            if (x2 > right)
                cc2 |= CC_RIGHT;
            else if (x2 < left)
                cc2 |= CC_LEFT;
        }
    }
}