/*
 * The Npic library
 *
 * Copyright (C) 2003 Edouard Thiel <Edouard.Thiel@lif.univ-mrs.fr>
 *
 * This library is free software under the terms of the 
 * GNU Lesser General Public License (LGPL) version 2.1.
*/

/* DO NOT EDIT !!!  Generated by npic-templa from "calc_bresenham.ct" */

/*
 * calc_bresenham.c - 05/01/2009
 *
 * Draw lines with generalized Bresenham's algorithm
 *
 * References:
 *
 * [1] J. E. Bresenham, 1965. "Algorithm for computer control of a digital 
 *     plotter", IBM Systems Journal, Vol. 4, No.1, pp. 25-30.
*/


#include <npic.h>


/*--------------------- P U B L I C - I N T E R F A C E ----------------------*/

/*
 * Draw a discrete line bounded by P and Q, with pixel values varying from
 * cP to cQ. P and Q are inside image, else nothing is done.
 *
 * Do not modify other pixels neither border pixels.
 * Do nothing if np is not ok. set not ok on error. Verbose.
*/


void NpicDrawLine_2c (Npic_image *np, 
        int yP, int xP, 
        int yQ, int xQ, 
        Npic_c cP, Npic_c cQ)
{
    Npic_image_2c *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        jm, j;

    /* We need signed integers for sc, etc */
    int c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_2C) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; 
    c = cP; p1->pix[y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_2l (Npic_image *np, 
        int yP, int xP, 
        int yQ, int xQ, 
        Npic_l cP, Npic_l cQ)
{
    Npic_image_2l *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        jm, j;

    Npic_l c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_2L) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; 
    c = cP; p1->pix[y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "cP = %"NPIC_PL"  cQ = %"NPIC_PL"  c = %"NPIC_PL"  dc = %"NPIC_PL"  sc = %"NPIC_PL"  mc = %"NPIC_PL"  k1c = %"NPIC_PL"  k2c = %"NPIC_PL"  kc = %"NPIC_PL"\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_2d (Npic_image *np, 
        int yP, int xP, 
        int yQ, int xQ, 
        Npic_d cP, Npic_d cQ)
{
    Npic_image_2d *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        jm, j;

    Npic_d c, dc = fabs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_2D) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; 
    c = cP; p1->pix[y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || fabs(c-cQ) > 1E-8) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "cP = %lf  cQ = %lf  c = %lf  dc = %lf  sc = %lf  mc = %lf  k1c = %lf  k2c = %lf  kc = %lf\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_3c (Npic_image *np, 
        int zP, int yP, int xP, 
        int zQ, int yQ, int xQ, 
        Npic_c cP, Npic_c cQ)
{
    Npic_image_3c *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        jm, j;

    /* We need signed integers for sc, etc */
    int c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_3C) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; 
    c = cP; p1->pix[z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_3l (Npic_image *np, 
        int zP, int yP, int xP, 
        int zQ, int yQ, int xQ, 
        Npic_l cP, Npic_l cQ)
{
    Npic_image_3l *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        jm, j;

    Npic_l c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_3L) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; 
    c = cP; p1->pix[z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "cP = %"NPIC_PL"  cQ = %"NPIC_PL"  c = %"NPIC_PL"  dc = %"NPIC_PL"  sc = %"NPIC_PL"  mc = %"NPIC_PL"  k1c = %"NPIC_PL"  k2c = %"NPIC_PL"  kc = %"NPIC_PL"\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_3d (Npic_image *np, 
        int zP, int yP, int xP, 
        int zQ, int yQ, int xQ, 
        Npic_d cP, Npic_d cQ)
{
    Npic_image_3d *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        jm, j;

    Npic_d c, dc = fabs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_3D) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; 
    c = cP; p1->pix[z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || fabs(c-cQ) > 1E-8) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "cP = %lf  cQ = %lf  c = %lf  dc = %lf  sc = %lf  mc = %lf  k1c = %lf  k2c = %lf  kc = %lf\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_4c (Npic_image *np, 
        int tP, int zP, int yP, int xP, 
        int tQ, int zQ, int yQ, int xQ, 
        Npic_c cP, Npic_c cQ)
{
    Npic_image_4c *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        jm, j;

    /* We need signed integers for sc, etc */
    int c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_4C) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; 
    c = cP; p1->pix[t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_4l (Npic_image *np, 
        int tP, int zP, int yP, int xP, 
        int tQ, int zQ, int yQ, int xQ, 
        Npic_l cP, Npic_l cQ)
{
    Npic_image_4l *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        jm, j;

    Npic_l c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_4L) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; 
    c = cP; p1->pix[t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "cP = %"NPIC_PL"  cQ = %"NPIC_PL"  c = %"NPIC_PL"  dc = %"NPIC_PL"  sc = %"NPIC_PL"  mc = %"NPIC_PL"  k1c = %"NPIC_PL"  k2c = %"NPIC_PL"  kc = %"NPIC_PL"\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_4d (Npic_image *np, 
        int tP, int zP, int yP, int xP, 
        int tQ, int zQ, int yQ, int xQ, 
        Npic_d cP, Npic_d cQ)
{
    Npic_image_4d *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        jm, j;

    Npic_d c, dc = fabs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_4D) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; 
    c = cP; p1->pix[t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || fabs(c-cQ) > 1E-8) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "cP = %lf  cQ = %lf  c = %lf  dc = %lf  sc = %lf  mc = %lf  k1c = %lf  k2c = %lf  kc = %lf\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_5c (Npic_image *np, 
        int sP, int tP, int zP, int yP, int xP, 
        int sQ, int tQ, int zQ, int yQ, int xQ, 
        Npic_c cP, Npic_c cQ)
{
    Npic_image_5c *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        jm, j;

    /* We need signed integers for sc, etc */
    int c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_5C) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; 
    c = cP; p1->pix[s][t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[s][t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_5l (Npic_image *np, 
        int sP, int tP, int zP, int yP, int xP, 
        int sQ, int tQ, int zQ, int yQ, int xQ, 
        Npic_l cP, Npic_l cQ)
{
    Npic_image_5l *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        jm, j;

    Npic_l c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_5L) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; 
    c = cP; p1->pix[s][t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[s][t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);
        fprintf (stderr, "cP = %"NPIC_PL"  cQ = %"NPIC_PL"  c = %"NPIC_PL"  dc = %"NPIC_PL"  sc = %"NPIC_PL"  mc = %"NPIC_PL"  k1c = %"NPIC_PL"  k2c = %"NPIC_PL"  kc = %"NPIC_PL"\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_5d (Npic_image *np, 
        int sP, int tP, int zP, int yP, int xP, 
        int sQ, int tQ, int zQ, int yQ, int xQ, 
        Npic_d cP, Npic_d cQ)
{
    Npic_image_5d *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        jm, j;

    Npic_d c, dc = fabs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_5D) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; 
    c = cP; p1->pix[s][t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[s][t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || fabs(c-cQ) > 1E-8) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);
        fprintf (stderr, "cP = %lf  cQ = %lf  c = %lf  dc = %lf  sc = %lf  mc = %lf  k1c = %lf  k2c = %lf  kc = %lf\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_6c (Npic_image *np, 
        int rP, int sP, int tP, int zP, int yP, int xP, 
        int rQ, int sQ, int tQ, int zQ, int yQ, int xQ, 
        Npic_c cP, Npic_c cQ)
{
    Npic_image_6c *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        r, dr = abs (rQ-rP), sr = (rQ >= rP) ? 1:-1, mr, kr, k1r, k2r, h1r, h2r,
        jm, j;

    /* We need signed integers for sc, etc */
    int c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_6C) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
      || rP < 0 || rP >= p1->rmax || rQ < 0 || rQ >= p1->rmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; r = rP; 
    c = cP; p1->pix[r][s][t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;
    if (jm < dr) jm = dr;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    mr = dr / jm; k1r = 2*dr - 2*mr*jm; kr = k1r - jm; k2r = kr - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1r = mr*sr ; h2r = (mr+1)*sr;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (kr < 0) { kr += k1r ; r += h1r; } else { kr += k2r; r += h2r; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[r][s][t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || r != rQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);
        fprintf (stderr, "rP = %d  rQ = %d  r = %d  dr = %d  sr = %d  mr = %d  k1r = %d  k2r = %d\n  kr = %d",
                rP, rQ, r, dr, sr, mr, k1r, k2r, kr);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_6l (Npic_image *np, 
        int rP, int sP, int tP, int zP, int yP, int xP, 
        int rQ, int sQ, int tQ, int zQ, int yQ, int xQ, 
        Npic_l cP, Npic_l cQ)
{
    Npic_image_6l *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        r, dr = abs (rQ-rP), sr = (rQ >= rP) ? 1:-1, mr, kr, k1r, k2r, h1r, h2r,
        jm, j;

    Npic_l c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_6L) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
      || rP < 0 || rP >= p1->rmax || rQ < 0 || rQ >= p1->rmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; r = rP; 
    c = cP; p1->pix[r][s][t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;
    if (jm < dr) jm = dr;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    mr = dr / jm; k1r = 2*dr - 2*mr*jm; kr = k1r - jm; k2r = kr - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1r = mr*sr ; h2r = (mr+1)*sr;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (kr < 0) { kr += k1r ; r += h1r; } else { kr += k2r; r += h2r; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[r][s][t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || r != rQ || c != cQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);
        fprintf (stderr, "rP = %d  rQ = %d  r = %d  dr = %d  sr = %d  mr = %d  k1r = %d  k2r = %d\n  kr = %d",
                rP, rQ, r, dr, sr, mr, k1r, k2r, kr);
        fprintf (stderr, "cP = %"NPIC_PL"  cQ = %"NPIC_PL"  c = %"NPIC_PL"  dc = %"NPIC_PL"  sc = %"NPIC_PL"  mc = %"NPIC_PL"  k1c = %"NPIC_PL"  k2c = %"NPIC_PL"  kc = %"NPIC_PL"\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_6d (Npic_image *np, 
        int rP, int sP, int tP, int zP, int yP, int xP, 
        int rQ, int sQ, int tQ, int zQ, int yQ, int xQ, 
        Npic_d cP, Npic_d cQ)
{
    Npic_image_6d *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        r, dr = abs (rQ-rP), sr = (rQ >= rP) ? 1:-1, mr, kr, k1r, k2r, h1r, h2r,
        jm, j;

    Npic_d c, dc = fabs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_6D) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
      || rP < 0 || rP >= p1->rmax || rQ < 0 || rQ >= p1->rmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; r = rP; 
    c = cP; p1->pix[r][s][t][z][y][x] = c;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;
    if (jm < dr) jm = dr;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    mr = dr / jm; k1r = 2*dr - 2*mr*jm; kr = k1r - jm; k2r = kr - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1r = mr*sr ; h2r = (mr+1)*sr;
    h1c = mc*sc ; h2c = (mc+1)*sc;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (kr < 0) { kr += k1r ; r += h1r; } else { kr += k2r; r += h2r; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }

        p1->pix[r][s][t][z][y][x] = c;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || r != rQ || fabs(c-cQ) > 1E-8) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);
        fprintf (stderr, "rP = %d  rQ = %d  r = %d  dr = %d  sr = %d  mr = %d  k1r = %d  k2r = %d\n  kr = %d",
                rP, rQ, r, dr, sr, mr, k1r, k2r, kr);
        fprintf (stderr, "cP = %lf  cQ = %lf  c = %lf  dc = %lf  sc = %lf  mc = %lf  k1c = %lf  k2c = %lf  kc = %lf\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);

   }
}


void NpicDrawLine_2q (Npic_image *np, 
        int yP, int xP, 
        int yQ, int xQ, 
        int aP, int bP, int cP, int dP, int aQ, int bQ, int cQ, int dQ)
{
    Npic_image_2q *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        jm, j;

    /* We need signed integers */
    int a, da = abs (aQ-aP), sa = (aQ >= aP) ? 1:-1, ma, ka, k1a, k2a, h1a, h2a,
        b, db = abs (bQ-bP), sb = (bQ >= bP) ? 1:-1, mb, kb, k1b, k2b, h1b, h2b,
        c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c,
        d, dd = abs (dQ-dP), sd = (dQ >= dP) ? 1:-1, md, kd, k1d, k2d, h1d, h2d;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_2Q) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; 
    a = aP; b = bP; c = cP; d = dP; 
    p1->pix[y][x].a = a; p1->pix[y][x].b = b;
    p1->pix[y][x].c = c; p1->pix[y][x].d = d;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    ma = da / jm; k1a = 2*da - 2*ma*jm; ka = k1a - jm; k2a = ka - jm;
    mb = db / jm; k1b = 2*db - 2*mb*jm; kb = k1b - jm; k2b = kb - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;
    md = dd / jm; k1d = 2*dd - 2*md*jm; kd = k1d - jm; k2d = kd - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1a = ma*sa ; h2a = (ma+1)*sa;
    h1b = mb*sb ; h2b = (mb+1)*sb;
    h1c = mc*sc ; h2c = (mc+1)*sc;
    h1d = md*sd ; h2d = (md+1)*sd;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (ka < 0) { ka += k1a ; a += h1a; } else { ka += k2a; a += h2a; }
        if (kb < 0) { kb += k1b ; b += h1b; } else { kb += k2b; b += h2b; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }
        if (kd < 0) { kd += k1d ; d += h1d; } else { kd += k2d; d += h2d; }

        p1->pix[y][x].a = a; p1->pix[y][x].b = b;
        p1->pix[y][x].c = c; p1->pix[y][x].d = d;
    }

    /* Debug code */
    if (x != xQ || y != yQ || a != aQ || b != bQ || c != cQ || d != dQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);

         fprintf (stderr, "aP = %d  aQ = %d  a = %d  da = %d  sa = %d  ma = %d  k1a = %d  k2a = %d  ka = %d\n",
                aP, aQ, a, da, sa, ma, k1a, k2a, ka);
        fprintf (stderr, "bP = %d  bQ = %d  b = %d  db = %d  sb = %d  mb = %d  k1b = %d  k2b = %d  kb = %d\n",
                bP, bQ, b, db, sb, mb, k1b, k2b, kb);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);
        fprintf (stderr, "dP = %d  dQ = %d  d = %d  dd = %d  sd = %d  md = %d  k1d = %d  k2d = %d  kd = %d\n",
                dP, dQ, d, dd, sd, md, k1d, k2d, kd);
   }
}


