Marker code of the 'scatter-3D' example

[QhelDIV]

Thanks for your book! The figures are amazing.
Since I am working in 3D and have to play around with point clouds a lot, I really want to know how the 'Scatter-3D' exmaple is made.
It seems the core part is the custom marker. The marker's fading shadow rings create a beautiful effect of 3D shading.
Can you share the code for creating this example? Thanks a lot!

[rougier]

Sorry for delay. This is actually a regular scatter plot but the trick is to render each point several times using the disc marker, from back to front. "Halo" is made of transparent black disc with larger radius and the highlight is made with a smaller shifted white disc. Quite slow and primite but it gives the illusion. I thought the code was available from https://github.com/rougier/matplotlib-3d but I did not find it. Where did you find the image ? I could give me some hints on where to search the code.

[QhelDIV]

Hi, thank you for your reply!

The image comes from the gallery show cases in this repo's README.md, the url is: https://github.com/rougier/scientific-visualization-book/blob/master/images/scatter-3d.png

Thank you again!

[rougier]

Cannot find the original but here is another similar on with a different technique:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection, LineCollection


def frustum(left, right, bottom, top, znear, zfar):
    M = np.zeros((4, 4))
    M[0, 0] = +2.0 * znear / (right - left)
    M[2, 0] = (right + left) / (right - left)
    M[1, 1] = +2.0 * znear / (top - bottom)
    M[2, 1] = (top + bottom) / (top - bottom)
    M[2, 2] = -(zfar + znear) / (zfar - znear)
    M[3, 2] = -2.0 * znear * zfar / (zfar - znear)
    M[2, 3] = -1.0
    return M.T

def perspective(fovy, aspect, znear, zfar):
    h = np.tan(fovy / 360.0 * np.pi) * znear
    w = h * aspect
    return frustum(-w, w, -h, h, znear, zfar)

def scale(x, y, z):
    return np.array([[x, 0, 0, 0],
                     [0, y, 0, 0],
                     [0, 0, z, 0],
                     [0, 0, 0, 1]], dtype=float)
def zoom(z):
    return scale(z,z,z)

def translate(x, y, z):
    return np.array([[1, 0, 0, x],
                     [0, 1, 0, y],
                     [0, 0, 1, z],
                     [0, 0, 0, 1]], dtype=float)

def xrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return np.array([[1, 0,  0, 0],
                     [0, c, -s, 0],
                     [0, s,  c, 0],
                     [0, 0,  0, 1]], dtype=float)

def yrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return  np.array([[ c, 0, s, 0],
                      [ 0, 1, 0, 0],
                      [-s, 0, c, 0],
                      [ 0, 0, 0, 1]], dtype=float)

def zrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return np.array([[ c, -s, 0, 0],
                     [ s,  c, 0, 0],
                     [ 0,  0, 1, 0],
                     [ 0,  0, 0, 1]], dtype=float)


n = 11
V = np.zeros((12,n,3))
V[0] = [(x,0,0) for x in np.linspace(0,1,n)] # (0,0,0) to (1,0,0)
V[1] = [(x,0,1) for x in np.linspace(0,1,n)] # (0,0,1) to (1,0,1)
V[2] = [(0,0,z) for z in np.linspace(0,1,n)] # (0,0,0) to (0,0,1)
V[3] = [(1,0,z) for z in np.linspace(0,1,n)] # (1,0,0) to (1,0,1)
V[4] = [(0,y,0) for y in np.linspace(0,1,n)] # (0,0,0) to (0,1,0)
V[5] = [(1,y,0) for y in np.linspace(0,1,n)] # (1,0,0) to (1,1,0)
V[6] = [(x,1,0) for x in np.linspace(0,1,n)] # (0,1,0) to (1,1,0)
V[7] = [(1,y,1) for y in np.linspace(0,1,n)] # (1,0,0) to (1,1,0)
V[8] = [(1,1,z) for z in np.linspace(0,1,n)] # (0,1,0) to (1,1,0)

V[9]  = [(x,0,1.025) for x in np.linspace(0,1,n)]
V[10] = [(-0.025,0,z) for z in np.linspace(0,1,n)]
V[11] = [(-0.025,y,0) for y in np.linspace(0,1,n)]

