There is No intersect point returned by Mesh.intersect_with_line().

[ZiguoAtGitHub]

Hi Marco,
I want to get the intersection point of a ray with a mesh by using Mesh.intersect_with_line() function.
However sometimes I can get the intersection point, sometimes I can not like this.

Here is my code. The data file is also uploaded. Looking forward to your reply.🙏🏻

from vedo import *
import re
import math as mmath
from scipy.interpolate import interp1d

#------------------------------------------------------------------------------------------Begin reading PMM Result File
fin=open(file='A000_pmm.txt',mode='r')
tlines=fin.readlines()
fin.close()
n_tlines=len(tlines)
Angle_N_M_k={}
Angle_N_M_c={}
iline=0
while (iline < n_tlines):
  if(tlines[iline].find('极限强度面')>1):
    jline=iline+2
    tlist=[iv for iv in re.split('[, ]+',tlines[jline].strip()) if len(iv)>3]
    for iangle,anglei in enumerate(tlist):Angle_N_M_k[iangle]=[float(anglei),[],[]]
    jline=jline+2
    nm=[]
    while(tlines[jline].find('设计强度面')<0):
      ttlist=[float(ival) for ival in re.split('[,]',tlines[jline])[:-1]]
      for iangle in Angle_N_M_k:
        Angle_N_M_k[iangle][2].append(ttlist[iangle*2+1])
        Angle_N_M_k[iangle][1].append(ttlist[iangle*2])
      jline=jline+1
    iline=jline
  if(tlines[iline].find('设计强度面')>1):
    jline=iline+1
    tlist=[iv for iv in re.split('[, ]+',tlines[jline].strip()) if len(iv)>3]
    for iangle,anglei in enumerate(tlist):Angle_N_M_c[iangle]=[float(anglei),[],[]]
    jline=jline+2
    while(jline<n_tlines):
      NM_list=tlines[jline].split(',')[:-1]
      for k in Angle_N_M_c:
        Angle_N_M_c[k][1].append(float(NM_list[2*k]))
        Angle_N_M_c[k][2].append(float(NM_list[2*k+1]))
      jline=jline+1
    iline=jline
if len(Angle_N_M_c)==1:
  for i in range(1,10):
    Angle_N_M_c[i]=[i*10.0,Angle_N_M_c[0][1],Angle_N_M_c[0][2]]
    Angle_N_M_k[i]=[i*10.0,Angle_N_M_k[0][1],Angle_N_M_k[0][2]]
#------------------------------------------------------------------------------------------Reading PMM Result File Finished
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Generate Fitting Functions based on PMM Result
if(len(Angle_N_M_k)>0): #--------------------------------Column & Brace PMMk
  itemp=10
  dist_lmt=abs(Angle_N_M_k[0][1][0]-Angle_N_M_k[0][1][-1])*0.001 #最大与最小轴力差的1/1000;Angle_N_M_k[iangle]=[angle,[N],[M]]
  new_Angle_N_M_k={}
  for iangle in Angle_N_M_k:
    new_Angle_N_M_k[iangle]=[Angle_N_M_k[iangle][0],[],[]]
    pts=[[Angle_N_M_k[iangle][1][0]],[Angle_N_M_k[iangle][2][0]]]				
    for irow in range(len(Angle_N_M_k[iangle][1])):
      dist=dist_lmt*5
      for ipt in range(len(pts[0])):dist=min(dist,mmath.sqrt((Angle_N_M_k[iangle][1][irow]-pts[0][ipt])**2+(Angle_N_M_k[iangle][2][irow]-pts[1][ipt])**2))
      if(dist>dist_lmt):
        pts[0].append(Angle_N_M_k[iangle][1][irow])
        pts[1].append(Angle_N_M_k[iangle][2][irow])
    new_Angle_N_M_k[iangle]=[Angle_N_M_k[iangle][0],pts[0],pts[1]]
    MaxNk=max(new_Angle_N_M_k[0][1])
    MinNk=min(new_Angle_N_M_k[0][1])
  #-------------------------fit the function of anglei
  fitfunc_k=[] #[(angle1,fitfunc1),(angle2,fitfunc2),........]
  for iangle in new_Angle_N_M_k:fitfunc_k.append((new_Angle_N_M_k[iangle][0],interp1d(new_Angle_N_M_k[iangle][1],new_Angle_N_M_k[iangle][2],kind='quadratic',fill_value="extrapolate")))