void NpicDrawLine_3q (Npic_image *np, 
        int zP, int yP, int xP, 
        int zQ, int yQ, int xQ, 
        int aP, int bP, int cP, int dP, int aQ, int bQ, int cQ, int dQ)
{
    Npic_image_3q *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        jm, j;

    /* We need signed integers */
    int a, da = abs (aQ-aP), sa = (aQ >= aP) ? 1:-1, ma, ka, k1a, k2a, h1a, h2a,
        b, db = abs (bQ-bP), sb = (bQ >= bP) ? 1:-1, mb, kb, k1b, k2b, h1b, h2b,
        c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c,
        d, dd = abs (dQ-dP), sd = (dQ >= dP) ? 1:-1, md, kd, k1d, k2d, h1d, h2d;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_3Q) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; 
    a = aP; b = bP; c = cP; d = dP; 
    p1->pix[z][y][x].a = a; p1->pix[z][y][x].b = b;
    p1->pix[z][y][x].c = c; p1->pix[z][y][x].d = d;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    ma = da / jm; k1a = 2*da - 2*ma*jm; ka = k1a - jm; k2a = ka - jm;
    mb = db / jm; k1b = 2*db - 2*mb*jm; kb = k1b - jm; k2b = kb - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;
    md = dd / jm; k1d = 2*dd - 2*md*jm; kd = k1d - jm; k2d = kd - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1a = ma*sa ; h2a = (ma+1)*sa;
    h1b = mb*sb ; h2b = (mb+1)*sb;
    h1c = mc*sc ; h2c = (mc+1)*sc;
    h1d = md*sd ; h2d = (md+1)*sd;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (ka < 0) { ka += k1a ; a += h1a; } else { ka += k2a; a += h2a; }
        if (kb < 0) { kb += k1b ; b += h1b; } else { kb += k2b; b += h2b; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }
        if (kd < 0) { kd += k1d ; d += h1d; } else { kd += k2d; d += h2d; }

        p1->pix[z][y][x].a = a; p1->pix[z][y][x].b = b;
        p1->pix[z][y][x].c = c; p1->pix[z][y][x].d = d;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || a != aQ || b != bQ || c != cQ || d != dQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);

         fprintf (stderr, "aP = %d  aQ = %d  a = %d  da = %d  sa = %d  ma = %d  k1a = %d  k2a = %d  ka = %d\n",
                aP, aQ, a, da, sa, ma, k1a, k2a, ka);
        fprintf (stderr, "bP = %d  bQ = %d  b = %d  db = %d  sb = %d  mb = %d  k1b = %d  k2b = %d  kb = %d\n",
                bP, bQ, b, db, sb, mb, k1b, k2b, kb);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);
        fprintf (stderr, "dP = %d  dQ = %d  d = %d  dd = %d  sd = %d  md = %d  k1d = %d  k2d = %d  kd = %d\n",
                dP, dQ, d, dd, sd, md, k1d, k2d, kd);
   }
}


