Rendering voxel functions in Blender

Hi all, I need to be able to render arbitrary functions easily. In other words, the following function would draw a sphere by checking a giant 3D grid of say… 100x100x100 voxels to check if the point is in the sphere:

// Check through multiple x, y and z voxels to draw a sphere:
bool sphere(float x, float y, float z) {    // x, y and z are normalized to 1 or less.
 float f = sqrt(x^2 + y^2 + z^2)
 if (f<1.0) return 1;	// Part of sphere
 else return 0;		// Not part of sphere

Is this easy to do in Blender or maybe some other 3D software? An equivalent to voxels would be to trace rays to see if it intercepts the function. Either way is fine as long as a sphere is produced from the above ‘code’/function.

Blender’s fluid simulation uses Voxels, but other than that Blender as of 2.49b does not have access to voxels. You can create a 3d object in Blender, export it out and into Voxel3d.
I’m not sure if Voxel3d will allow you to use calculations or not. It is here:

There was talk of 2.5 using volumetrics. There is some versions of Blender that you may find that use volumetrics at I have not used those though, but they are there to explore.



Okay, even if it’s not done directly through voxels, it’s still possible in theory to trace rays to see if they intercept the function’s criteria. Is this possible to do - create a function like the above and have Blender render it? (preferably through Indigo or similar for max realism).

I’m not sure about Blender rendering a function… Well, as long as you’re going for a sphere thre’s no problem, just add a sphere, Subsurf, set smooth etc. If you wanted a more complex function, the way would be to make a python script working like yours, checking through a grid of verts and creating them where the function works, then keeping the ones on the limits, skin them, render, and you should see the result of the function. I’m not a good python scripter, but still this could work, couldn’t it? Hope it helped…

what are you trying to do? Its very easy to render a 3d grid of volume in a object even in 2.49b. Just add a particle system to the object Suzanne for ex. then change the settings to emit from volume and choose grid and up the res. then chose the emit to object and change the particle size to fit.

Thanks for the replies. I would like to render arbitrary 3D fractals, so a sphere would be too simple for my needs.

What approach do you think would be easier / more flexible out of both of yours? I am a beginner, so any tutorial on the net regarding this kind of thing would be great.

Are isosurfaces planned for a future version of Blender?

sorry I not sure what you mean by 3d fractals… I know what fractals but I am not sure what you are trying to create. Arrays modifer might help?? there is another program too that is open source

I’m also interested in rendering functions as densities in 3 dimensions (4D rendering, would have numbers between 1 and 0 too). I was thinking about rendering a quantum mechanical wave function as a volumetric object, or at least the probability density. It’s not too hard to generate the data using python or matlab, but I’m not sure how to load it into the volumetrics of 2.5. It lists ‘image sequence’ as a type, but nothing comes up to choose one. Would be nice to have function support built in…

I’ll write the script to do it if I can find out how to load it as voxels.

Actual isosurfaces would be trickier…


Here you can find is a python script by Matt Ebb which might be of interest?

With the newer 2.5 build I cannot get Voxel Data image sequences to load exactly
as HenryIII describes, but you can use Image sequences and 8-bit Raw with the
following build on

Here is a screencast by Matt Ebb showing how to setup using an image sequence:

Hope this Helps?

Kindest Regards,


Hi Again,

Almost forgot a link to some datasets:



Thank you, pixeltwist! I knew about some of those, but I didn’t know about that last link. Raw data does work in 2.5, so I assume that image sequence is just missing some operators at the moment. I’m looking at those raw files to see if I can output from python to that format, otherwise I’ll have to try that build. Oh, and thank you for letting me know about that build, if I can’t use the RAW that should work great!

Hi HenryIII,

Glad to be of help.

Rgds P.

Okay, here’s my code. I managed to use the raw format, so either an older build or a current one will render it. I tested it with sin(x) + sin(y) + sin(z), and with spherical functions. To use, just take the output file (I keep all of my scripts in C:\Python31, so my file ended up there) and load it as an 8 bit raw, and set the three resolution parameters to match ‘res’. Play with the density till it’s pretty. :slight_smile:

Oh, and I ran it in IDLE, but I would suppose that blender could do it, too.

To get a isosurface, just play with the color ramp. It should work, depending on the resolutions.

# 4D Function Plotter
# For use with Blender Volumetrics (2.5 prerelease version)
# Henry Schreiner
# Python 3.1

from math import *

class mat4d:
    '''Make a 4D matrix with a given function'''
    def __init__(self,
                 func: 'Has to have 3 variables'=lambda x,y,z: x+y+z,
                 res=9, rang=5, sph=False):
        '''This uses x,y,z; or it uses r, phi, and theta if shp is True'''


        self.M = [0]*res**3
        ipos = lambda a: (a/(res-1)-.5)*2*rang

        for n in range(res**3):
            if sph: (x,y,z)=self.rec2sph(x,y,z)
            try: self.M[n]=func(x,y,z)
            except ZeroDivisionError: self.M[n]=0

    def save(self, filename):
        '''This will save a .raw.gz file for Blender RAW volumetrics'''
        fob = open(filename + '.raw.gz','bw')

    def num2ind(self,a):
        return (int(a/self.res**2),

    def ind2num(self,i,j,k):
       return i*self.res**2 + j*self.res + k

    def normalize(self, absolute=False, scale =255):
        '''Normalizes to 0-scl, uses absolute values or total'''
        minimum = min(map(abs,self.M) if absolute else self.M)
        rang = max(map(abs,self.M) if absolute else self.M) - minimum
        for n in range(self.res**3):
            self.M[n]=((abs(self.M[n]) if absolute else self.M[n])

    def rec2sph(x,y,z):
        '''Returns (r, phi, theta) for x,y,x'''
        r= sqrt(x**2 + y**2 + z**2)
        p = atan2(y,x)
        t = (acos(z/r) if r else 0.0) # quick and dirty fix for r=0
        return (r, p, t)

## Example usage
# Setting up wavefunction arrays
psi = [lambda r: 0,                                           #n=0
       lambda r: 1/(r*sqrt(2*pi))*sin(pi*r),                  #n=1
       lambda r: 1/(r*sqrt(2*pi))*sin(2*pi*r)]                #n=2

Y = [[lambda p,t: (1/(4*pi))**(1/2)],                         #lm=00
     [lambda p,t: (3/(4*pi))**(1/2)*cos(t),                   #lm=10
      lambda p,t: -(3/(8*pi))**(1/2)*sin(t)*e**(p*1j)],       #lm=11
     [lambda p,t: (5/(16*pi))**(1/2)*(3*cos(t)**2-1),         #lm=20
      lambda p,t: -(15/(8*pi))**(1/2)*sin(t)*cos(t)*e**(p*1j),#lm=21
      lambda p,t: (15/(32*pi))**(1/2)*sin(t)**2*e**(p*2j)]]   #lm=22

# Example using nlm=222
def wavefn(r,p,t):
    if r>1: return 0.0
    else: return (Y[l][m](p,t)*Y[l][m](p,t).conjugate()).real * psi[n](r)**2

if __name__ == '__main__':
    import time
    t = time.time()

    #  !!! This is how you use it!!!
    q=mat4d(wavefn,res=32,rang=1,sph=True) # 1st parameter = function
    q.normalize() # scales numbers to 0-255, for 8 bit use'wvfn222') # saves it as the name 
    t = time.time()-t
    print('The time taken was:', round(t,3), 'seconds.')


Don’t know if you saw this: