[Solved] Intersections between shell and lines (or other geometric objects)

Need help, or want to share a macro? Post here!
Forum rules
Be nice to others! Respect the FreeCAD code of conduct!
Post Reply
User avatar
onekk
Veteran
Posts: 6146
Joined: Sat Jan 17, 2015 7:48 am
Contact:

[Solved] Intersections between shell and lines (or other geometric objects)

Post by onekk »

I'm facing a problem:
inters_f1.png
inters_f1.png (5.76 KiB) Viewed 1684 times
inters_f2.png
inters_f2.png (5.81 KiB) Viewed 1684 times
This simple problem is very difficult to sort out.

I have these lines, that are axis of solids, (see post https://forum.freecadweb.org/viewtopic.php?f=22&t=65308 for the real thing)

I have to obtain some planes to cut these solids, at a predefined distance, determined by the shell. (Real problem have a number of concentric shells that have to be used derive profiles, for a complex loft).

I have thought that it will be simple to obtain intersections from a shell and a line, but it is not.

see the code:
20220127-intersect.py
(4.43 KiB) Downloaded 116 times

Code: Select all

    int_pts = []
    for l_idx in range(0, n_sp):
        for face in shell.Faces:
            points, curves = s_lines[l_idx].intersectCS(face.Surface)
            if len(points) != 0:
                int_pts.append((l_idx, points[0]))

    for p_idx, points in enumerate(int_pts):
        print(points)
Function intersectCS is expecting a surface, so I take face.Surface property (each shell has multiple faces).

But it seems that Surface is using the whole surface, so in my case the two circles and the two infinite lines, matching point outside the "real Shell".
  • There is a way to avoid this?
  • There is a technique to decimate obtained points to keep only real points. (I could think of checking if a point is in contact with the Shell, or similar approach)
Some consideration about speed would be appreciated, as in complex cases, as in complex cases it will be very useful to have a fast way.

TIA and Regards

Carlo D.
Last edited by onekk on Sat Jan 29, 2022 6:47 pm, edited 1 time in total.
GitHub page: https://github.com/onekk/freecad-doc.
- In deep articles on FreeCAD.
- Learning how to model with scripting.
- Various other stuffs.

Blog: https://okkmkblog.wordpress.com/
User avatar
Chris_G
Veteran
Posts: 2579
Joined: Tue Dec 31, 2013 4:10 pm
Location: France
Contact:

Re: Intersections between shell and lines (or other geometric objects)

Post by Chris_G »

Here is another solution, working on topology :

Code: Select all

diag_len = shell.BoundBox.DiagonalLength
edges = []
for line in s_lines:
    edges.append(line.toShape(0.0, diag_len))  # we assume lines are starting somewhere inside the shell
comp = Part.Compound(edges)
dist, ipts, info = comp.distToShape(shell)
# if dist == 0.0 (real intersection where found)
# then the 2 lists of points below will be identical
points_on_lines, points_on_shell = list(zip(*ipts))
Part.show(Part.Compound([Part.Vertex(p) for p in points_on_lines]))
Part.show(Part.Compound([Part.Vertex(p) for p in points_on_shell]))
All edges are grouped into a compound in order to perform a single distToShape operation.
But if desired, you can also do it on each edge individually.
Attachments
20220127-intersect_distToShape.py
(4.5 KiB) Downloaded 21 times
User avatar
onekk
Veteran
Posts: 6146
Joined: Sat Jan 17, 2015 7:48 am
Contact:

Re: Intersections between shell and lines (or other geometric objects)

Post by onekk »

Chris_G wrote: Thu Jan 27, 2022 3:07 pm Here is another solution, working on topology :
...
Many Thanks, @Chris_G

Very briiliant solution.

I have to dig in the definition of distToShape() to full see the working of your solution.

So for the casual user I will report here a copy of the documention found in the Python Console.
Find the minimum distance to another shape.

distToShape(shape) -> (dist, vectors, infos)

--

dist is the minimum distance, in mm (float value).

vectors is a list of pairs of App.Vector. Each pair corresponds to solution.


Example: [(Vector (2.0, -1.0, 2.0), Vector (2.0, 0.0, 2.0)), (Vector (2.0,-1.0, 2.0), Vector (2.0, -1.0, 3.0))]

First vector is a point on self

Second vector is a point on s.

infos contains additional info on the solutions. It is a list of tuples:

(topo1, index1, params1, topo2, index2, params2)

topo1, topo2 are strings identifying type of BREP element:

'Vertex',
'Edge',
or 'Face'.

index1, index2 are indexes of the elements (zero-based).

params1, params2 are parameters of internal space of the elements.

For vertices, params is None.
For edges, params is one float, u.
For faces, params is a tuple (u,v)
I have to digest the returned infos distToShape().

In the project above I have to put a sectionings plane orthogonal to each line at intersection point, to section properly the shape.

I could use some information returned from infos?

If I guess correctly the returned infos are:
  • the line parameter, so the distance from line start
  • the "face index" on which the line intersect the shell (u, v) are some "black magic" for me, (could you advise me some beginner text to learn something about them).
If you would and want spend some more time in helping me.

However many, many thanks again.

Carlo D.
GitHub page: https://github.com/onekk/freecad-doc.
- In deep articles on FreeCAD.
- Learning how to model with scripting.
- Various other stuffs.

Blog: https://okkmkblog.wordpress.com/
edwilliams16
Veteran
Posts: 3109
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: Intersections between shell and lines (or other geometric objects)

Post by edwilliams16 »

This was interesting to me, so I poked around, adding some diagnostics to the code.

Code: Select all

"""Sample code.

This code was written as an sample code

Name: filename.py

Author: Carlo Dormeletti
Copyright: 2022
Licence: CC BY-NC-ND 4.0 IT
"""
import os # noqa
import sys  # noqa
import math # noqa
from math import cos, pi, sin, sqrt  # noqa

import FreeCAD
import FreeCADGui # noqa
from FreeCAD import Placement, Rotation, Vector # noqa
import Part # noqa

DOC = FreeCAD.activeDocument()

DOC_NAME = "test_offset"


def clear_doc():
    """Clear active document deleting all the objects."""
    for obj in DOC.Objects:
        DOC.removeObject(obj.Name)


def setview():
    """Rearrange View."""
    FreeCAD.Gui.SendMsgToActiveView("ViewFit")
    FreeCAD.Gui.activeDocument().activeView().viewAxometric()


if DOC is None:
    FreeCAD.newDocument(DOC_NAME)
    FreeCAD.setActiveDocument(DOC_NAME)
    DOC = FreeCAD.activeDocument()

else:

    clear_doc()

# EPS= tolerance to use to cut the parts
EPS = 0.10
EPS_C = EPS * -0.5
ROT0 = Rotation(0, 0, 0)

# CODE START HERE


def tangent2C(c1, c2, rad1, rad2):
    """Calculate four tangent point to join two circles.

    Courtesy of edwilliams16

    Parameters:
    c1    Vector  Center of circle1
    c2    Vector  Center of circle2
    rad1  float   Radius of circle1
    rad2  float   Radius of circle2

    Return:
    [tg1, tg2, tg3, tg4] tangent points.
    If one circle strictly encloses the other or are they are concentric
    return []
    """
    lg = (c1 - c2).Length

    if abs(rad1 - rad2) > lg or lg == 0:
        return []
    else:
        c = (rad1 - rad2) / lg
        s = sqrt(1 - c * c)
        n = (c2 - c1) / lg  # unit vector along line of centers
        m = Vector(0, 0, 1).cross(n)  # orthogonal vector
        vp = c * n + s * m  # unit vectors from center to tangent points
        vm = c * n - s * m
        return [c1 + rad1 * vp, c1 + rad1 * vm, c2 + rad2 * vp, c2 + rad2 * vm]


def dest_pt(pt, angle, length):
    """Calculate destination point.

    Vector math courtesy of edwilliams16

    Parameters:
      pt    = starting point
      angle =  rad
      length = units
    """
    # return Vector(pt.x + math.cos(angle) * length, pt.y + math.sin(angle) * length, pt.z)

    lv = Vector(length, 0, 0)
    p_rot = Rotation(math.degrees(angle), 0, 0)
    dept = Placement(pt, p_rot).multVec(lv)

    return dept


def mk_shell():
    """Make an shell for test."""
    c1cen = Vector(0, 0, 0)
    c1rad = 100
    c2cen = Vector(0, 100, 0)
    c2rad = 50

    tgp1, tgp2, tgp3, tgp4 = tangent2C(c1cen, c2cen, c1rad, c2rad)

    cir1 = Part.makeCircle(c1rad, c1cen)
    cir2 = Part.makeCircle(c2rad, c2cen)
    edg1 = Part.LineSegment(tgp3, tgp1).toShape()
    edg2 = Part.LineSegment(tgp2, tgp4).toShape()

    # obtain cir1 section between tangents
    par1 = cir1.Curve.parameter(tgp1)
    par2 = cir1.Curve.parameter(tgp2)

    tgv1 = cir1.Curve.value(par1)
    tgv2 = cir1.Curve.value(par2)

    c1_3pos = tgv1 - (tgv2 - tgv1) * 0.5
    par3 = cir1.Curve.parameter(c1_3pos)

    tgv3 = cir1.Curve.value(par3)

    edg3 = Part.ArcOfCircle(tgv2, tgv3, tgv1).toShape()

    # obtain cir2 section between tangents
    par4 = cir2.Curve.parameter(tgp3)
    par5 = cir2.Curve.parameter(tgp4)

    tgv4 = cir2.Curve.value(par4)
    tgv5 = cir2.Curve.value(par5)

    c2_3pos = tgv4 + (tgv5 - tgv4) * 0.5
    par6 = cir2.Curve.parameter(c2_3pos)

    tgv6 = cir2.Curve.value(par6)

    edg4 = Part.ArcOfCircle(tgv4, tgv6, tgv5).toShape()

    # edge are ordered in counterclockwise
    obj_f = Part.Wire((edg2, edg3, edg1, edg4))

    # Part.show(obj_f, "wire")

    return obj_f


def find_inters():
    """Find intersection between lines and a shell."""
    n_sp = 21
    ang_adv = 360 / n_sp
    ln_lg = 200
    sh_thi = 80

    s_lines = []
    for idx in range(0, n_sp):
        end_pt = dest_pt(
            Vector(0, 0, 0),
            pi / 2.0 + math.radians(ang_adv * idx),
            ln_lg)
        line = Part.LineSegment(
            Vector(0, 0, sh_thi * 0.5),
            Vector(end_pt.x, end_pt.y, sh_thi * 0.5))
        s_lines.append(line)
        Part.show(line.toShape(), "line{}".format(idx))

    wire = mk_shell()
    shell = wire.extrude(Vector(0, 0, sh_thi))
    Part.show(shell, "shell")

    diag_len = shell.BoundBox.DiagonalLength
    edges = []
    for line in s_lines:
        edges.append(line.toShape(0.0, diag_len))
    comp = Part.Compound(edges)
    d, ipts, info = comp.distToShape(shell)
    #print(f'd[0]= {d[0]} ipts[0]= {ipts[0]} info = {info[0]}')
    points_on_lines, points_on_shell = list(zip(*ipts))
    Part.show(Part.Compound([Part.Vertex(p) for p in points_on_lines]))
    Part.show(Part.Compound([Part.Vertex(p) for p in points_on_shell]))
    return (d, ipts, info, comp)

d, ipts, info, comp = find_inters()

setview()

#look at lineN intersection with the shell
N = 5
lN = App.ActiveDocument.getObject('line'+str(N)).Shape
shl= App.ActiveDocument.getObject('shell').Shape
shp=App.ActiveDocument.getObject('Shape').Shape
infoN = info[N]
print(f'Intersection from Shape list at {shp.Vertexes[N].Point}')
print(f'Intersection on line at {lN.valueAt(infoN[2])}')
face = shl.Faces[infoN[4]]
faceuv = infoN[5]  #u,v coordinate of point on the shell face
print(f'Intersection on shell Face at {face.valueAt(*faceuv)}')
edge = comp.Edges[infoN[1]]
edgeu = infoN[2] # u coordinate along edge
print(f'Intersection on edge at {edge.valueAt(edgeu)}')
# in addition to valueAt() we have normalAt() froviding a normal vector
# and derivative1At() providing tangent vectors
The info return locates the intersections (Actually nearest points) on the two shapes - in this case lines and a shell, giving the coordinate(s) on the edge or face.
The sample code shows the numerous ways of extracting the location of the point. By the same route, you can find the normal and tangent vectors at the point, which you can use to construct further geometry.

There were two puzzles that came up:
(1) d was returned a single scalar. Shouldn't is be a list of distances for the items in the compound?
(2) In the case of the face, derivative2At() provides (I think) d2x/du^2 and d2x/dv^2. If this is so, the cross-derivative d2x/dudv is missing, and would be required to characterize the surface curvature.

I've not followed up on these yet.
User avatar
Chris_G
Veteran
Posts: 2579
Joined: Tue Dec 31, 2013 4:10 pm
Location: France
Contact:

Re: Intersections between shell and lines (or other geometric objects)

Post by Chris_G »

edwilliams16 wrote: Thu Jan 27, 2022 9:40 pm (1) d was returned a single scalar. Shouldn't is be a list of distances for the items in the compound?
The distToShape help says :
Find the minimum distance to another shape
So, logically, this can only be a single value.
But this value (especially when there is a real contact, and d=0) can happen in multiple locations.
This is why d is always a single float, but points and info are lists of solutions.
edwilliams16 wrote: Thu Jan 27, 2022 9:40 pm (2) In the case of the face, derivative2At() provides (I think) d2x/du^2 and d2x/dv^2. If this is so, the cross-derivative d2x/dudv is missing, and would be required to characterize the surface curvature.
The surface can provide that kind of information :

Code: Select all

aFace.Surface.getDN(u, v, u_order, v_order)
So the dudv derivative will be :

Code: Select all

aFace.Surface.getDN(u, v, 1, 1)
Post Reply