Commit 809b4f67 by Falko

### improve triangle to point distance caclulation, improve plotting

parent 9140a4ce
 # Tests distance between point and triangle in 3D. Aligns and uses 2D technique. # # Was originally some code on mathworks from numpy import dot from math import sqrt def pointTriangleDistance(TRI, P): # function [dist,PP0] = pointTriangleDistance(TRI,P) # calculate distance between a point and a triangle in 3D # SYNTAX # dist = pointTriangleDistance(TRI,P) # [dist,PP0] = pointTriangleDistance(TRI,P) # # DESCRIPTION # Calculate the distance of a given point P from a triangle TRI. # Point P is a row vector of the form 1x3. The triangle is a matrix # formed by three rows of points TRI = [P1;P2;P3] each of size 1x3. # dist = pointTriangleDistance(TRI,P) returns the distance of the point P # to the triangle TRI. # [dist,PP0] = pointTriangleDistance(TRI,P) additionally returns the # closest point PP0 to P on the triangle TRI. # # Author: Gwolyn Fischer # Release: 1.0 # Release date: 09/02/02 # Release: 1.1 Fixed Bug because of normalization # Release: 1.2 Fixed Bug because of typo in region 5 20101013 # Release: 1.3 Fixed Bug because of typo in region 2 20101014 # Possible extention could be a version tailored not to return the distance # and additionally the closest point, but instead return only the closest # point. Could lead to a small speed gain. # Example: # %% The Problem # P0 = [0.5 -0.3 0.5] # # P1 = [0 -1 0] # P2 = [1 0 0] # P3 = [0 0 0] # # vertices = [P1; P2; P3] # faces = [1 2 3] # # %% The Engine # [dist,PP0] = pointTriangleDistance([P1;P2;P3],P0) # # %% Visualization # [x,y,z] = sphere(20) # x = dist*x+P0(1) # y = dist*y+P0(2) # z = dist*z+P0(3) # # figure # hold all # patch('Vertices',vertices,'Faces',faces,'FaceColor','r','FaceAlpha',0.8) # plot3(P0(1),P0(2),P0(3),'b*') # plot3(PP0(1),PP0(2),PP0(3),'*g') # surf(x,y,z,'FaceColor','b','FaceAlpha',0.3) # view(3) # The algorithm is based on # "David Eberly, 'Distance Between Point and Triangle in 3D', # Geometric Tools, LLC, (1999)" # http:\\www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf # # ^t # \ | # \reg2| # \ | # \ | # \ | # \| # *P2 # |\ # | \ # reg3 | \ reg1 # | \ # |reg0\ # | \ # | \ P1 # -------*-------*------->s # |P0 \ # reg4 | reg5 \ reg6 # rewrite triangle in normal form B = TRI[0, :] E0 = TRI[1, :] - B # E0 = E0/sqrt(sum(E0.^2)); %normalize vector E1 = TRI[2, :] - B # E1 = E1/sqrt(sum(E1.^2)); %normalize vector D = B - P a = dot(E0, E0) b = dot(E0, E1) c = dot(E1, E1) d = dot(E0, D) e = dot(E1, D) f = dot(D, D) #print "{0} {1} {2} ".format(B,E1,E0) det = a * c - b * b s = b * e - c * d t = b * d - a * e # Terible tree of conditionals to determine in which region of the diagram # shown above the projection of the point into the triangle-plane lies. if (s + t) <= det: if s < 0.0: if t < 0.0: # region4 if d < 0: t = 0.0 if -d >= a: s = 1.0 sqrdistance = a + 2.0 * d + f else: s = -d / a sqrdistance = d * s + f else: s = 0.0 if e >= 0.0: t = 0.0 sqrdistance = f else: if -e >= c: t = 1.0 sqrdistance = c + 2.0 * e + f else: t = -e / c sqrdistance = e * t + f # of region 4 else: # region 3 s = 0 if e >= 0: t = 0 sqrdistance = f else: if -e >= c: t = 1 sqrdistance = c + 2.0 * e + f else: t = -e / c sqrdistance = e * t + f # of region 3 else: if t <= 0: # region 5 t = 0 if d >= 0: s = 0 sqrdistance = f else: if -d >= a: s = 1 sqrdistance = a + 2.0 * d + f; # GF 20101013 fixed typo d*s ->2*d else: s = -d / a sqrdistance = d * s + f else: # region 0 invDet = 1.0 / det s = s * invDet t = t * invDet sqrdistance = s * (a * s + b * t + 2.0 * d) + t * (b * s + c * t + 2.0 * e) + f else: if s < 0.0: # region 2 tmp0 = b + d tmp1 = c + e if tmp1 > tmp0: # minimum on edge s+t=1 numer = tmp1 - tmp0 denom = a - 2.0 * b + c if numer >= denom: s = 1.0 t = 0.0 sqrdistance = a + 2.0 * d + f; # GF 20101014 fixed typo 2*b -> 2*d else: s = numer / denom t = 1 - s sqrdistance = s * (a * s + b * t + 2 * d) + t * (b * s + c * t + 2 * e) + f else: # minimum on edge s=0 s = 0.0 if tmp1 <= 0.0: t = 1 sqrdistance = c + 2.0 * e + f else: if e >= 0.0: t = 0.0 sqrdistance = f else: t = -e / c sqrdistance = e * t + f # of region 2 else: if t < 0.0: # region6 tmp0 = b + e tmp1 = a + d if tmp1 > tmp0: numer = tmp1 - tmp0 denom = a - 2.0 * b + c if numer >= denom: t = 1.0 s = 0 sqrdistance = c + 2.0 * e + f else: t = numer / denom s = 1 - t sqrdistance = s * (a * s + b * t + 2.0 * d) + t * (b * s + c * t + 2.0 * e) + f else: t = 0.0 if tmp1 <= 0.0: s = 1 sqrdistance = a + 2.0 * d + f else: if d >= 0.0: s = 0.0 sqrdistance = f else: s = -d / a sqrdistance = d * s + f else: # region 1 numer = c + e - b - d if numer <= 0: s = 0.0 t = 1.0 sqrdistance = c + 2.0 * e + f else: denom = a - 2.0 * b + c if numer >= denom: s = 1.0 t = 0.0 sqrdistance = a + 2.0 * d + f else: s = numer / denom t = 1 - s sqrdistance = s * (a * s + b * t + 2.0 * d) + t * (b * s + c * t + 2.0 * e) + f # account for numerical round-off error if sqrdistance < 0: sqrdistance = 0 dist = sqrt(sqrdistance) PP0 = B + s * E0 + t * E1 return dist, PP0 \ No newline at end of file
No preview for this file type
 splipy k3d #k3d scipy ipyvolume seaborn \ No newline at end of file tqdm
 %% Cell type:code id: tags: ``` python import itertools import random from tqdm import tqdm import ipyvolume as ipv import ipywidgets as widgets import numpy as np import trimesh from scipy import spatial import matplotlib.pyplot as plt import seaborn as sns ``` from distance_triangle import pointTriangleDistance %% Cell type:code id: tags: num_base_points = 200 # number of voronoi cells resize_base_points_factor = 2 # scaling factor for base points of cells (so the boundary of the unit cube is uniformly splitted into voronoi cells) ``` python num_base_points = 1000 # number of voronoi cells resize_base_points_factor = 2 # scaling factor for base points of cells (so the boundary of the unit cube is uniformly split into voronoi cells) join_cell_ratio = 0.95 # the ratio of cells which are to be joined with another joined_cells_size_gap = 10 # the largest cell cluster is allowed 5 times the size of the smallest num_pixels = 25 # number of voxels per side of the unit cube noise = 0.1 # noise in the output voxels num_pixels = 100 # number of voxels per side of the unit cube noise = 0.04 # noise in the output voxels