[3d] 3d 구를 DSM으로 만들기

TIN to DSM

TIN? Triangular Irregular Network

DSM? Digital Surface Model

보통 지형 모델을 2.5D TIN으로 만든 후 이를 DSM으로 만든다. 이 작업을 구현했는데 Terratin(지형)으로 바로하기엔 고난이도라 간단한 구로 먼저 해보았다.

 

In [1]:

from pyhull.convex_hull import ConvexHull
from pyhull.delaunay import DelaunayTri
from Delaunator.Delaunator import Delaunator
from time import time
import pyvista

In [2]:

import parmap
import numpy as np

create grid¶

In [212]:

mesh_sphere = o3d.geometry.TriangleMesh.create_sphere(radius=5)

In [213]:

o3d.visualization.draw_geometries([mesh_sphere],mesh_show_wireframe=True)

In [214]:

bbox = mesh_sphere.get_axis_aligned_bounding_box()

In [215]:

RB = bbox.max_bound; LT = bbox.min_bound

In [216]:

RB[1] - LT[0]

Out[216]:

10.0

In [217]:

import math
pixellength = 0.5
width = RB[0] - LT[0]
height = RB[1] - LT[1]
x = np.linspace(LT[0], RB[0], math.ceil(width/pixellength))
y = np.linspace(LT[1], RB[1], math.ceil(height/pixellength))

In [218]:

xx, yy = np.meshgrid(x, y)

In [220]:

xyz = np.concatenate([xx[:,:, np.newaxis], yy[:,:, np.newaxis], np.zeros_like(xx)[:, :, np.newaxis]], axis=-1)

In [221]:

xyz2d = np.concatenate(xyz)

In [222]:

xyz2d.shape

Out[222]:

(400, 3)

point cloud and faces data¶

In [223]:

pcd_sphere = np.asarray(mesh_sphere.vertices)

In [224]:

face_sphere = np.asarray(mesh_sphere.triangles)

projected terratin on z = 0¶

In [225]:

tmesh = trimesh.Trimesh(pcd_sphere, face_sphere)

In [226]:

tmesh_projectioned = trimesh.Trimesh(np.concatenate([pcd_sphere[:, :2], np.zeros_like(pcd_sphere[:, :1])], axis=1), face_sphere)

In [227]:

tmesh_query = trimesh.proximity.ProximityQuery(tmesh_projectioned)

find points on the triangle¶

In [229]:

candidates = tmesh_query.on_surface(np.array(xyz2d, dtype=np.float32))

if signed distance < 0 -> out of the mesh, write 0¶

In [231]:

drawed_region = tmesh_query.signed_distance(xyz2d) >= 0

In [232]:

drawed_faces = candidates[-1][drawed_region]
intersect_pts = candidates[0][drawed_region]

find normal_Vector of original mesh¶

In [233]:

normals = tmesh.face_normals[candidates[2]]
normal_candidates = normals[drawed_region]
normal_candidates.shape
# [candidates[2]]

Out[233]:

(276, 3)

find pt on the original mesh plane¶

In [234]:

plane_origin = tmesh.triangles_center[candidates[2]][drawed_region]

In [235]:

plane_origin.shape

Out[235]:

(276, 3)

find endpoints of lines¶

In [236]:

temp = xyz2d[drawed_region]
temp[:, 2] = 1

In [237]:

endpoints = np.concatenate([xyz2d[drawed_region][np.newaxis, :, :], temp[np.newaxis, :, :]])

find z coordinates of which intersection point lies onn the plane (intersetion btw line and plane)¶

In [238]:

intersections = []
for i in range(len(plane_origin)):
    it = trimesh.intersections.plane_lines(plane_origin[i], normal_candidates[i], endpoints[:, i:i+1, :], line_segments=False)
    intersections.append(it)

In [239]:

intersections_pts = np.asarray(list(map(lambda x:x[0][0], intersections)))

draw DSM¶

In [256]:

from copy import deepcopy

In [257]:

dsm = deepcopy(xyz2d)

In [258]:

temp2 = dsm[drawed_region]

In [259]:

temp2[:, 2] = intersections_pts[:, 2]

In [260]:

dsm[drawed_region] = temp2

replace empty zone to 0 pixel¶

In [261]:

empty_zone= dsm[~drawed_region]
empty_zone[:, 2] = 0

In [262]:

dsm[~drawed_region] = empty_zone

dsm reshape to w,h,3 for image processing¶

In [263]:

dsm = dsm.reshape(xyz.shape)

other z pixel value smaller than 0 makes 0 pixel¶

In [264]:

mask = dsm[:, :, 2] < 0

In [265]:

temp = dsm[mask]
temp[:, 2] = np.abs(temp[:, 2])
dsm[mask] = temp

In [266]:

dsm_last = dsm[:, :, 2][:, :, np.newaxis]

In [267]:

import cv2

In [268]:

color_img = cv2.normalize(dsm_last, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)
print(color_img.shape, color_img.dtype)

# convert to 3 channel
color_img = cv2.cvtColor(color_img, cv2.COLOR_GRAY2BGR)
print(color_img.shape, color_img.dtype)
print(np.amin(color_img),np.amax(color_img))

(20, 20) uint8
(20, 20, 3) uint8
0 255

In [269]:

cv2.imwrite('color_img.png', color_img)

Out[269]:

True