if(len(Angle_N_M_c)>=1): #--------------------------------Column & Brace PMMc
  dist_lmt=abs(Angle_N_M_c[0][1][0]-Angle_N_M_c[0][1][-1])*0.001
  new_Angle_N_M_c={}
  for iangle in Angle_N_M_c:
    new_Angle_N_M_c[iangle]=[Angle_N_M_c[iangle][0],[],[]]
    pts=[[Angle_N_M_c[iangle][1][0]],[Angle_N_M_c[iangle][2][0]]]				
    for irow in range(len(Angle_N_M_c[iangle][1])):
      dist=dist_lmt*5
      for ipt in range(len(pts[0])):dist=min(dist,mmath.sqrt((Angle_N_M_c[iangle][1][irow]-pts[0][ipt])**2+(Angle_N_M_c[iangle][2][irow]-pts[1][ipt])**2))
      if(dist>dist_lmt):
        pts[0].append(Angle_N_M_c[iangle][1][irow])
        pts[1].append(Angle_N_M_c[iangle][2][irow])
    new_Angle_N_M_c[iangle]=[Angle_N_M_c[iangle][0],pts[0],pts[1]]
    MaxNc=max(new_Angle_N_M_c[0][1])
    MinNc=min(new_Angle_N_M_c[0][1])
  #-------------------------fit the function of anglei
  fitfunc_c=[]#interpolate function of N-M
  for iangle in new_Angle_N_M_c:fitfunc_c.append((new_Angle_N_M_c[iangle][0],interp1d(new_Angle_N_M_c[iangle][1],new_Angle_N_M_c[iangle][2],kind='quadratic',fill_value="extrapolate")))
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Fitting Functions based on PMM Result Finished
#~~
#------------------------------------------------------------------------------------------Create Meshes
tN=np.linspace(MinNk,MaxNk,51)
tpts=[]
cmeshs=[[],[]]
tangles=len(fitfunc_k)
if tangles>2:
  #---------------------------------Create PMM Surface Mesh
  for iangle,ifunc in fitfunc_k:
    cmeshs[0]+=[(mmath.cos(mmath.radians(iangle))*ifunc(tN[inz])[None][0],mmath.sin(mmath.radians(iangle))*ifunc(tN[inz])[None][0],tN[inz]) for inz in range(51)]
    tpts.append([(mmath.cos(mmath.radians(iangle))*ifunc(tN[inz])[None][0],mmath.sin(mmath.radians(iangle))*ifunc(tN[inz])[None][0],tN[inz]) for inz in range(51)])
  for ipts in tpts[-2::-1]:cmeshs[0]+=[(-iipt[0],iipt[1],iipt[2]) for iipt in ipts]
  for ipts in tpts[1:]:cmeshs[0]+=[(-iipt[0],-iipt[1],iipt[2]) for iipt in ipts]
  for ipts in tpts[-2:0:-1]:cmeshs[0]+=[(iipt[0],-iipt[1],iipt[2]) for iipt in ipts]
  #~
  for i in range((tangles-1)*4):
    j=0 if i==(tangles-1)*4-1 else i+1
    cmeshs[1]+=[(i*51+k,i*51+k+1,j*51+k+1) for k in range(50)]
    cmeshs[1]+=[(j*51+k,i*51+k,j*51+k+1) for k in range(50)]
  #~~
  shs=[]
  nl_g=[4242.1,5124,16191.3]
  nl_pt=[iv*15 for iv in [2238.3,4709,212.5]]
  pmm=Mesh(cmeshs,c='cyan').alpha(0.2).lw(0.2).lc('k')
  #---------------------------------PMM Surface Mesh Created
  shs.append(pmm)
  nl_gps = Points([nl_g], r=10, c='g')
  nl_eps=Points([nl_pt],r=10,c='r').render_points_as_spheres(True)
  ln = Line([0,0,0], nl_g, c='blue')
  ln2 = Line(nl_g,nl_pt, c='red')
  shs.append(nl_gps)
  shs.append(nl_eps)
  shs.append(ln)
  shs.append(ln2)
  if pmm.is_inside(nl_g)==1:# nl_g is inside pmm
    int_pts = pmm.intersect_with_line(nl_g,nl_pt)
    intps = Points(int_pts, r=10, c='b')
    shs.append(int_pts)
    print(int_pts)
  show(shs, bg='white',axes=7).close()

A000_pmm.txt

[marcomusy]

Hi Ziguo, unfortunately I cannot run your script, I get:

Traceback (most recent call last):
  File "/home/musy/Projects/vedo/v5.py", line 8, in <module>
    tlines = fin.readlines()
  File "/home/musy/soft/miniconda3/lib/python3.10/codecs.py", line 322, in decode
    (result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xbc in position 17: invalid start byte

what format is A000_pmm.txt?

[ZiguoAtGitHub]

Hi Marco,
The problem should be caused by those chinese characters in the file.
The file codec is ANSI/OEM GBK.
I am not familiar with file encoding/decoding related knowledge but I found one solution by searching on web that is using 'ignore' option when decoding, like str.decode('utf-8','ignore'). Hope the problem could be solved.
To be honest, it is not a very tough work to get the intersection of a ray with plane in space and should not bother you.
Maybe I should take a try first.😀
Your Response is so quickly always. Thank you again.

[ZiguoAtGitHub]

Hi Marco,
Here is my solution:
😊

  def findIntersection(Ray, PolyMesh):
    tpoints=PolyMesh.points()
    ray0=np.array(Ray[0])
    ray1=np.array(Ray[1])
    ray_n1=ray1-ray0
    int_pt=None
    for itriangle in PolyMesh.faces():
      p0,p1,p2=[tpoints[itriangle[iv]] for iv in range(3)]
      plane_normal=np.cross(p1-p0,p2-p0)
      tdot=np.dot(plane_normal,ray_n1)
      if abs(tdot) !=0:
        ray_scalar=np.dot(p0-ray0,plane_normal)/np.dot(ray_n1,plane_normal)
        int_pt=ray0+ray_scalar*ray_n1
        # check the intersection falls inside the itriangle
        area0=np.linalg.norm(plane_normal) #itriangle area
        area1=np.linalg.norm(np.cross(p0-int_pt,p1-int_pt))+np.linalg.norm(np.cross(p1-int_pt,p2-int_pt))+np.linalg.norm(np.cross(p2-int_pt,p0-int_pt))
        if abs(area1-area0)/area0<0.01 and np.dot(int_pt-ray0,ray_n1)>0:return int_pt
  #~~~~

[marcomusy]

Hi Ziguo sorry I could not run it (i also tried to change characters manually but the program hangs)... anyways, consider that you can play with tol=... in intersect_with_line()

Also consider:
vedo.utils.intersection_ray_triangle(P0, P1, V0, V1, V2)

[ZiguoAtGitHub]

OK, Thanks.