Skip to main content Skills Marketplace Descubra e explore skills de IA criadas pela comunidade.
Instalar com Codex ou Claude Copie este prompt, cole no Codex, Claude ou outro assistente e deixe que ele revise a página da skill e instale para você.
Mostrar detalhes do prompt Um comando direto ignora o prompt de revisão. Verifique a origem antes de executá-lo.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill photogrammetry-3dO comando permanece em uma só linha. Role horizontalmente para revisá-lo antes de copiar.
Prefere uma cópia local? Baixe os arquivos disponíveis atualmente no SkillsMP.
Baixar Zip Baixando... Ocupações relacionadas SOC
Baseado na classificação ocupacional SOC
name photogrammetry-3d description Use this Skill for archaeological photogrammetry: 3D point cloud processing, mesh analysis, cross-section profiles, and volume measurement with Open3D.
tags ["archaeology","photogrammetry","3d-modeling","point-cloud","open3d"] version 1.0.0 authors [{"name":"Rosetta Skills Contributors","github":"@xjtulyc"}] license MIT platforms ["claude-code","codex","gemini-cli","cursor"] dependencies {"python":["open3d>=0.18","numpy>=1.24","scipy>=1.11","matplotlib>=3.7","pandas>=2.0"]} last_updated 2026-03-17 status stable
Archaeological Photogrammetry and 3D Analysis
One-line summary : Process archaeological 3D models from photogrammetry: point cloud cleaning, surface normal estimation, cross-section extraction, volume computation, and color analysis with Open3D.
When to Use This Skill
When processing photogrammetric point clouds of archaeological sites or objects
When computing surface areas and volumes from 3D meshes
When extracting cross-sectional profiles from artifacts or excavations
When performing outlier removal and surface smoothing
When analyzing color distribution on ceramic or painted surfaces
When computing roughness metrics for lithic use-wear analysis
Trigger keywords : photogrammetry, 3D model, point cloud, Open3D, mesh, SfM, structure from motion, cross-section, volume measurement, surface normal, artifact 3D, site survey, DEM, digital elevation model
Background & Key Concepts
Structure from Motion (SfM) Photogrammetry
SfM reconstructs 3D geometry from overlapping 2D photographs:
Feature detection and matching (SIFT/ORB)
Camera pose estimation
Bundle adjustment
Dense point cloud generation
Mesh reconstruction (Poisson or Delaunay)
Common software: Agisoft Metashape, COLMAP (free), RealityCapture.
Point Cloud Processing Pipeline
Load PLY/LAS/XYZ file
Statistical outlier removal (SOR filter)
Normal estimation (PCA over k-NN)
Surface reconstruction (Ball-Pivoting or Poisson)
Mesh simplification and hole filling
Roughness (Use-Wear Proxy)
$$
R_a = \frac{1}{N}\sum_{i=1}^N |z_i - \bar{z}|
$$
Higher roughness on lithic edges indicates use-wear; smooth facets = unused or polished.
Environment Setup
Install Dependencies
pip install open3d>=0.18 numpy>=1.24 scipy>=1.11 matplotlib>=3.7 pandas>=2.0
pip install laspy>=2.5
Verify Installation
import open3d as o3d
import numpy as np
print (f"Open3D version: {o3d.__version__} " )
1000
3
print
f"Test point cloud: {len (pcd.points)} points"
Core Workflow
Step 1: Point Cloud Loading, Cleaning, and Normal Estimation import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt
def create_ceramic_vessel_pcd (n_points=5000 , noise_level=0.002 ):
"""
Simulate a pottery vessel point cloud (truncated ellipsoid + rim).
"""
phi = np.random.uniform(0 , np.pi * 0.75 , n_points)
theta = np.random.uniform(0 , 2 *np.pi, n_points)
r_belly = 0.12
r_neck = 0.06
h_total = 0.25
t = phi / (np.pi * 0.75 )
r_profile = r_belly * np.sin(np.pi * t) + r_neck * (1 - np.sin(np.pi * t))
x = r_profile * np.sin(phi) * np.cos(theta)
y = r_profile * np.sin(phi) * np.sin(theta)
z = h_total * (1 - phi / (np.pi * 0.75 ))
xyz = np.column_stack([x, y, z])
xyz += np.random.randn(*xyz.shape) * noise_level
base_color = np.array([0.78 , 0.47 , 0.33 ])
colors = base_color + np.random.randn(n_points, 3 ) * 0.05
colors = np.clip(colors, 0 , 1 )
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz)
pcd.colors = o3d.utility.Vector3dVector(colors)
return pcd
pcd_raw = create_ceramic_vessel_pcd(n_points=5000 , noise_level=0.003 )
print (f"Raw point cloud: {len (pcd_raw.points)} points" )
pcd_cleaned, ind = pcd_raw.remove_statistical_outlier(nb_neighbors=20 , std_ratio=2.0 )
outlier_count = len (pcd_raw.points) - len (pcd_cleaned.points)
print (f"After SOR filtering: {len (pcd_cleaned.points)} points ({outlier_count} outliers removed)" )
pcd_cleaned.estimate_normals(
search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.02 , max_nn=30 )
)
pcd_cleaned.orient_normals_towards_camera_location(camera_location=[0 , 0 , 1 ])
print ("Surface normals estimated and oriented" )
pcd_down = pcd_cleaned.voxel_down_sample(voxel_size=0.005 )
print (f"After voxel downsampling: {len (pcd_down.points)} points" )
pts = np.asarray(pcd_down.points)
colors = np.asarray(pcd_down.colors) if pcd_down.has_colors() else np.ones((len (pts), 3 )) * 0.7
fig, axes = plt.subplots(1 , 3 , figsize=(14 , 5 ))
for ax, (i, j, title) in zip (axes, [(0 ,2 ,'Top view (XZ)' ), (0 ,1 ,'Side view (XY)' ), (1 ,2 ,'Front view (YZ)' )]):
ax.scatter(pts[:,i], pts[:,j], c=colors, s=2 , alpha=0.6 )
ax.set_xlabel(['X' ,'X' ,'Y' ][axes.tolist().index(ax)] + ' (m)' )
ax.set_ylabel(['Z' ,'Y' ,'Z' ][axes.tolist().index(ax)] + ' (m)' )
ax.set_title(title); ax.set_aspect('equal' ); ax.grid(True , alpha=0.3 )
plt.suptitle(f"Ceramic Vessel Point Cloud ({len (pts)} points after cleaning)" )
plt.tight_layout()
plt.savefig("point_cloud_views.png" , dpi=150 )
plt.show()
Step 2: Surface Reconstruction and Volume Computation import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt
pcd_down = create_ceramic_vessel_pcd(n_points=3000 , noise_level=0.001 )
pcd_down, _ = pcd_down.remove_statistical_outlier(nb_neighbors=20 , std_ratio=2.0 )
pcd_down.estimate_normals(
search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.02 , max_nn=30 )
)
pcd_down.orient_normals_towards_camera_location([0 , 0 , 1 ])
print ("Running Poisson surface reconstruction..." )
mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(
pcd_down, depth=8 , width=0 , scale=1.1 , linear_fit=False
)
densities_arr = np.asarray(densities)
density_threshold = np.percentile(densities_arr, 10 )
vertices_to_remove = densities_arr < density_threshold
mesh.remove_vertices_by_mask(vertices_to_remove)
mesh.remove_degenerate_triangles()
mesh.remove_non_manifold_edges()
print (f"Mesh: {len (mesh.vertices)} vertices, {len (mesh.triangles)} triangles" )
surface_area_m2 = mesh.get_surface_area()
print (f"Surface area: {surface_area_m2*1e4 :.2 f} cm²" )
if mesh.is_watertight():
volume_m3 = mesh.get_volume()
volume_ml = volume_m3 * 1e6
print (f"Volume: {volume_ml:.0 f} mL ({volume_ml/1000 :.3 f} L)" )
else :
print ("Mesh not watertight — attempting to fill holes..." )
mesh.fill_holes()
if mesh.is_watertight():
volume_m3 = mesh.get_volume()
print (f"Volume (after fill): {volume_m3*1e6 :.0 f} mL" )
else :
print ("Volume: mesh still not watertight; use alpha shape or convex hull approximation" )
hull, _ = pcd_down.compute_convex_hull()
vol_approx = hull.get_volume()
print (f"Convex hull volume (upper bound): {vol_approx*1e6 :.0 f} mL" )
pts = np.asarray(pcd_down.points)
z_min, z_max = pts[:,2 ].min (), pts[:,2 ].max ()
n_slices = 30
z_levels = np.linspace(z_min, z_max, n_slices)
radii = []
for z_lev in z_levels:
mask = np.abs (pts[:,2 ] - z_lev) < 0.005
if mask.sum () > 5 :
slice_pts = pts[mask, :2 ]
radii.append(np.linalg.norm(slice_pts, axis=1 ).mean())
else :
radii.append(np.nan)
radii = np.array(radii)
fig, axes = plt.subplots(1 , 2 , figsize=(12 , 5 ))
axes[0 ].plot(radii * 100 , z_levels * 100 , 'b-' , linewidth=2 )
axes[0 ].plot(-radii * 100 , z_levels * 100 , 'b-' , linewidth=2 )
axes[0 ].fill_betweenx(z_levels*100 , -radii*100 , radii*100 , alpha=0.2 , color='brown' )
axes[0 ].set_xlabel("Radius (cm)" ); axes[0 ].set_ylabel("Height (cm)" )
axes[0 ].set_title("Vessel Cross-Section Profile" ); axes[0 ].grid(True , alpha=0.3 )
axes[0 ].set_aspect('equal' )
axes[1 ].plot(z_levels * 100 , radii * 100 , 'ro-' , linewidth=1.5 , markersize=5 )
axes[1 ].set_xlabel("Height (cm)" ); axes[1 ].set_ylabel("Radius (cm)" )
axes[1 ].set_title("Radius Profile — Vessel Morphology" )
axes[1 ].grid(True , alpha=0.3 )
plt.tight_layout()
plt.savefig("vessel_profile.png" , dpi=150 )
plt.show()
print (f"\nVessel metrics:" )
print (f" Height: {(z_max-z_min)*100 :.1 f} cm" )
print (f" Maximum diameter: {np.nanmax(radii)*200 :.1 f} cm" )
print (f" Rim diameter: {radii[-1 ]*200 :.1 f} cm (at top)" )
print (f" Base diameter: {radii[0 ]*200 :.1 f} cm (at bottom)" )
Step 3: Color Analysis and Surface Roughness import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d
pcd = create_ceramic_vessel_pcd(n_points=5000 )
pcd, _ = pcd.remove_statistical_outlier(20 , 2.0 )
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.02 , max_nn=30 ))
pts = np.asarray(pcd.points)
colors = np.asarray(pcd.colors)
normals = np.asarray(pcd.normals)
z_min, z_max = pts[:,2 ].min (), pts[:,2 ].max ()
lower_third = pts[:,2 ] < z_min + (z_max-z_min)/3
upper_third = pts[:,2 ] > z_min + 2 *(z_max-z_min)/3
middle = ~lower_third & ~upper_third
zones = {'Base zone' : lower_third, 'Body zone' : middle, 'Neck/rim' : upper_third}
fig, axes = plt.subplots(2 , 2 , figsize=(13 , 9 ))
hue = colors[:,0 ]
sc = axes[0 ][0 ].scatter(pts[:,0 ]*100 , pts[:,1 ]*100 , c=colors, s=4 , alpha=0.7 )
axes[0 ][0 ].set_xlabel("X (cm)" ); axes[0 ][0 ].set_ylabel("Y (cm)" )
axes[0 ][0 ].set_title("Color Distribution (Top View)" ); axes[0 ][0 ].set_aspect('equal' )
for ax, (zone_name, mask) in zip ([axes[0 ][1 ], axes[1 ][0 ], axes[1 ][1 ]], zones.items()):
if mask.sum () > 0 :
zone_colors = colors[mask]
ax.hist(zone_colors[:,0 ], bins=30 , alpha=0.5 , color='red' , density=True , label='R' )
ax.hist(zone_colors[:,1 ], bins=30 , alpha=0.5 , color='green' , density=True , label='G' )
ax.hist(zone_colors[:,2 ], bins=30 , alpha=0.5 , color='blue' , density=True , label='B' )
ax.set_title(f"RGB Distribution — {zone_name} \n(n={mask.sum ()} points)" )
ax.set_xlabel("Intensity" ); ax.set_ylabel("Density" )
ax.legend(fontsize=8 ); ax.grid(True , alpha=0.3 )
plt.tight_layout()
plt.savefig("color_analysis.png" , dpi=150 )
plt.show()
normal_z = normals[:,2 ]
roughness_proxy = 1 - np.abs (normal_z)
print ("\nSurface roughness (normal deviation proxy):" )
print (f" Mean roughness: {roughness_proxy.mean():.4 f} " )
print (f" Std roughness: {roughness_proxy.std():.4 f} " )
print (f" Max roughness: {roughness_proxy.max ():.4 f} " )
print ("\nRoughness by zone:" )
for zone_name, mask in zones.items():
if mask.sum () > 0 :
r = roughness_proxy[mask].mean()
print (f" {zone_name} : mean roughness = {r:.4 f} " )
Advanced Usage
ICP Registration (Aligning Fragments) import open3d as o3d
import numpy as np
pcd1 = create_ceramic_vessel_pcd(500 )
pts2 = np.asarray(pcd1.points)[:250 ] + np.array([0.05 , 0.02 , 0.01 ])
pcd2 = o3d.geometry.PointCloud()
pcd2.points = o3d.utility.Vector3dVector(pts2)
pcd2.paint_uniform_color([0.8 , 0.3 , 0.3 ])
pcd1.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(0.02 , 30 ))
pcd2.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(0.02 , 30 ))
threshold = 0.02
trans_init = np.eye(4 )
reg = o3d.pipelines.registration.registration_icp(
pcd2, pcd1, threshold, trans_init,
o3d.pipelines.registration.TransformationEstimationPointToPoint(),
)
print (f"ICP registration fitness: {reg.fitness:.4 f} " )
print (f"ICP RMSE: {reg.inlier_rmse*1000 :.2 f} mm" )
print (f"Translation: {reg.transformation[:3 ,3 ]*100 } cm" )
Troubleshooting
Open3D visualization doesn't open on headless server
import open3d as o3d
o3d.visualization.rendering.OffscreenRenderer
pts = np.asarray(pcd.points)
import matplotlib.pyplot as plt
fig = plt.figure(); ax = fig.add_subplot(111 , projection='3d' )
ax.scatter(pts[::10 ,0 ], pts[::10 ,1 ], pts[::10 ,2 ], s=1 )
plt.savefig("pcd.png" )
Poisson reconstruction creates artifacts Fix : Increase depth or add more points:
mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(
pcd, depth=9 , width=0 , scale=1.1 , linear_fit=True
)
densities_arr = np.asarray(densities)
keep = densities_arr > np.percentile(densities_arr, 15 )
mesh.remove_vertices_by_mask(~keep)
Error: AttributeError: module 'open3d' has no attribute 'geometry' pip install open3d --upgrade
Version Compatibility Package Tested versions Notes open3d 0.18 API stable; 0.18 added better Poisson API numpy 1.24, 1.26 Required by open3d scipy 1.11, 1.12 For supplementary spatial analysis
External Resources
Official Documentation
Key Papers
Schindler, G. & Dellaert, F. (2012). 4D cities . 3DV Workshop.
Opitz, R. & Limp, W.F. (2015). Recent developments in high-density survey and measurement (HDSM) for archaeology . Annual Review of Anthropology.
Examples
Example 1: Flat Surface Extraction for Profile Drawing import numpy as np
import matplotlib.pyplot as plt
pts = np.asarray(pcd.points)
slice_mask = np.abs (pts[:,1 ]) < 0.002
slice_pts = pts[slice_mask]
fig, ax = plt.subplots(figsize=(5 , 8 ))
ax.scatter(slice_pts[:,0 ]*100 , slice_pts[:,2 ]*100 , s=3 , c='brown' , alpha=0.7 )
ax.set_xlabel("X (cm)" ); ax.set_ylabel("Z/Height (cm)" )
ax.set_title("Vessel Cross-Section Profile\n(Y ≈ 0 plane)" )
ax.set_aspect('equal' ); ax.grid(True , alpha=0.3 )
plt.tight_layout(); plt.savefig("cross_section_profile.png" , dpi=150 ); plt.show()
Example 2: Comparing Two Vessel Forms import numpy as np
import matplotlib.pyplot as plt
def vessel_profile (belly_r, neck_r, height, n_pts=50 ):
z = np.linspace(0 , height, n_pts)
t = z / height
r = belly_r * np.sin(np.pi * t) + neck_r * (1 - np.sin(np.pi * t))
return z, r
fig, ax = plt.subplots(figsize=(7 , 7 ))
for (br, nr, h, label, color) in [
(12 , 5 , 25 , 'Wide-mouth jar' , 'blue' ),
(8 , 7 , 20 , 'Storage jar' , 'red' ),
]:
z, r = vessel_profile(br, nr, h)
ax.plot(r, z, color=color, linewidth=2.5 , label=label)
ax.plot(-r, z, color=color, linewidth=2.5 )
ax.fill_betweenx(z, -r, r, alpha=0.1 , color=color)
ax.axhline(0 , color='gray' , linewidth=0.5 )
ax.set_xlabel("Radius (cm)" ); ax.set_ylabel("Height (cm)" )
ax.set_title("Comparative Vessel Profile Analysis" ); ax.legend()
ax.set_aspect('equal' ); ax.grid(True , alpha=0.3 )
plt.tight_layout(); plt.savefig("vessel_comparison.png" , dpi=150 ); plt.show()
Last updated: 2026-03-17 | Maintainer: @xjtulyc
Issues: GitHub Issues
