"Draw" a 3d line into an array

Question

I am looking for a solution do "draw" a 3d line (filled cylinder) into a 3d array using python just like the skimage.draw.line function does for 2 arrays.

The line should have a starting point (x1, y1, z1), an end point (x2, y2, z2), and a radius R.

I had a look at an example for a 2d line but I was not able to modify it to work in the 3d case.

I was thinking about drawing successive ellipses into the array, but I was not able to figure out how to calculate the two axes and the rotation angle.

Maybe there is a much simpler approach to this problem?

Answer 1

Let assume this GLSL volumetric back raytracer as a start point. To make a filled 3D line like this:

You need:

1. endpoints as spheres

   see the add_sphere in the link above.

2. discs cut at the endpoints

   for that we need U,V basis vectors (perpendicular vectors to each other and to line itself). With those we can simply use any 2D circle pixels acquisition and convert them to 3D voxel positions with ease. So if u,v are coordinates in some 2D circle centered at (0,0) then:

   x = x0 + u*U.x + v*V.x
   y = y0 + u*U.y + v*V.y
   z = z0 + u*U.z + v*V.z
   

   (x,y,z) are corresponding to 3D circle voxel coordinate with center (x0,y0,z0). For more info see my C++ glCircle3D implementation.

3. line body

   As wee got all the voxel positions in the disc around x0,y0,z0 endpoint of the line just cast a line from each of it with the same slope as line (x0,y0,z0),(x1,y1,z1) which are the endpoints of your line.

When put together in C++ (sorry I do not code in python) I got this:

void volume::add_line(int x0,int y0,int z0,int x1,int y1,int z1,int r,GLuint col)
    {
    if (!_init) return;
    int i,n,x,y,z,cx,cy,cz,dx,dy,dz,kx,ky,kz;
    // endpoints are (half)spheres
    add_sphere(x0,y0,z0,r,col);
    add_sphere(x1,y1,z1,r,col);
    // DDA constants
    kx=0; dx=x1-x0; if (dx>0) kx=+1; if (dx<0) { kx=-1; dx=-dx; } dx++;           n=dx;
    ky=0; dy=y1-y0; if (dy>0) ky=+1; if (dy<0) { ky=-1; dy=-dy; } dy++; if (n<dy) n=dy;
    kz=0; dz=z1-z0; if (dz>0) kz=+1; if (dz<0) { kz=-1; dz=-dz; } dz++; if (n<dz) n=dz;
    // basis vectors
    double U[3],V[3],N[3]={x1-x0,y1-y0,z1-z0},u,v,rr=r*r;
    vector_one(N,N); // unit vector
    vector_ld(U,1.0,0.0,0.0); if (fabs(vector_mul(U,N))>=0.75) vector_ld(U,0.0,1.0,0.0); // |dot(U,N)|<0.75 means (1.0,0.0,0.0) is nearly parallel to N so chose (0.0,1.0,0.0) instead
    vector_mul(U,U,N); // U = U x N
    vector_mul(V,U,N); // V = U x N
    vector_one(U,U);   // U /= |U|
    vector_one(V,V);   // V /= |V|
    // disc
    for (u=-r;u<=+r;u++)
     for (v=-r;v<=+r;v++)
      if (u*u+v*v<=rr)
        {
        x=x0+double((u*U[0])+(v*V[0]));
        y=y0+double((u*U[1])+(v*V[1]));
        z=z0+double((u*U[2])+(v*V[2]));
        // DDA line
        for (cx=cy=cz=n,i=0;i<n;i++)
            {
            if ((x>=0)&&(x<size)&&(y>=0)&&(y<size)&&(z>=0)&&(z<size)) data[z][y][x]=col;
            cx-=dx; if (cx<=0) { cx+=n; x+=kx; }
            cy-=dy; if (cy<=0) { cy+=n; y+=ky; }
            cz-=dz; if (cz<=0) { cz+=n; z+=kz; }
            }
        }
    }

The vector_xxx functions are just my 3D vector math and just dot,cross product and normalize to unit size are used which is easy to implement. You can see them here:

• Understanding 4x4 homogenous transform matrices

There are still things that can be improved like the spheres could be just half spheres and their generation can be joined with the disc stuff ... as the dot between normal and not offseted 3D sphere coordinate is either positive/zero/negative which distinct endpoint half-sphere and disc ... that would also fully eliminate the need for U,V.

Also depending on used HW and circumstances there might be also faster approaches like analytical (filling BBOX based on distance from line) if fast vector math is combined with massive parallelism like on GPU.

After some tweaking in my engine (added zoom and handle some accuracy problem) I got this result:

for 128x128x128 volume inited like this:

// init volume raytracer
vol.gl_init();
vol.beg();
int r,a,b,c;
r=10.0; a=r+1; b=vol.size-r-2; c=vol.size>>1;
                             //BBGGRR
vol.add_line(a,a,a,b,a,a,r,0x00FF2020);
vol.add_line(a,b,a,b,b,a,r,0x00FF2020);
vol.add_line(a,a,a,a,b,a,r,0x00FF2020);
vol.add_line(b,a,a,b,b,a,r,0x00FF2020);

vol.add_line(a,a,b,b,a,b,r,0x00FF2020);
vol.add_line(a,b,b,b,b,b,r,0x00FF2020);
vol.add_line(a,a,b,a,b,b,r,0x00FF2020);
vol.add_line(b,a,b,b,b,b,r,0x00FF2020);

vol.add_line(a,a,a,a,a,b,r,0x00FF2020);
vol.add_line(a,b,a,a,b,b,r,0x00FF2020);
vol.add_line(b,a,a,b,a,b,r,0x00FF2020);
vol.add_line(b,b,a,b,b,b,r,0x00FF2020);

vol.add_sphere(c,c,c,c>>1,0x00FF8040);
vol.add_sphere(a,c,c,r,0x004080FF);
vol.add_sphere(b,c,c,r,0x0080FF40);
vol.add_sphere(c,a,c,r,0x00FF4080);
vol.add_sphere(c,b,c,r,0x00AAAAAA);
vol.add_box(c,c,a,r,r,r,0x0060FF60);
vol.add_box(c,c,b,r,r,r,0x00FF2020);

vol.end();

Answer 2

I finally managed to implement this in python using skimage.draw.ellipse:

import numpy as np
from numpy.linalg import norm
import skimage.draw

c = lambda *x: np.array(x, dtype=float)
vector_angle = lambda V, U: np.arccos(norm(np.dot(V, U)) / (norm(V) + norm(U)))

r = 10 # radius of cylinder
C0 = c(10, 10, 10) # first (x,y,z) point of cylinder
C1 = c(99, 90, 15) # second (x,y,z) point of cylinder

C = C1 - C0

X, Y, Z = np.eye(3)

theta = vector_angle(Z, C)
print('theta={} deg'.format(theta / np.pi * 180))

minor_axis = r
major_axis = r / np.cos(theta)
print('major_axis', major_axis)

alpha = vector_angle(X, C0 + C)
print('alpha={} deg'.format(alpha / np.pi * 180))

data = np.zeros([100, 100, 100])
nz, ny, nx = data.shape

for z in range(nz):
    lam = - (C0[2] - z)/C[2]
    P = C0 + C * lam
    y, x = skimage.draw.ellipse(P[1], P[0], major_axis, minor_axis, shape=(ny, nx), rotation=alpha)
    data[z, y, x] = 1