Rotating a 3D body in python results in holes in the body

Question

I have a 3D body, let's say a cuboid, which I want to rotate. For simplicity, let's assume it's just rotation around the x-axis. So I use the corresponding rotation matrix R and multiply it with the coordinate vector v to get the new coordinate vector v': v'=R*v. As a visualization tool, I use mayavi.

While the rotation does work, it has some annoying side effect: some values inside of the rotated cuboid are missing. You can see that in the following snapshot, the blue cuboid is the rotated body, and it has some "holes" in it.

My question now is, what am I doing wrong with my rotation?

Here is the corresponding python code:

# import standard modules
import numpy as np

# import modules for 3D visualization of data
from mayavi import mlab

def plot_simple( data2plot ):

    # define contour levels
    contLevels  = np.linspace(0, np.amax(data2plot), 5)[1:].tolist()

    # create figure with white background and black axes and labels
    fig1    = mlab.figure( bgcolor=(1,1,1), fgcolor=(0,0,0), size=(800,600) )

    # make the contour plot
    cont_plt    = mlab.contour3d( data2plot, contours=contLevels,
                                  transparent=True, opacity=.4,
                                  figure=fig1 )

    # create axes instance to modify some of its properties
    ax1 = mlab.axes( nb_labels=4, extent=[1, data2plot.shape[0],
                                          1, data2plot.shape[1],
                                          1, data2plot.shape[2] ] )
    mlab.outline(ax1)
    ax1.axes.label_format   = '%.0f'
    ax1.axes.x_label    = 'x'
    ax1.axes.y_label    = 'y'
    ax1.axes.z_label    = 'z'

    # set initial viewing angle
    mlab.view( azimuth=290, elevation=80 )

    mlab.show()


def Rx(alpha):
    # rotation matrix for rotation around x-axis
    return np.matrix([[ 1, 0            , 0             ],
                      [ 0, np.cos(alpha), -np.sin(alpha)],
                      [ 0, np.sin(alpha), np.cos(alpha) ]])


def make_rotated_cube( Nx=100, Ny=70, Nz=40 ):
    arr = np.zeros( [Nx, Ny, Nz] )

    # define center of cuboid
    xc  = Nx/2
    yc  = Ny/2
    zc  = Nz/2

    # define width of cuboid in each direction
    dx  = Nx/4
    dy  = Ny/4
    dz  = Nz/4

    # rotation angle in degrees
    alpha   = 20
    alpha   = np.radians(alpha)

    # loop through arr and define cuboid and rotated cuboid
    # note that this is a very inefficient way to define a cuboid
    # (the actual thing to rotate is different, this is just to make it simple)
    for ii in range(Nx):
        for jj in range(Ny):
            for kk in range(Nz):
                # check if coordinate is inside original cuboid
                if (    (ii > (xc-dx/2) and ii < (xc+dx/2))
                    and (jj > (yc-dy/2) and jj < (yc+dy/2))
                    and (kk > (zc-dz/2) and kk < (zc+dz/2)) ):
                    # set density of original cuboid
                    arr[ii,jj,kk]   = 5.

                    # apply rotation
                    new_coords  = Rx(alpha)*np.array([[ii],[jj],[kk]])
                    # set density of rotated cuboid to different value
                    arr[ round(new_coords[0,0]),
                         round(new_coords[1,0]),
                         round(new_coords[2,0]) ] = 2

    return arr


def main():
    cubes   = make_rotated_cube()
    plot_simple(cubes)


if __name__ == '__main__':
    main()

Answer 1

Thanks to @Julien and his comment, I tested the skewing method. As can be seen in the following figures, the holes are indeed gone! One can now clearly see, however, that the edges are quite ugly, but that should be due to the bad resolution of the grid in my understanding (as I am not "loosing" any pixels with this method as I was before).

I just changed the code in the for loop (following the comment # apply rotation):

# translations
skew_1  = np.tan(alpha/2)
skew_2  = -np.sin(alpha)
skew_3  = skew_1
# calculating new coords
kk_new   = kk + round(skew_1*jj)
jj_new   = jj + round(skew_2*kk_new)
kk_new   = kk_new + round(skew_1*jj_new)
arr[ ii, jj_new, kk_new ] = 2

I would still like to see a solution involving an interpolation (but that skewing method is pretty cool).