void NpicDrawLine_4q (Npic_image *np, 
        int tP, int zP, int yP, int xP, 
        int tQ, int zQ, int yQ, int xQ, 
        int aP, int bP, int cP, int dP, int aQ, int bQ, int cQ, int dQ)
{
    Npic_image_4q *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        jm, j;

    /* We need signed integers */
    int a, da = abs (aQ-aP), sa = (aQ >= aP) ? 1:-1, ma, ka, k1a, k2a, h1a, h2a,
        b, db = abs (bQ-bP), sb = (bQ >= bP) ? 1:-1, mb, kb, k1b, k2b, h1b, h2b,
        c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c,
        d, dd = abs (dQ-dP), sd = (dQ >= dP) ? 1:-1, md, kd, k1d, k2d, h1d, h2d;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_4Q) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; 
    a = aP; b = bP; c = cP; d = dP; 
    p1->pix[t][z][y][x].a = a; p1->pix[t][z][y][x].b = b;
    p1->pix[t][z][y][x].c = c; p1->pix[t][z][y][x].d = d;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ma = da / jm; k1a = 2*da - 2*ma*jm; ka = k1a - jm; k2a = ka - jm;
    mb = db / jm; k1b = 2*db - 2*mb*jm; kb = k1b - jm; k2b = kb - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;
    md = dd / jm; k1d = 2*dd - 2*md*jm; kd = k1d - jm; k2d = kd - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1a = ma*sa ; h2a = (ma+1)*sa;
    h1b = mb*sb ; h2b = (mb+1)*sb;
    h1c = mc*sc ; h2c = (mc+1)*sc;
    h1d = md*sd ; h2d = (md+1)*sd;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ka < 0) { ka += k1a ; a += h1a; } else { ka += k2a; a += h2a; }
        if (kb < 0) { kb += k1b ; b += h1b; } else { kb += k2b; b += h2b; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }
        if (kd < 0) { kd += k1d ; d += h1d; } else { kd += k2d; d += h2d; }

        p1->pix[t][z][y][x].a = a; p1->pix[t][z][y][x].b = b;
        p1->pix[t][z][y][x].c = c; p1->pix[t][z][y][x].d = d;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || a != aQ || b != bQ || c != cQ || d != dQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);

         fprintf (stderr, "aP = %d  aQ = %d  a = %d  da = %d  sa = %d  ma = %d  k1a = %d  k2a = %d  ka = %d\n",
                aP, aQ, a, da, sa, ma, k1a, k2a, ka);
        fprintf (stderr, "bP = %d  bQ = %d  b = %d  db = %d  sb = %d  mb = %d  k1b = %d  k2b = %d  kb = %d\n",
                bP, bQ, b, db, sb, mb, k1b, k2b, kb);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);
        fprintf (stderr, "dP = %d  dQ = %d  d = %d  dd = %d  sd = %d  md = %d  k1d = %d  k2d = %d  kd = %d\n",
                dP, dQ, d, dd, sd, md, k1d, k2d, kd);
   }
}


void NpicDrawLine_5q (Npic_image *np, 
        int sP, int tP, int zP, int yP, int xP, 
        int sQ, int tQ, int zQ, int yQ, int xQ, 
        int aP, int bP, int cP, int dP, int aQ, int bQ, int cQ, int dQ)
{
    Npic_image_5q *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        jm, j;

    /* We need signed integers */
    int a, da = abs (aQ-aP), sa = (aQ >= aP) ? 1:-1, ma, ka, k1a, k2a, h1a, h2a,
        b, db = abs (bQ-bP), sb = (bQ >= bP) ? 1:-1, mb, kb, k1b, k2b, h1b, h2b,
        c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c,
        d, dd = abs (dQ-dP), sd = (dQ >= dP) ? 1:-1, md, kd, k1d, k2d, h1d, h2d;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_5Q) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; 
    a = aP; b = bP; c = cP; d = dP; 
    p1->pix[s][t][z][y][x].a = a; p1->pix[s][t][z][y][x].b = b;
    p1->pix[s][t][z][y][x].c = c; p1->pix[s][t][z][y][x].d = d;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    ma = da / jm; k1a = 2*da - 2*ma*jm; ka = k1a - jm; k2a = ka - jm;
    mb = db / jm; k1b = 2*db - 2*mb*jm; kb = k1b - jm; k2b = kb - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;
    md = dd / jm; k1d = 2*dd - 2*md*jm; kd = k1d - jm; k2d = kd - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1a = ma*sa ; h2a = (ma+1)*sa;
    h1b = mb*sb ; h2b = (mb+1)*sb;
    h1c = mc*sc ; h2c = (mc+1)*sc;
    h1d = md*sd ; h2d = (md+1)*sd;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (ka < 0) { ka += k1a ; a += h1a; } else { ka += k2a; a += h2a; }
        if (kb < 0) { kb += k1b ; b += h1b; } else { kb += k2b; b += h2b; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }
        if (kd < 0) { kd += k1d ; d += h1d; } else { kd += k2d; d += h2d; }

        p1->pix[s][t][z][y][x].a = a; p1->pix[s][t][z][y][x].b = b;
        p1->pix[s][t][z][y][x].c = c; p1->pix[s][t][z][y][x].d = d;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || a != aQ || b != bQ || c != cQ || d != dQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);

         fprintf (stderr, "aP = %d  aQ = %d  a = %d  da = %d  sa = %d  ma = %d  k1a = %d  k2a = %d  ka = %d\n",
                aP, aQ, a, da, sa, ma, k1a, k2a, ka);
        fprintf (stderr, "bP = %d  bQ = %d  b = %d  db = %d  sb = %d  mb = %d  k1b = %d  k2b = %d  kb = %d\n",
                bP, bQ, b, db, sb, mb, k1b, k2b, kb);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);
        fprintf (stderr, "dP = %d  dQ = %d  d = %d  dd = %d  sd = %d  md = %d  k1d = %d  k2d = %d  kd = %d\n",
                dP, dQ, d, dd, sd, md, k1d, k2d, kd);
   }
}


