Volume rendering comparison

ghost commented 2 years ago

I'm trying out volume rendering in different libraries, namely this one and yt (https://yt-project.org/doc/visualizing/volume_rendering.html)

When rendering the same dataset with the same transfer function, I get different results.

Dataset is created manually, and is supposed to represent a spherical volume:

import k3d
import yt
import numpy as np
import ipywidgets

nx, ny, nz = 110, 110, 110
volume = np.zeros((nx, ny, nz))
center = np.array([55, 55, 55])
for i in range(nx):
    for j in range(ny):
        for k in range(nz):
            if np.linalg.norm(np.array([i, j ,k]) - center) < 40:
                volume[i, j, k] = 1.0

Transfer function is initialized via the yt libary:

tf = yt.visualization.volume_rendering.transfer_functions.ColorTransferFunction(
    x_bounds=(np.min(volume), np.max(volume) + 0.1) # 0.1 offset needed to prevent rendering issues in yt
)
tf.map_to_colormap(
    tf.x_bounds[0],
    tf.x_bounds[1],
    colormap='Blues',
    scale_func=lambda vals, minval, maxval: (vals - vals.min()) / (vals.max() - vals.min()),
)

Which results in the following transfer function:

Where the opacity increases linearly from 0.0 to 1.0 and the end color is blue.

I convert the transfer function into a colormap and opacity function suitable for k3d

# create colormap from tf
colormap = np.zeros((256, 4), np.float32)
colormap[:, 0] = tf.alpha.x  /  max(tf.alpha.x) # rescale back between 0.0 and 1.0
colormap[:, 1] = tf.red.y
colormap[:, 2] = tf.green.y
colormap[:, 3] = tf.blue.y
# create opacity func
opacity_func = np.zeros((256, 2), dtype=np.float32)
opacity_func[:, 0] = tf.alpha.x /  max(tf.alpha.x) # rescale back between 0.0 and 1.0
opacity_func[:, 1] = np.linspace(0, 1, 256)
color_range = (np.min(volume), np.max(volume))

And feed to k3d.volume

out = ipywidgets.Output(layout={'width': '600px', 'height': '600px'})
plot = k3d.plot(background_color=0, grid_visible=False, lighting=0)
plot += k3d.volume(
    volume=volume,
    color_range=color_range,
    opacity_function=opacity_func,
    color_map=colormap,
    alpha_coef=1.0
)
with out:
    display(plot)
display(out)
plot.camera_reset(factor=0.5)

I also feed the same data into yt (setting their camera lens to 'perspective' to match k3d)

data = dict(density = (volume, "g"))
ds = yt.load_uniform_grid(data, volume.shape, length_unit="m", nprocs=1)
sc = yt.create_scene(ds, 'density')
dd = ds.all_data()
sc[0].log_field = False
sc[0].set_transfer_function(tf)
cam = sc.add_camera(ds, lens_type="perspective")
cam.resolution = [600, 600]
sc.show(sigma_clip=0)

The images look very different, and I have some questions:

Why does k3d not saturate the rays that pass through the middle of the sphere? I would expect the rays to be oversaturated with blue there, and mixing many contributions of the darker blue eventually leads to light blue?
Why is there almost no gradient in color between the edges of the sphere (where the rays pass only through a small section of the volume) and the middle
Why does the volume seem to consist of little cubes? Even if I set the lighting to 0, the volume seems to be made of little cubes that have a 'dark' and a 'bright' side

I found out that if I use the sigma_clip option in yt, I can make the renders look more similar:

sigma_clip = N can address this by removing values that are more than N standard deviations brighter than the mean of your image

Is this something that is built-in to k3d as well?

I hope this question is not to long and complicated. I included a notebook that can render the examples as well. Any help would be much appreciated! k3d_vs_yt.ipynb.gz

ghost commented 2 years ago

The maintainers over at yt suspect the difference is k3d uses a 'max intensity transfer function', meaning that it will take the highest value along the ray, while yt integrates over the ray path, leading to saturation.

Is this indeed the case?

artur-trzesiok commented 2 years ago

Hi!

Thanks for analysis. @MironLeon k3d use both method but for maximum intensity projection we have k3d.mip like here: https://github.com/K3D-tools/K3D-jupyter/blob/main/examples/volume_renderer.ipynb

embryo = k3d.mip(volume_data.astype(np.float16), 
                 color_map=np.array(k3d.basic_color_maps.BlackBodyRadiation, dtype=np.float32), 
                 bounds=bounds)

plot = k3d.plot(background_color=0, grid_visible=False)
plot += embryo
plot.display()

Why does k3d not saturate the rays that pass through the middle of the sphere? I would expect the rays to be oversaturated with blue there, and mixing many contributions of the darker blue eventually leads to light blue?