V = V - 0.5


model = xrotate(25) @ yrotate(45) 
view  = translate(0, 0,-3.5)
proj  = perspective(30, 1, 1, 100) 
MVP   = proj  @ view  @ model 

def transform(V, depth=False):
    V = np.asarray(V).reshape(-1,3)
    VH = np.c_[V, np.ones(len(V))]     # Homogenous coordinates
    VT = (VH @ MVP.T)                  # Transformed coordinates
    VN = VT/abs(VT[:,3].reshape(-1,1)) # Normalization
    if depth:
        return VN[:,:2], VN[:,2]
    else:
        return VN[:,:2]
    
V = transform(V).reshape((12,n,2))

# 3 axis + grid (at once)
segments, linewidths = [], []
for i in range(0,n):
    if 0 < i < (n-1): linewidth = 0.25
    else:             linewidth = 1.00
    for (j,k) in [(0,1), (2,3), (4,5), (6,0), (7,5), (8,3)]:
        triangle = [V[j,i], V[k,i]]
        segments.append(triangle)
        linewidths.append(linewidth)

# Ticks
for i in range(0,n):
    triangle = [V[1,i], V[9,i]]
    segments.append(triangle)
    linewidths.append(1.0)

    triangle = [V[2,i], V[10,i]]
    segments.append(triangle)
    linewidths.append(1.0)

    triangle = [V[4,i], V[11,i]]
    segments.append(triangle)
    linewidths.append(1.0)



    
fig = plt.figure(figsize=(6,6))
ax = fig.add_axes([0,0,1,1], xlim=[-1,1], ylim=[-1,1], aspect=1)
ax.axis("off")
collection = PolyCollection(segments, linewidths=linewidths, 
                            facecolors="None", edgecolor="black")
ax.add_collection(collection)

# Tick labels
P = [(x, -0.5,+0.55) for x in np.linspace(-0.5,0.5,n)]
for i,pos in enumerate(transform(P)):
    ax.text(pos[0], pos[1], "%d" % (i-10),
            ha="center", va="center", size="x-small")

P = [(-0.55, -0.5, z) for z in np.linspace(-0.5,0.5,n)]
for i,pos in enumerate(transform(P)):
    ax.text(pos[0], pos[1], "%d" % i,
            ha="center", va="center", size="x-small")

P = [(-0.55, y, -0.50) for y in np.linspace(-0.5,0.5,n)]
for i,pos in enumerate(transform(P)):
    ax.text(pos[0], pos[1], "%d" % i,
            ha="center", va="center", size="x-small")


normalize = lambda x: np.asarray(x)/np.linalg.norm(x)
n = 256
I = np.zeros((n,n,4))
X,Y = np.linspace(-1,1,n), np.linspace(-1,1,n)
color = np.array([1,0,0])
white = np.array([1,1,1])
direction = normalize([1,1,1])
for i,x in enumerate(X):
    for j,y in enumerate(Y):
        z = 1 - np.sqrt(x*x+y*y)
        if z < 0: continue
        normal = normalize([x,y,z])
        diffuse = max(0,(direction*normal).sum())
        specular = np.power(diffuse,24)
        I[j,i,:3] = np.maximum(diffuse*color, specular)
        I[j,i,3] = 1


from matplotlib.offsetbox import OffsetImage, AnnotationBbox


mol = np.load("protein.npy")
V = mol["position"]
Z = mol["radius"]
V = (2*(V-V.min())/(V.max()-V.min()) - 1 )*.75
V = np.random.normal(0, .2, (1500,3))


V,Z = transform(V, depth=True)
zmin,zmax = Z.min(), Z.max()
Z = (Z-zmin)/(zmax-zmin)
V = V[np.argsort(-Z)]
for i in range(len(V)):
    image = OffsetImage(I, zoom=0.05, #*Z[i],
                        origin="lower", interpolation="nearest")
    box = AnnotationBbox(image, V[i], frameon=False)
    ax.add_artist(box)

plt.show()

[rougier]

Found it: https://github.com/rougier/scientific-visualization-book/blob/master/code/unsorted/3d/scatter.py

[QhelDIV]

Thank you so much! I can directly render new point sets without any problem :)