Getting Isosurfaces from a data grid into POV-Ray

POV-Ray makes for fantastic scientific visualizations but if you want to render isosurfaces from a 3D scalar field, there isn’t a very good option. In the past, I used a meta-ball method to make smooth membranes for isosurfaces. There is a isosurface object that allows you to specify a mathematical function that is evaluated to compute an isosurface, but if your data is just a 3D grid, than this isn’t a very good option. Or is it? My new method makes a Fourier series plane wave reconstruction of the data using an FFT. Since plane waves can be evaluated as function, POV-Ray is effectively doing Fourier interpolation! In this scheme, the discrete Fourier transform value at each k-point is a complex value where the real component is the amplitude of the cosine component and the imaginary component is the amplitude of the sine component. You can implement this easily but immediately there will be a problem for any data grid with more than about 10,000 plane waves. POV-Ray won’t parse it, returning the error: ‘Parse Error: Function too complex!’. Since we the isosurface object is only designed real valued functions, we can use the Hermitian symmetry to effectively increase this upper limit to about 20,000. In other words, because the data is real valued, forwards plane waves and backwards plane waves are the same so the backwards ones can be ignored if the forwards ones have their amplitudes doubled. In my code, any plane wave with a negative x component in k-space is ignored and any with a positive x component is amplitude doubled. We can double the plane wave limit again to 40,000 by merging the sine and cosine components into one cosine wave with a phase shift! Quadrupling our capability is great, but 40,000 plane waves still only represents the FFT of a cube of data 34 cell in edge length. The real trick that makes this method possible is sorting by plane wave amplitude and taking only the strongest amplitudes. For my example system, I am using the electron density of a subphthalocyanine crystal. Below I am showing the plane wave amplitudes sorted in order of decreasing amplitude. I am also showing the cumulative distribution function of these amplitudes, while not directly meaningful since waves have phase too, it helps demonstrate how complete a plane wave reconstruction is for a selection of the strongest waves. What does this look like? For 100 plane waves: For 1000 plane waves For 10000 plane waves:
I think that looks nice enough for publication.