About color question. Volume rendering is about transparency in dense data. let's simplify our case to 5 dark blue translucent glass plates stacked on top of each other. When you put them on, will you suddenly see a light blue color? In my opinion you will not see.

YT must use the summation of the color values multiplied by the correction factor for the redering, so that the color appears light blue (additive model).

Please run simple code to export data:

import SimpleITK as sitk

writer = sitk.ImageFileWriter()
writer.SetFileName("test.mhd")
writer.Execute(sitk.GetImageFromArray(volume))

right now We can load mhd file intro paraview (https://www.paraview.org/):

I don't see bright blue here. So from my perspective its look like we have good volume rendering but image that you expect to see is another type of algo - We definetly can think about add it into k3d :)

Why is there almost no gradient in color between the edges of the sphere (where the rays pass only through a small section of the volume) and the middle

Please check alpha_coef parameter to k3d.volume

Why does the volume seem to consist of little cubes? Even if I set the lighting to 0, the volume seems to be made of little cubes that have a 'dark' and a 'bright' side

It is come from resolution 110^3 . There is also chance that k3d cannot run data interpolation on your computer. Pleaes notice that Paraview have "even more cubes"

Some graphics card have limitation of features if texture resolution are not power of 2. Please check it with 128^3 volume data :)

Thanks for very good question - I like it! hope I put a little bit of light here. Please don't hesitate to ask any questions!

Best regards,

artur-trzesiok commented 2 years ago

Btw outcome from k3d.mip is:

ghost commented 2 years ago

Thanks for the detailed explanation.

I think the cubes come from the type of data I created, where all values are either 0 or 1. If I change the data cube to be more smooth:

import k3d
import yt
import numpy as np
import ipywidgets

nx, ny, nz = 128, 128, 128
volume = np.zeros((nx, ny, nz))
center = np.array([64, 64, 64])
for i in range(nx):
    for j in range(ny):
        for k in range(nz):
            d = np.linalg.norm(np.array([i, j ,k]) - center)
            if d < 60:
                volume[i, j, k] = 1.0 - d / 60

Both k3d and paraview produces smooth renders, where the cubes are not visible.

k3d use both method but for maximum intensity projection we have k3d.mip like here

So you mean to say the regular k3d.volume object does not use maximum intensity projection, only k3d.mip?

The main reason I am looking at these libraries is that I want to preview volume data in an interactive view, and have the possibility to render a timeseries into a movie. So far I have found only paraview can do both, but it is not web based. Do you maybe have a suggestion? Or is there a way for k3d to do headless rendering?

artur-trzesiok commented 2 years ago

Hi!

Yes, K3D have support to headless rendering - https://github.com/K3D-tools/K3D-jupyter/blob/main/examples/headless.ipynb . This is one of most recent features :)

Best regards

ghost commented 2 years ago

For the headless rendering, I am getting this error quite often:

~/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/k3d/headless.py in sync(self, hold_until_refreshed)
    108 
    109     def sync(self, hold_until_refreshed=False):
--> 110         self.browser.execute_script("k3dRefresh()")
    111 
    112         if hold_until_refreshed:

~/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/selenium/webdriver/remote/webdriver.py in execute_script(self, script, *args)
    870         command = Command.W3C_EXECUTE_SCRIPT
    871 
--> 872         return self.execute(command, {
    873             'script': script,
    874             'args': converted_args})['value']

~/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/selenium/webdriver/remote/webdriver.py in execute(self, driver_command, params)
    416         response = self.command_executor.execute(driver_command, params)
    417         if response:
--> 418             self.error_handler.check_response(response)
    419             response['value'] = self._unwrap_value(
    420                 response.get('value', None))

~/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/selenium/webdriver/remote/errorhandler.py in check_response(self, response)
    241                 alert_text = value['alert'].get('text')
    242             raise exception_class(message, screen, stacktrace, alert_text)  # type: ignore[call-arg]  # mypy is not smart enough here
--> 243         raise exception_class(message, screen, stacktrace)
    244 
    245     def _value_or_default(self, obj: Mapping[_KT, _VT], key: _KT, default: _VT) -> _VT:

JavascriptException: Message: javascript error: k3dRefresh is not defined
  (Session info: headless chrome=95.0.4638.69)
