如何用R生成对照3D曲面并随机布点以评估点集空间分布

Hey Andrea, let's tackle this problem step by step—you're on the right track with extracting coordinate ranges, so let's build out the full workflow for evaluating your 3D point distribution and generating control random points that match your original surface's area:

一、完成三维坐标范围提取 #

First, let's finish that coordinate range extraction you started. Assuming you're using R (given the snippet you shared), here's a robust version that handles missing values:

# Extract min/max for x, y, z coordinates
x_range <- range(as.numeric(df$x), na.rm = TRUE)
y_range <- range(as.numeric(df$y), na.rm = TRUE)
z_range <- range(as.numeric(df$z), na.rm = TRUE)

Adding na.rm = TRUE ensures missing values don't mess up your range calculations.

二、生成面积匹配的空曲面与随机布点 #

Since your original surface is an approximate rectangular 3D surface from image segmentation, we have two solid approaches to generate control points that match its area:

Approach 1: Fit an Approximate Plane (for near-flat surfaces) #

If your surface is close to planar, first fit a plane to your original points, then generate random x/y values within your extracted ranges and calculate corresponding z values to match the plane:

# Fit a linear plane model (z as a function of x + y)
plane_model <- lm(z ~ x + y, data = df)

# Generate the same number of random points as your original dataset
n_points <- nrow(df)
rand_x <- runif(n_points, min = x_range[1], max = x_range[2])
rand_y <- runif(n_points, min = y_range[1], max = y_range[2])

# Predict z values using the fitted plane to match the surface's geometry
rand_z <- predict(plane_model, newdata = data.frame(x = rand_x, y = rand_y))

# Combine into your control dataset
rand_df <- data.frame(x = rand_x, y = rand_y, z = rand_z)

This method ensures the control surface's area (both projected and approximate curved area) matches your original surface closely.

Approach 2: Density-Constrained Random Points (for curved surfaces) #

If your surface has noticeable curvature, use kernel density estimation to restrict random points to the same spatial region as your original surface:

library(MASS)

# Estimate the 2D kernel density of your original x/y coordinates
kd <- kde2d(df$x, df$y, n = 100)

# Generate random points, filtering to keep only those in high-density regions (matching your surface's bounds)
rand_points <- data.frame()
target_n <- nrow(df)

while(nrow(rand_points) < target_n) {
  # Generate batch of candidate points
  temp_x <- runif(100, x_range[1], x_range[2])
  temp_y <- runif(100, y_range[1], y_range[2])
  temp_z <- runif(100, z_range[1], z_range[2])
  
  # Keep points where density is above a threshold (adjust quantile to tweak tightness)
  dens <- interp.surface(kd, cbind(temp_x, temp_y))
  valid_points <- data.frame(x = temp_x, y = temp_y, z = temp_z)[dens > quantile(dens, 0.2), ]
  
  rand_points <- rbind(rand_points, valid_points)
}

# Trim to exactly the number of points in your original dataset
rand_df <- rand_points[1:target_n, ]

This method preserves the original surface's spatial footprint, ensuring the control surface's area matches.

三、评估空间分布（聚类检测与量化） #

With your original and control datasets ready, here are three ways to assess clustering and spatial patterns:

1. Density Visualization & Comparison #

Use kernel density plots to visually compare point clustering between your original and control data:

library(ggplot2)
library(gridExtra)

# Original data density plot
original_kd <- kde2d(df$x, df$y, n = 100)
p1 <- ggplot(df, aes(x, y)) +
  geom_tile(data = as.data.frame(original_kd), aes(x = x, y = y, fill = z)) +
  geom_point(size = 0.5, color = "white") +
  ggtitle("Original Point Density") +
  scale_fill_viridis_c()

# Control data density plot
rand_kd <- kde2d(rand_df$x, rand_df$y, n = 100)
p2 <- ggplot(rand_df, aes(x, y)) +
  geom_tile(data = as.data.frame(rand_kd), aes(x = x, y = y, fill = z)) +
  geom_point(size = 0.5, color = "white") +
  ggtitle("Random Control Point Density") +
  scale_fill_viridis_c()

# Display side-by-side
grid.arrange(p1, p2, ncol = 2)

If your original data has distinct high-density patches that aren't present in the control, that's clear evidence of clustering.

2. DBSCAN Clustering Detection #

Use the DBSCAN algorithm to explicitly detect clusters in your 3D point cloud:

library(dbscan)
library(plotly)

# Calculate 3D distance matrix
dist_matrix <- dist(df[, c("x", "y", "z")])

# Run DBSCAN (adjust eps and minPts based on your point scale; test different values!)
db_results <- dbscan(dist_matrix, eps = 0.5, minPts = 5)

# Print cluster counts (0 = noise points)
cat("Cluster counts:\n")
print(table(db_results$cluster))

# 3D visualization of clusters
plot_ly(df, x = ~x, y = ~y, z = ~z, 
        color = ~as.factor(db_results$cluster),
        type = "scatter3d", mode = "markers",
        marker = list(size = 3)) %>%
  layout(title = "3D Clustering of Original Points")

Compare this to the control dataset—you should see no meaningful clusters (only noise) in the random points if your original data is clustered.

3. Spatial Autocorrelation (Moran's I) #

Quantify spatial clustering using Moran's I, which measures how much points are correlated with their neighbors:

library(spdep)

# Calculate Moran's I for original data
original_coords <- df[, c("x", "y", "z")]
original_knn <- knn2nb(knearneigh(original_coords, k = 5))
original_weights <- nb2listw(original_knn, style = "W")
original_moran <- moran.test(df$x, original_weights)

# Calculate Moran's I for control data
rand_coords <- rand_df[, c("x", "y", "z")]
rand_knn <- knn2nb(knearneigh(rand_coords, k = 5))
rand_weights <- nb2listw(rand_knn, style = "W")
rand_moran <- moran.test(rand_df$x, rand_weights)

# Compare results
cat("Original Moran's I:", round(original_moran$estimate[1], 3), "\n")
cat("Control Moran's I:", round(rand_moran$estimate[1], 3), "\n")

A significantly higher Moran's I value for your original data indicates positive spatial autocorrelation (i.e., clustering).

内容的提问来源于stack exchange，提问作者Andrea Paterlini