知乎上的帖子，说MATLAB画的shit都比Python好看：
我想可以翻译成R来画一下，先照着写一个产生数据的函数：
shitshape <- function(points = 400, circle_pts = 40) {
linspace <- function(a, b, n) seq(a, b, length.out = n)
cross3 <- function(a, b) {
c(
a[2] * b[3] - a[3] * b[2],
a[3] * b[1] - a[1] * b[3],
a[1] * b[2] - a[2] * b[1]
)
}
normalize <- function(v) {
v / sqrt(sum(v * v))
}
z <- linspace(0, 1, points)
z[1:50] <- z[1:50] + z[50] * linspace(1, 0, 50)
z[201:400] <- z[201:400] - 0.12 * linspace(0, 1, 200)
s <- linspace(0, 1, points)
path_radius <- linspace(0.7, 0.02, points) * (cos(linspace(0, pi / 2, points)) + 1.5) / 3.5
tube_radius <- 0.18 * (1 - s^2.5) + 0.001
tube_radius[1:50] <- sin(linspace(0, pi / 2, 50)) * tube_radius[50]
theta <- linspace(0, 2 * pi * 4, points)
circle_theta <- linspace(0, 2 * pi, circle_pts)
X <- matrix(0, nrow = circle_pts, ncol = points)
Y <- matrix(0, nrow = circle_pts, ncol = points)
Z <- matrix(0, nrow = circle_pts, ncol = points)
for (i in 1:(points - 1)) {
R <- path_radius[i]
center_x <- R * cos(theta[i])
center_y <- R * sin(theta[i])
center_z <- z[i]
dx <- path_radius[i + 1] * cos(theta[i + 1]) - center_x
dy <- path_radius[i + 1] * sin(theta[i + 1]) - center_y
dz <- z[i + 1] - center_z
tangent <- normalize(c(dx, dy, dz))
ref <- c(0, 0, 1)
if (abs(sum(ref * tangent)) > 0.99) ref <- c(1, 0, 0)
normal1 <- normalize(cross3(tangent, ref))
normal2 <- cross3(tangent, normal1)
r_tube <- tube_radius[i]
for (j in 1:circle_pts) {
offset <- r_tube * (cos(circle_theta[j]) * normal1 + sin(circle_theta[j]) * normal2)
X[j, i] <- center_x + offset[1]
Y[j, i] <- center_y + offset[2]
Z[j, i] <- center_z + offset[3]
}
}
Z[, points] <- Z[, points - 1]
list(X = X, Y = Y, Z = Z)
}然后就是3D画图了：
library(rgl)
xyz <- shitshape()
open3d()
bg3d("white")
surface3d(xyz$X, xyz$Y, xyz$Z, color = rgb(0.4, 0.2, 0), lit = TRUE)
aspect3d(1, 1, 1)
light3d(theta = -60, phi = 10)
view3d(theta = -60, phi = 10)
确实是比Python的好看。
评论想看进阶版本，我翻一下，这不就是MATLAB的原创作者嘛，必须安排上。
把他/她的MATLAB 版本给翻译成R，结构与参数基本一一对应，MATLAB我还是用过的，虽然经历比较短。
library(rgl)
rot_x <- function(theta) {
matrix(
c(
1, 0, 0,
0, cos(theta), -sin(theta),
0, sin(theta), cos(theta)
),
nrow = 3,
byrow = TRUE
)
}
rot_y <- function(theta) {
matrix(
c(
cos(theta), 0, -sin(theta),
0, 1, 0,
sin(theta), 0, cos(theta)
),
nrow = 3,
byrow = TRUE
)
}
rot_z <- function(theta) {
matrix(
c(
cos(theta), -sin(theta), 0,
sin(theta), cos(theta), 0,
0, 0, 1
),
nrow = 3,
byrow = TRUE
)
}
rotate_xyz <- function(X, Y, Z, R) {
pts <- cbind(as.vector(X), as.vector(Y), as.vector(Z))
npts <- pts %*% R
list(
X = matrix(npts[, 1], nrow = nrow(X), ncol = ncol(X)),
Y = matrix(npts[, 2], nrow = nrow(Y), ncol = ncol(Y)),
Z = matrix(npts[, 3], nrow = nrow(Z), ncol = ncol(Z))
)
}
rotate_xyz_col <- function(X, Y, Z, R) {
pts <- cbind(as.vector(X), as.vector(Y), as.vector(Z))
npts <- t(R %*% t(pts))
list(
X = matrix(npts[, 1], nrow = nrow(X), ncol = ncol(X)),
Y = matrix(npts[, 2], nrow = nrow(Y), ncol = ncol(Y)),
Z = matrix(npts[, 3], nrow = nrow(Z), ncol = ncol(Z))
)
}
bezier_curve <- function(ctrl, n = 50) {
t <- seq(0, 1, length.out = n)
p <- nrow(ctrl) - 1
basis <- vapply(
0:p,
function(k) choose(p, k) * (t^k) * ((1 - t)^(p - k)),
numeric(length(t))
)
basis %*% ctrl
}
draw_straw <- function(X, Y, Z, color = rgb(88, 130, 126, maxColorValue = 255), lwd = 2) {
idx <- which(Z == min(Z, na.rm = TRUE), arr.ind = TRUE)[1, ]
x1 <- X[idx[1], idx[2]]
y1 <- Y[idx[1], idx[2]]
z1 <- Z[idx[1], idx[2]] + 0.03
ctrl <- rbind(
c(x1, y1, z1),
c(0, 0, -0.7),
c((x1 * cos(pi / 3) - y1 * sin(pi / 3)) / 3, (y1 * cos(pi / 3) + x1 * sin(pi / 3)) / 3, -1.5)
)
p <- bezier_curve(ctrl, n = 50)
lines3d(p[, 1], p[, 2], p[, 3], color = color, lwd = lwd)
}
set_color_by_h <- function(H, c_list, n = 256) {
hmin <- min(H, na.rm = TRUE)
hmax <- max(H, na.rm = TRUE)
base_cols <- grDevices::rgb(c_list[, 1], c_list[, 2], c_list[, 3])
pal <- grDevices::colorRampPalette(base_cols)(n)
if (!is.finite(hmin) || !is.finite(hmax) || hmax <= hmin) {
return(matrix(pal[n], nrow = nrow(H), ncol = ncol(H)))
}
x <- (H - hmin) / (hmax - hmin)
idx <- pmin(n, pmax(1, floor(x * (n - 1)) + 1))
matrix(pal[as.vector(idx)], nrow = nrow(H), ncol = ncol(H))
}
make_waffle_cone <- function(N = 500, n = 30) {
xs <- seq(-1, 1, length.out = N)
ys <- seq(-1, 1, length.out = N)
X <- matrix(rep(xs, each = N), nrow = N, ncol = N)
Y <- matrix(rep(ys, times = N), nrow = N, ncol = N)
Z <- X * 0
for (k in 1:2) {
Z[seq(k, N, by = n), ] <- 0.01
Z[, seq(k, N, by = n)] <- 0.01
}
Z[X^2 + Y^2 > 1] <- NA_real_
R_x <- rot_x(-pi / 2.9)
nXYZ <- cbind(as.vector(X), as.vector(Y), as.vector(Z)) %*% R_x
nX <- nXYZ[, 1]
nY <- nXYZ[, 2]
nZ <- nXYZ[, 3]
nZ <- nZ - min(nZ, na.rm = TRUE)
nY <- nY - min(nY, na.rm = TRUE)
denom <- pmax(nZ / 2.5, 1e-6)
nT <- nX / denom
nR <- nY
nX2 <- cos(nT) * nR
nY2 <- sin(nT) * nR
list(
X = matrix(nX2, nrow = N, ncol = N),
Y = matrix(nY2, nrow = N, ncol = N),
Z = matrix(nZ, nrow = N, ncol = N)
)
}
make_cylinder <- function(radius = 0.04, segments = 100, height = 1) {
theta <- seq(0, 2 * pi, length.out = segments + 1)
X <- outer(cos(theta), c(radius, radius))
Y <- outer(sin(theta), c(radius, radius))
Z <- outer(rep(1, length(theta)), c(0, height))
list(X = X, Y = Y, Z = Z)
}
poop2 <- function() {
open3d()
bg3d("white")
material3d(specular = "black")
poop <- shitshape()
surface3d(poop$X + 0.1, poop$Y, poop$Z + 1.06, color = rgb(0.4, 0.2, 0), lit = TRUE)
surface3d(poop$X / 1.3, poop$Y / 1.3, -poop$Z + 1.1, color = rgb(0.4, 0.2, 0), lit = TRUE)
cone <- make_waffle_cone(N = 500, n = 30)
surface3d(cone$X, cone$Y, cone$Z, color = rgb(228, 200, 142, maxColorValue = 255), lit = TRUE)
cyl <- make_cylinder(radius = 0.04, segments = 100, height = 1)
r1 <- rotate_xyz(cyl$X, cyl$Y, cyl$Z, rot_x(pi / 5))
surface3d(r1$X, r1$Y, r1$Z + 1.3, color = rgb(0.8, 0.6, 0.4), lit = TRUE)
r2 <- rotate_xyz(cyl$X, cyl$Y, cyl$Z, rot_x(pi / 6))
r2 <- rotate_xyz(r2$X, r2$Y, r2$Z, rot_y(pi / 7))
surface3d(r2$X, r2$Y, r2$Z + 1.2, color = rgb(0.8, 0.6, 0.4), lit = TRUE)
aspect3d(1, 1, 1)
light3d(theta = -60, phi = 10)
view3d(theta = -60, phi = 10)
}
poop3 <- function() {
open3d()
bg3d("white")
material3d(specular = "black")
a <- shitshape()
Xa <- a$X
Ya <- a$Y
Za <- a$Z
rb <- seq(0, 1, by = 0.01)
tb <- seq(0, 2, length.out = 151)
f <- (abs(1 - ((tb * 5) %% 2)) / 2) + 0.3
wb <- outer(rb, f, "*")
xb <- wb * matrix(cos(pi * tb), nrow = length(rb), ncol = length(tb), byrow = TRUE)
yb <- wb * matrix(sin(pi * tb), nrow = length(rb), ncol = length(tb), byrow = TRUE)
Zb <- (-cos(pi * wb * 1.2) + 1)^0.2
color_list <- rbind(
c(1.00, 0.92, 0.96),
c(1.00, 0.76, 0.90),
c(1.00, 0.60, 0.84),
c(0.98, 0.44, 0.76),
c(0.90, 0.26, 0.62),
c(0.78, 0.12, 0.48)
)
col_b <- set_color_by_h(Zb, color_list)
yaw_z <- 72 * pi / 180
roll_x_1 <- pi / 8
roll_x_2 <- pi / 9
R_z_2 <- rot_z(yaw_z)
R_z_1 <- rot_z(yaw_z / 2)
R_z_3 <- rot_z(yaw_z / 3)
R_x_1 <- rot_x(roll_x_1)
R_x_2 <- rot_x(roll_x_2)
surface3d(Xa, Ya, Za + 0.7, color = rgb(0.4, 0.2, 0), lit = TRUE)
r <- rotate_xyz_col(Xa, Ya, Za + 0.7, R_x_1)
r$Y <- r$Y - 0.4
surface3d(r$X, r$Y, r$Z - 0.1, color = rgb(0.4, 0.2, 0), lit = TRUE)
draw_straw(r$X, r$Y, r$Z - 0.1)
for (k in 1:4) {
r <- rotate_xyz_col(r$X, r$Y, r$Z, R_z_2)
surface3d(r$X, r$Y, r$Z - 0.1, color = rgb(0.4, 0.2, 0), lit = TRUE)
draw_straw(r$X, r$Y, r$Z - 0.1)
}
b <- rotate_xyz_col(xb / 2.5, yb / 2.5, Zb / 2.5 + 0.32, R_x_2)
b$Y <- b$Y - 1.35
for (k in 1:5) {
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_2)
surface3d(b$X, b$Y, b$Z, color = as.vector(col_b), lit = TRUE)
draw_straw(b$X, b$Y, b$Z)
}
b <- rotate_xyz_col(xb / 2.5, yb / 2.5, Zb / 2.5 + 0.32, R_x_2)
b$Y <- b$Y - 1.15
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_1)
for (k in 1:5) {
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_2)
surface3d(b$X, b$Y, b$Z, color = as.vector(col_b), lit = TRUE)
draw_straw(b$X, b$Y, b$Z)
}
b <- rotate_xyz_col(xb / 2.5, yb / 2.5, Zb / 2.5 + 0.32, R_x_2)
b$Y <- b$Y - 1.25
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_3)
for (k in 1:5) {
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_2)
surface3d(b$X, b$Y, b$Z, color = as.vector(col_b), lit = TRUE)
draw_straw(b$X, b$Y, b$Z)
}
b <- rotate_xyz_col(xb / 2.5, yb / 2.5, Zb / 2.5 + 0.32, R_x_2)
b$Y <- b$Y - 1.25
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_3)
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_3)
for (k in 1:5) {
b <- rotate_xyz_col(b$X, b$Y, b$Z, R_z_2)
surface3d(b$X, b$Y, b$Z, color = as.vector(col_b), lit = TRUE)
draw_straw(b$X, b$Y, b$Z)
}
aspect3d(1, 1, 1)
light3d(theta = -15, phi = 35)
light3d(theta = 120, phi = 20)
view3d(theta = -15, phi = 35)
}我把shit这个词改成用poop，毕竟shit是脏话，poop就不是了，是幼儿友好的，等同于中文的粑粑。
我们只要运行一个poop2()就可以出来巧克力味的冰激凌了，还带两饼干🍪：
然后运行poop3()，就能够得到花朵簇拥的一堆粑粑：