Stacktrace:
#0 0x563a25f7e463 <unknown>
#1 0x563a25a56678 <unknown>
#2 0x563a25a594fc <unknown>
#3 0x563a25a592f6 <unknown>
#4 0x563a25a59eb2 <unknown>
#5 0x563a25abd6a3 <unknown>
#6 0x563a25aa9792 <unknown>
#7 0x563a25abc9b1 <unknown>
#8 0x563a25aa9683 <unknown>
#9 0x563a25a7fb64 <unknown>
#10 0x563a25a80b55 <unknown>
#11 0x563a25fae0fe <unknown>
#12 0x563a25fc39d0 <unknown>
#13 0x563a25faf055 <unknown>
#14 0x563a25fc4e18 <unknown>
#15 0x563a25fa382b <unknown>
#16 0x563a25fdff98 <unknown>
#17 0x563a25fe0118 <unknown>
#18 0x563a25ffb59d <unknown>
#19 0x7f0188732609 <unknown>

Any idea what this relates to?

ghost commented 2 years ago

ERROR:k3d.headless:Exception on / [POST]
Traceback (most recent call last):
  File "/home/miron/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/flask/app.py", line 2073, in wsgi_app
    response = self.full_dispatch_request()
  File "/home/miron/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/flask/app.py", line 1518, in full_dispatch_request
    rv = self.handle_user_exception(e)
  File "/home/miron/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/flask/app.py", line 1516, in full_dispatch_request
    rv = self.dispatch_request()
  File "/home/miron/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/flask/app.py", line 1502, in dispatch_request
    return self.ensure_sync(self.view_functions[rule.endpoint])(**req.view_args)
  File "/home/miron/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/k3d/headless.py", line 101, in generate
    return Response(msgpack.packb(diff, use_bin_type=True),
  File "/home/miron/.pyenv/versions/3.8.7/envs/susipop/lib/python3.8/site-packages/msgpack/__init__.py", line 35, in packb
    return Packer(**kwargs).pack(o)
  File "msgpack/_packer.pyx", line 292, in msgpack._cmsgpack.Packer.pack
  File "msgpack/_packer.pyx", line 298, in msgpack._cmsgpack.Packer.pack
  File "msgpack/_packer.pyx", line 295, in msgpack._cmsgpack.Packer.pack
  File "msgpack/_packer.pyx", line 231, in msgpack._cmsgpack.Packer._pack
  File "msgpack/_packer.pyx", line 231, in msgpack._cmsgpack.Packer._pack
  File "msgpack/_packer.pyx", line 231, in msgpack._cmsgpack.Packer._pack
  File "msgpack/_packer.pyx", line 264, in msgpack._cmsgpack.Packer._pack
  File "msgpack/_packer.pyx", line 289, in msgpack._cmsgpack.Packer._pack
TypeError: can not serialize 'numpy.float32' object

This one also shows up quite often, but not consistently

tgandor commented 2 years ago

Interesting. I'm completely new to headless.py, but it looks like, when diff contains numpy arrays, it can fail. In a different project we used:

https://pypi.org/project/msgpack-numpy/

... for "transporting" ndarrays over msgpack.

(I'm not a great fan of msgpack; it's fast, but it's not like e.g. pickle -- it doesn't distinguish str and bytes, nor list from tuple).

Depending on the situation, this may be easy to handle by just adding this dependency, or harder because of corner cases (like mix of numpy with other data types).

tgandor commented 2 years ago

About missing k3dRefresh() problem I have no idea ;) @artur-trzesiok ? BTW, @MironLeon, if you could provide min working examples where these things happen, it may help us track this down. Anyway, thanks for reporting and the stack traces!

artur-trzesiok commented 2 years ago

Hi!

about: TypeError: can not serialize 'numpy.float32' object

msgpack have problems with serialization of numpy-float like values e.g: np.sum(np.arange(5)) is not a float/int but:

So if you pass to k3d any "number" not as python plain int/float it will crash :(. you have to do: int(np.sum(np.arange(5)))

It is annoying. @tgandor maybe on traitlets level we can handle that on k3d library level on setter?

About k3dRefresh() - that mean that headless js are not ready for some reasons. If K3D loading failed (for example because "TypeError: can not serialize 'numpy.float32' object") k3dRefresh will be not exposed.

I think that you spot two different glitch that comes from the same reason.

basilevh commented 2 years ago

@artur-trzesiok Could you please elaborate on the k3dRefresh() error? My plot works perfectly when I display and interact with it normally within the notebook, but as soon as I call headless.sync(), I get:

JavascriptException: Message: javascript error: k3dRefresh is not defined
  (Session info: headless chrome=100.0.4896.75)

artur-trzesiok commented 2 years ago

@basilevh do you have still that issue or new version resolved that?

basilevh commented 2 years ago

@artur-trzesiok I updated to k3d 2.14.1 but still have the same error:

JavascriptException: Message: javascript error: k3dRefresh is not defined
  (Session info: headless chrome=101.0.4951.41)

(PS the same happens when I run https://github.com/K3D-tools/K3D-jupyter/blob/main/examples/headless.ipynb verbatim, so let me open a new issue instead)

K3D-tools / K3D-jupyter

Volume rendering comparison #318