I have made a python script using this amazing addon Node2Python from the geometry nodes setup I made when I asked about this question. Shoutout to lemon thank you so much for your help! The geometry nodes to voxelize worked perfectly. At first when I generated the python script equivalent of the node setup and then executed it, you will see that the setup is exactly the same with the previous setup I have (The blend file is also available in the previous thread)
But when I apply this python generated setup, the colors output don't match anymore. I compared both the original GN version and the python generated GN version and I don't see any difference yet the output is not right. What's going on? I'm not sure if I reproduced a bug, can anyone confirm?
UPDATE: It looks like a bug as confirmed by lemon here
When you apply the modifier notice that there is an empty material added. This new empty material is added at the end in the original setup, hence not interfering with the setup. The python-generated node setup however inserts this empty material at index zero and thus disrupting the node setup by shifting all the material indices. So this is definitely a bug because of the inconsistent behavior despite identical GN setups.
Now the question is, does anyone know of a workaround to this?