How-to use Blender's geometry nodes for computation on CMB
In this note I’m going to introduce Blender’s power to scientists, mainly cosmologists. Blender is an open-source cross-platform 3D content creation software and supports the entire 3D pipeline (see Blender’s about for more information).
We will see how to visialize cmb maps (see visualization nodes) and how do we can deal with cmb maps like a geometry and benefit from (implicit) gpu parallel computing. I will show you how to import your data into Blender and how to work with it step by step. The methods explained here are for local computers but could also be applied elsewhere (e.g. google colab). Each step could be a separate post but the order could be problematic in case of separation.
Creating HEALPix mesh
First, we need to create a 3d mesh out of healpix data. There are several ways to do it but the easiest way is to use healpy package.
Since at the moment, healpy is only available on linux & mac, windows users can use wsl(windows subsystem for linux) that allows them to run linux terminals in windows environment.
In the following, we will create a handmade OBJ file that can be imported to any 3d DCC. Now in this section we will use system’s python to create our mesh, not Blender’s, since Blender uses its own internal python which is a totally separate one (Insha’allah I will write a post on how to install external modules on Blender’s python; but even if we install healpy on Blender’s python it is not still useable in windows).
First, we create a function to generate texts needed for making OBJ file.
1
2
3
4
5
6
7
8
9
import numpy as np
import healpy as hp
def generate_obj_txt(element_arr, str_key):
txt = ""
space, next_line = " ", "\n"
for elem in element_arr:
txt += str_key + "".join(space + str(component) for component in elem) + next_line
return txt
The above function creates simple texts like this:
1
2
3
4
5
v 0.0 1.0 -0.0
v 0.0 1.0 1.0
f 1 2 3 4
f 5 6 7 8
...
That defines vertex coordinates and certain indices that make a face(forget about any additional data; edges, skinning data etc. since we don’t need those in this tutorial).
Then we store pixels’ corners(boundaries) in a well sorted way like bellow:
1
2
3
4
5
6
7
8
9
10
nside = 64
npix = 12 * nside ** 2
v_coords = np.zeros((npix * 4, 3)) # 4 corners for each pixel and 3 coords for each corner
for i in range(npix):
pbounds = hp.boundaries(nside = nside, pix = i, step = 1, nest = False) # we use ring indexing
v_coords[4 * i : 4 * (i + 1)] = np.transpose(pbounds)
v_indices = np.arange(len(v_coords)) + 1
faces = v_indices.reshape((npix, 4))
One can optionally merge duplicate vertices here, by utilizing merging algorithms, but it is much easier and faster to do this process inside Blender.
Now its time to write the OBJ file using above arrays of vertices and faces. Note that we have to change coordinates from (right handed)xyz to (left handed)xzy (y_lh = z , z_lh = -y) since its the most usual used case for OBJ file importers:
1
2
3
4
5
6
7
8
# mesh creation
f3dname = "healpix_sphere.obj"
with open(f3dname,'w') as f3d:
# before saving, change coords
v_coords[:, 1], v_coords[:, 2] = np.copy(v_coords[:, 2]), -np.copy(v_coords[:, 1])
ftxt = generate_obj_txt(v_coords, "v")
ftxt += generate_obj_txt(faces, "f")
f3d.write(ftxt)
For high resolution maps(nside>1024), mesh creation is not a good approach for processing, because only the mesh by itself, takes a huge amount of memmory and needs a high-end graphics card to manage. Although Blender can handle a 50M polygons mesh with a high-end graphics card, but it’s not really an optimized way for data processing.
For high nside values (>256) it is recommended to create the mesh part-by-part; that means dividing npix array into several sections and making a mesh for each part (you can merge all these parts in Blender). It may take up to an hour for generating a nside=2048 mesh but it can also become very heavy(about 12 GB while with vertices merged it would be ~2.5GB at most). So divide it and merge part-by-part. we can write a script to do all that.
I would like to say thanks to Andrea Zonca who helped me with finding boundaries of pixels.
Importing the mesh into Blender (Intro to Geometry Nodes)
This section is for complete begginners in Geometry Nodes. If you already know about that, skip this part.
Open Blender, and now you see default objects. Click on the 3dView and press A to select them all and press X and Enter to confirm deletion.
Now click on File > import > Wavefront(.obj) and head to the directory of your OBJ file you created and import it. The imege below shows a HEALPix mesh with nside = 8:
As explained before the file contains lots of redundant vertices since neighboring pixels have duplicate vertices. To this end we could use several different tools but for getting familiar to Blender’s Geometry Nodes(from now on geo-nodes), we use this tool to merge those redundant vertices.
At the top bar in the far right find Geometry Nodes tab and click it.
If you don’t find Geometry Nodes tab, hover your cursor on the top bar and hold MMB to pan the bar to the left to find Geometry Nodes. It is worth noting that anywhere in Blender(in 3dView, panels etc.), MMB is used for panning, Ctrl+MMB for zooming and Ctrl+Space for maximizing view.
Now select the HEALPix sphere you imported, and click on the “+ New” button on top of the geo-nodes panel. Now you will see an input mesh and an output that are the same at the moment:
Now we want to add a node to merge duplicate vertices. Press Shift+A to open add-pannel.
Shift+A everywhere opens add menu.
At the top click on search and search merge to find Merge By Distance and press Enter. Now move it to the middle of the process and click to place it like this:
On the far right you will find Properties panel. Click on the modifier’s dropdown button and select apply.
Now you have merged all redundant vertices. Save the file somewhere for later use.
Importing map data and setting mesh attributes
In this section we are going to extract data from CMB maps and import them as face data of the mesh. First you have to extract data from e.g. Planck maps. To do this, first download one of the provided maps e.g. Commander.
Now in system’s python, read Plack map and extract its temperature, stokes_q and stokes_u data:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Extracting data from .fits files
import numpy as np
import healpy as hp
nside = 64
fpath = "E:/Your/Path/" # change it for your own case
fpath += "COM_CMB_IQU-commander_2048_R3.00_full.fits"
_fields = (5, 6, 7) # inpainted fields
_names = ("temp", "stokes_q", "stokes_u")
for name, field in zip(_names, _fields):
map = hp.read_map(fpath, field = field, nest=False)
map = hp.ud_grade(map, nside_out=nside, order_in='RING')
map *= 10**6
np.savetxt("./" + name + ".txt", map) # format doesn't matter that much
Having data extracted, open Blender > stored mesh file, pan the top bar and select Scripting. In the Text Editor create a new script and enter the bellow code into it:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Assigning attributes to mesh
import numpy as np
import bpy
path = "E:/Your/Path/To/" # change it for your own case
obj_attr = bpy.context.object.data.attributes # active object's attributes
_names = ("temp", "stokes_q", "stokes_u")
_attributes = [np.loadtxt(path + name +".txt") for name in _names]
for name, attr in zip(_names, _attributes):
obj_attr.new(name=name, type='FLOAT', domain='FACE')
obj_attr[name].data.foreach_set('value', attr)
The above code assigns your each pixel’s data into its analogous face in the created mesh. Assigning attributes to the mesh lets you use the power of geo-nodes in parallel computations and also visualization(in realtime!).
Visualizing assigned attributes
In this section we want to achieve a visualization mechanism. The result is like the image bellow:
Visualization nodes
Get back to Geometry Nodes tab and create a new one. Create a geo-nodes graph like the one below(I will explain what all these do):
We proceed from left to right. First, the Named Attributes nodes takes the data values we assigned to mesh(in previous section) and makes it visible for geo-nodes. The Attributes Statistics node reads these values and gives the useful outputs you see (Note that it needs the Geometry wire to know which mesh we are looking for). Then we map minimum & maximum values of the temp attribute into the range (0, 1) by Map Range node. The ColorRamp node maps these values to a set of color data(attributes) that can be read in materials as we will see.
For creating this color ramp, click on the ‘+’ button to create a new tick and enter 0.5 for its position. Then change the color of the first tick to RGB = (0,0,1), the middle one to RGB = (0,1,0) and the last tick to RGB = (1,0,0). This gives a similar gradient to the jet coloring.
connect the output of the color ramp to the free output socket. Press the N key to open the right panel, go to Group tab and rename the output to what you want e.g. temp_out_ch. Note that the type of this output is color.
Now in the Properties > Modifier (the far right panel) put a name for the attribute that this output creates e.g. temp_col.
Material
Having the object selected, go to Properties > Material, Click + New button to create a new material. Change the surface material type to Emission.
Then, by clicking on the circle button near the color property, change the color type to Attribute. In the Name field, type temp_col(or whatever you named the attribute).
Now in the 3dView hold the Z key and in the pop-up menu select Material Preview and you have to see the CMB sphere.