3D slices using matplotlib and python.

Share on FacebookShare on Google+Tweet about this on TwitterShare on LinkedInPin on PinterestShare on StumbleUponShare on Tumblr

This post is based on a figure I needed to make for a publication. I had to make a figure outlining a flow of data in an analysis. Data was extracted from PET images and was analyzed parallelly. I wanted to make a detailed enough figure that shows how the data is analyzed so I used ‘matplotlib’. The final figure had several different stages, but all stages were based on this image or the components of this image. So, if someone want to draw something similar, you can get help from the code I am posting.

Also, this post contains a python method directly taken from stack overflow. Thanks to the original author.
http://stackoverflow.com/a/31364297/480551

And of course, inputs/better implementations are welcome.

result

import numpy as np
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D

def set_axes_equal(ax):
    #Taken from http://stackoverflow.com/a/31364297/480551
    '''Make axes of 3D plot have equal scale so that spheres appear as spheres,
    cubes as cubes, etc..  This is one possible solution to Matplotlib's
    ax.set_aspect('equal') and ax.axis('equal') not working for 3D.

    Input
      ax: a matplotlib axis, e.g., as output from plt.gca().
    '''

    x_limits = ax.get_xlim3d()
    y_limits = ax.get_ylim3d()
    z_limits = ax.get_zlim3d()

    x_range = x_limits[1] - x_limits[0]; x_mean = np.mean(x_limits)
    y_range = y_limits[1] - y_limits[0]; y_mean = np.mean(y_limits)
    z_range = z_limits[1] - z_limits[0]; z_mean = np.mean(z_limits)

    # The plot bounding box is a sphere in the sense of the infinity
    # norm, hence I call half the max range the plot radius.
    plot_radius = 0.5*max([x_range, y_range, z_range])

    ax.set_xlim3d([x_mean - plot_radius, x_mean + plot_radius])
    ax.set_ylim3d([y_mean - plot_radius, y_mean + plot_radius])
    ax.set_zlim3d([z_mean - plot_radius, z_mean + plot_radius])

def getFaceColors(colorsk, col):
    '''This method return the color matrix for the surface. '''
    (sx, sy, sc) = colorsk.shape
    maxGridSize_d = sx
    maxGridSize_h = sy
    for i in range(sx):
        for j in range(sy):
            total = np.floor((np.floor(maxGridSize_h*1))-abs(np.floor(maxGridSize_h*0.6)-i)*maxGridSize_d/maxGridSize_h)
            firstidx = np.floor((maxGridSize_h-total)/2)
            lastidx = firstidx+total

            if firstidx < j & j < lastidx:
                if col == 'r':
                    colorsk[i, j, :] = np.array([np.random.random(size=1), 0, 0, 1], dtype=np.float64)
                if col == 'g':
                    colorsk[i, j, :] = np.array([0, np.random.random(size=1), 0, 1], dtype=np.float64)
                if col == 'b':
                    colorsk[i, j, :] = np.array([0, 0, np.random.random(size=1), 1], dtype=np.float64)

            else:
                colorsk[i, j, :] = np.array([1, 1, 1, 1], dtype=np.float64)
    return colorsk

def drawImage(xsize, ysize, zsize):

    '''This method draws the main figure. This method draws the meshes one
    adjusting location'''

    fig = pyplot.figure(figsize=(5, 5))
    ax = fig.add_subplot(111, projection='3d')

    u = np.linspace(0, xsize, xsize + 1)
    v = np.linspace(0, ysize, ysize + 1)

    x, y = np.meshgrid(u, v)
    z = np.ones_like(x)
    
    # This section is commented as I only need x surfaces. But if someone need y and z surfaces, uncomment. 
    # for i in range(zsize + 1):
    #     u = np.linspace(0, xsize, xsize + 1)
    #     v = np.linspace(0, ysize, ysize + 1)
    #     x, y = np.meshgrid(u, v)
    #     z = np.ones((ysize + 1, xsize + 1), dtype=np.int)
    #
    #     # Draw surface on this grid. Included if needed to draw a structur with cubes.
    #     surf = ax.plot_surface(x, y, z*i, rstride=1, cstride=1, color='w', alpha=0.7, shade=False)
    #
    # for i in range(ysize + 1):
    #     u = np.linspace(0, xsize, xsize + 1)
    #     v = np.linspace(0, zsize, zsize + 1)
    #     x, z = np.meshgrid(u, v)
    #     y = np.ones((zsize + 1, xsize + 1), dtype=np.int)
    #
    #     # Draw surface on this grid. Included if needed to draw a structur with cubes.
    #     surf = ax.plot_surface(x, y*i, z, rstride=1, cstride=1, color='w', alpha=0.7, shade=False)

    for i in range(xsize + 1):
        u = np.linspace(0, zsize, zsize + 1)
        v = np.linspace(0, ysize, ysize + 1)
        z, y = np.meshgrid(u, v)
        x = np.ones((ysize + 1, zsize + 1), dtype=np.int)
        colors = np.ones((ysize, zsize, 4), np.float64)

        # Draw surface on the grid.s
        surf = ax.plot_surface(x * i, y, z, rstride=1, cstride=1, facecolors=getFaceColors(colors, 'r'), alpha=0.7,
                               shade=False)
        surf.set_edgecolor('k')
        surf = ax.plot_surface(x * i + (xsize + 5), y, z, rstride=1, cstride=1, facecolors=getFaceColors(colors, 'g'),
                               alpha=0.7, shade=False)
        surf.set_edgecolor('k')
        surf = ax.plot_surface(x * i + 2 * (xsize + 5), y, z, rstride=1, cstride=1,
                               facecolors=getFaceColors(colors, 'b'), alpha=0.7, shade=False)
        surf.set_edgecolor('k')

    ax.axis('off')
    set_axes_equal(ax)
    ax.view_init(elev=20, azim=-45)

    # Modify these values for the align the plot on the figure.
    pyplot.subplots_adjust(left=0, right=1.2, top=1, bottom=-0.05)
    pyplot.savefig('all_imagesp_2.png', dpi=600, pad_inches=-1, bbox_inches='tight')


if __name__ == '__main__':
    drawImage(10, 25, 20)

Changing the colouring method, you can achieve multiple results. Changing the colouring method to

def getFaceColors (colorsk, col):
    '''This method return the color matrix for the surface. '''
    (sx, sy, sc) = colorsk.shape
    maxGridSize_d = sx
    maxGridSize_h = sy
    for i in range(sx):
        for j in range(sy):
            total = np.floor((np.floor(maxGridSize_h*1))-abs(np.floor(maxGridSize_h*0.6)-i)*maxGridSize_d/maxGridSize_h)
            firstidx = np.floor((maxGridSize_h-total)/2)
            lastidx = firstidx+total

            if col == 'r':
                colorsk[i, j, :] = np.array([np.random.random(size=1), 0, 0, 1], dtype=np.float64)
            if col == 'g':
                colorsk[i, j, :] = np.array([0, np.random.random(size=1), 0, 1], dtype=np.float64)
            if col == 'b':
                colorsk[i, j, :] = np.array([0, 0, np.random.random(size=1), 1], dtype=np.float64)
    return colorsk

will result in result_2

Hope this helps.

Share on FacebookShare on Google+Tweet about this on TwitterShare on LinkedInPin on PinterestShare on StumbleUponShare on Tumblr

Leave a Reply

Your email address will not be published. Required fields are marked *