void NpicDrawLine_6q (Npic_image *np, 
        int rP, int sP, int tP, int zP, int yP, int xP, 
        int rQ, int sQ, int tQ, int zQ, int yQ, int xQ, 
        int aP, int bP, int cP, int dP, int aQ, int bQ, int cQ, int dQ)
{
    Npic_image_6q *p1 = NpicCastImage(np);
    int x, dx = abs (xQ-xP), sx = (xQ >= xP) ? 1:-1, mx, kx, k1x, k2x, h1x, h2x,
        y, dy = abs (yQ-yP), sy = (yQ >= yP) ? 1:-1, my, ky, k1y, k2y, h1y, h2y,
        z, dz = abs (zQ-zP), sz = (zQ >= zP) ? 1:-1, mz, kz, k1z, k2z, h1z, h2z,
        t, dt = abs (tQ-tP), st = (tQ >= tP) ? 1:-1, mt, kt, k1t, k2t, h1t, h2t,
        s, ds = abs (sQ-sP), ss = (sQ >= sP) ? 1:-1, ms, ks, k1s, k2s, h1s, h2s,
        r, dr = abs (rQ-rP), sr = (rQ >= rP) ? 1:-1, mr, kr, k1r, k2r, h1r, h2r,
        jm, j;

    /* We need signed integers */
    int a, da = abs (aQ-aP), sa = (aQ >= aP) ? 1:-1, ma, ka, k1a, k2a, h1a, h2a,
        b, db = abs (bQ-bP), sb = (bQ >= bP) ? 1:-1, mb, kb, k1b, k2b, h1b, h2b,
        c, dc = abs (cQ-cP), sc = (cQ >= cP) ? 1:-1, mc, kc, k1c, k2c, h1c, h2c,
        d, dd = abs (dQ-dP), sd = (dQ >= dP) ? 1:-1, md, kd, k1d, k2d, h1d, h2d;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return;
    
    if (np->type != NPIC_IMAGE_6Q) {
        np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
        return;
    }

    if ( xP < 0 || xP >= p1->xmax || xQ < 0 || xQ >= p1->xmax
      || yP < 0 || yP >= p1->ymax || yQ < 0 || yQ >= p1->ymax
      || zP < 0 || zP >= p1->zmax || zQ < 0 || zQ >= p1->zmax
      || tP < 0 || tP >= p1->tmax || tQ < 0 || tQ >= p1->tmax
      || sP < 0 || sP >= p1->smax || sQ < 0 || sQ >= p1->smax
      || rP < 0 || rP >= p1->rmax || rQ < 0 || rQ >= p1->rmax
    ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": extrem points outside image");
        return;
    }

    /* Draw first point */
    x = xP; y = yP; z = zP; t = tP; s = sP; r = rP; 
    a = aP; b = bP; c = cP; d = dP; 
    p1->pix[r][s][t][z][y][x].a = a; p1->pix[r][s][t][z][y][x].b = b;
    p1->pix[r][s][t][z][y][x].c = c; p1->pix[r][s][t][z][y][x].d = d;

    /* Search jm = max (dx, dy, ..) but not on dc */
    jm = dx;
    if (jm < dy) jm = dy;
    if (jm < dz) jm = dz;
    if (jm < dt) jm = dt;
    if (jm < ds) jm = ds;
    if (jm < dr) jm = dr;

    if (jm == 0) return;  /* P == Q */

    /* Init control variables */
    mx = dx / jm; k1x = 2*dx - 2*mx*jm; kx = k1x - jm; k2x = kx - jm;
    my = dy / jm; k1y = 2*dy - 2*my*jm; ky = k1y - jm; k2y = ky - jm;
    mz = dz / jm; k1z = 2*dz - 2*mz*jm; kz = k1z - jm; k2z = kz - jm;
    mt = dt / jm; k1t = 2*dt - 2*mt*jm; kt = k1t - jm; k2t = kt - jm;
    ms = ds / jm; k1s = 2*ds - 2*ms*jm; ks = k1s - jm; k2s = ks - jm;
    mr = dr / jm; k1r = 2*dr - 2*mr*jm; kr = k1r - jm; k2r = kr - jm;
    ma = da / jm; k1a = 2*da - 2*ma*jm; ka = k1a - jm; k2a = ka - jm;
    mb = db / jm; k1b = 2*db - 2*mb*jm; kb = k1b - jm; k2b = kb - jm;
    mc = dc / jm; k1c = 2*dc - 2*mc*jm; kc = k1c - jm; k2c = kc - jm;
    md = dd / jm; k1d = 2*dd - 2*md*jm; kd = k1d - jm; k2d = kd - jm;

    /* Init coordinates increments */
    h1x = mx*sx ; h2x = (mx+1)*sx;
    h1y = my*sy ; h2y = (my+1)*sy;
    h1z = mz*sz ; h2z = (mz+1)*sz;
    h1t = mt*st ; h2t = (mt+1)*st;
    h1s = ms*ss ; h2s = (ms+1)*ss;
    h1r = mr*sr ; h2r = (mr+1)*sr;
    h1a = ma*sa ; h2a = (ma+1)*sa;
    h1b = mb*sb ; h2b = (mb+1)*sb;
    h1c = mc*sc ; h2c = (mc+1)*sc;
    h1d = md*sd ; h2d = (md+1)*sd;

    /* Draw jm next points */
    for (j = 1; j <= jm; j++) {
        if (kx < 0) { kx += k1x ; x += h1x; } else { kx += k2x; x += h2x; }
        if (ky < 0) { ky += k1y ; y += h1y; } else { ky += k2y; y += h2y; }
        if (kz < 0) { kz += k1z ; z += h1z; } else { kz += k2z; z += h2z; }
        if (kt < 0) { kt += k1t ; t += h1t; } else { kt += k2t; t += h2t; }
        if (ks < 0) { ks += k1s ; s += h1s; } else { ks += k2s; s += h2s; }
        if (kr < 0) { kr += k1r ; r += h1r; } else { kr += k2r; r += h2r; }
        if (ka < 0) { ka += k1a ; a += h1a; } else { ka += k2a; a += h2a; }
        if (kb < 0) { kb += k1b ; b += h1b; } else { kb += k2b; b += h2b; }
        if (kc < 0) { kc += k1c ; c += h1c; } else { kc += k2c; c += h2c; }
        if (kd < 0) { kd += k1d ; d += h1d; } else { kd += k2d; d += h2d; }

        p1->pix[r][s][t][z][y][x].a = a; p1->pix[r][s][t][z][y][x].b = b;
        p1->pix[r][s][t][z][y][x].c = c; p1->pix[r][s][t][z][y][x].d = d;
    }

    /* Debug code */
    if (x != xQ || y != yQ || z != zQ || t != tQ || s != sQ || r != rQ || a != aQ || b != bQ || c != cQ || d != dQ) {
        np->gen.ok = NpicError (__func__, NPIC_ERROR, ": dest point not reached");

        fprintf (stderr, "xP = %d  xQ = %d  x = %d  dx = %d  sx = %d  mx = %d  k1x = %d  k2x = %d  kx = %d\n",
                xP, xQ, x, dx, sx, mx, k1x, k2x, kx);
        fprintf (stderr, "yP = %d  yQ = %d  y = %d  dy = %d  sy = %d  my = %d  k1y = %d  k2y = %d  ky = %d\n",
                yP, yQ, y, dy, sy, my, k1y, k2y, ky);
        fprintf (stderr, "zP = %d  zQ = %d  z = %d  dz = %d  sz = %d  mz = %d  k1z = %d  k2z = %d  kz = %d\n",
                zP, zQ, z, dz, sz, mz, k1z, k2z, kz);
        fprintf (stderr, "tP = %d  tQ = %d  t = %d  dt = %d  st = %d  mt = %d  k1t = %d  k2t = %d\n  kt = %d",
                tP, tQ, t, dt, st, mt, k1t, k2t, kt);
        fprintf (stderr, "sP = %d  sQ = %d  s = %d  ds = %d  ss = %d  ms = %d  k1s = %d  k2s = %d\n  ks = %d",
                sP, sQ, s, ds, ss, ms, k1s, k2s, ks);
        fprintf (stderr, "rP = %d  rQ = %d  r = %d  dr = %d  sr = %d  mr = %d  k1r = %d  k2r = %d\n  kr = %d",
                rP, rQ, r, dr, sr, mr, k1r, k2r, kr);

         fprintf (stderr, "aP = %d  aQ = %d  a = %d  da = %d  sa = %d  ma = %d  k1a = %d  k2a = %d  ka = %d\n",
                aP, aQ, a, da, sa, ma, k1a, k2a, ka);
        fprintf (stderr, "bP = %d  bQ = %d  b = %d  db = %d  sb = %d  mb = %d  k1b = %d  k2b = %d  kb = %d\n",
                bP, bQ, b, db, sb, mb, k1b, k2b, kb);
        fprintf (stderr, "cP = %d  cQ = %d  c = %d  dc = %d  sc = %d  mc = %d  k1c = %d  k2c = %d  kc = %d\n",
                cP, cQ, c, dc, sc, mc, k1c, k2c, kc);
        fprintf (stderr, "dP = %d  dQ = %d  d = %d  dd = %d  sd = %d  md = %d  k1d = %d  k2d = %d  kd = %d\n",
                dP, dQ, d, dd, sd, md, k1d, k2d, kd);
   }
}


/*----------------------------------------------------------------------------*/