源代码:
##Copyright 2009-2016 Jelle Feringa (jelleferinga@gmail.com)
##
##This file is part of pythonOCC.
##
##pythonOCC is free software: you can redistribute it and/or modify
##it under the terms of the GNU Lesser General Public License as published by
##the Free Software Foundation, either version 3 of the License, or
##(at your option) any later version.
##
##pythonOCC is distributed in the hope that it will be useful,
##but WITHOUT ANY WARRANTY; without even the implied warranty of
##MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
##GNU Lesser General Public License for more details.
##
##You should have received a copy of the GNU Lesser General Public License
##along with pythonOCC. If not, see .
from __future__ import print_functionimport os
import sys
import timefrom OCC.Core.BRep import BRep_Tool
from OCC.Core.BRepAdaptor import BRepAdaptor_HCurve
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakePolygon
from OCC.Core.BRepFill import BRepFill_CurveConstraint
from OCC.Display.SimpleGui import init_display
from OCC.Core.GeomAbs import GeomAbs_C0
from OCC.Core.GeomLProp import GeomLProp_SLProps
from OCC.Core.GeomPlate import (GeomPlate_BuildPlateSurface, GeomPlate_PointConstraint,GeomPlate_MakeApprox)
from OCC.Core.ShapeAnalysis import ShapeAnalysis_Surface
from OCC.Core.gp import gp_Pnt
from OCC.Core.BRepFill import BRepFill_Fillingfrom OCC.Extend.TopologyUtils import TopologyExplorer, WireExplorer
from OCC.Extend.ShapeFactory import make_face, make_vertex
from OCC.Extend.DataExchange import read_iges_filedisplay, start_display, add_menu, add_function_to_menu = init_display()try:HAS_SCIPY = Truefrom scipy.optimize import fsolve
except ImportError:print('scipy not installed, will not be able to run the geomplate example')HAS_SCIPY = False# TODO:
# - need examples where the tangency to constraining faces is respecteddef make_n_sided(edges, points, continuity=GeomAbs_C0):"""builds an n-sided patch, respecting the constraints defined by *edges*and *points*a simplified call to the BRepFill_Filling classits simplified in the sense that to all constraining edges and pointsthe same level of *continuity* will be applied*continuity* represents:GeomAbs_C0 : the surface has to pass by 3D representation of the edgeGeomAbs_G1 : the surface has to pass by 3D representation of the edgeand to respect tangency with the given faceGeomAbs_G2 : the surface has to pass by 3D representation of the edgeand to respect tangency and curvature with the given face.NOTE: it is not required to set constraining points.just leave the tuple or list empty:param edges: the constraining edges:param points: the constraining points:param continuity: GeomAbs_0, 1, 2:return: TopoDS_Face"""n_sided = BRepFill_Filling()for edg in edges:n_sided.Add(edg, continuity)for pt in points:n_sided.Add(pt)n_sided.Build()face = n_sided.Face()return facedef make_closed_polygon(*args):poly = BRepBuilderAPI_MakePolygon()for pt in args:if isinstance(pt, list) or isinstance(pt, tuple):for i in pt:poly.Add(i)else:poly.Add(pt)poly.Build()poly.Close()result = poly.Wire()return resultdef geom_plate(event=None):display.EraseAll()p1 = gp_Pnt(0, 0, 0)p2 = gp_Pnt(0, 10, 0)p3 = gp_Pnt(0, 10, 10)p4 = gp_Pnt(0, 0, 10)p5 = gp_Pnt(5, 5, 5)poly = make_closed_polygon([p1, p2, p3, p4])edges = [i for i in TopologyExplorer(poly).edges()]face = make_n_sided(edges, [p5])display.DisplayShape(edges)display.DisplayShape(make_vertex(p5))display.DisplayShape(face, update=True)# ============================================================================
# Find a surface such that the radius at the vertex is n
# ============================================================================def build_plate(polygon, points):'''build a surface from a constraining polygon(s) and point(s)@param polygon: list of polygons ( TopoDS_Shape)@param points: list of points ( gp_Pnt )'''# plate surfacebpSrf = GeomPlate_BuildPlateSurface(3, 15, 2)# add curve constraintsfor poly in polygon:for edg in WireExplorer(poly).ordered_edges():c = BRepAdaptor_HCurve()c.ChangeCurve().Initialize(edg)constraint = BRepFill_CurveConstraint(c, 0)bpSrf.Add(constraint)# add point constraintfor pt in points:bpSrf.Add(GeomPlate_PointConstraint(pt, 0))bpSrf.Perform()maxSeg, maxDeg, critOrder = 9, 8, 0tol = 1e-4dmax = max([tol, 10 * bpSrf.G0Error()])srf = bpSrf.Surface()plate = GeomPlate_MakeApprox(srf, tol, maxSeg, maxDeg, dmax, critOrder)uMin, uMax, vMin, vMax = srf.Bounds()return make_face(plate.Surface(), uMin, uMax, vMin, vMax, 1e-4)def radius_at_uv(face, u, v):'''returns the mean radius at a u,v coordinate@param face: surface input@param u,v: u,v coordinate'''h_srf = BRep_Tool().Surface(face)#uv_domain = GeomLProp_SurfaceTool().Bounds(h_srf)curvature = GeomLProp_SLProps(h_srf, u, v, 1, 1e-6)try:_crv_min = 1. / curvature.MinCurvature()except ZeroDivisionError:_crv_min = 0.try:_crv_max = 1. / curvature.MaxCurvature()except ZeroDivisionError:_crv_max = 0.return abs((_crv_min + _crv_max) / 2.)def uv_from_projected_point_on_face(face, pt):'''returns the uv coordinate from a projected point on a face'''srf = BRep_Tool().Surface(face)sas = ShapeAnalysis_Surface(srf)uv = sas.ValueOfUV(pt, 1e-2)print('distance ', sas.Value(uv).Distance(pt))return uv.Coord()class RadiusConstrainedSurface():'''returns a surface that has `radius` at `pt`'''def __init__(self, display, poly, pnt, targetRadius):self.display = displayself.targetRadius = targetRadiusself.poly = polyself.pnt = pntself.plate = self.build_surface()def build_surface(self):'''builds and renders the plate'''self.plate = build_plate([self.poly], [self.pnt])self.display.EraseAll()self.display.DisplayShape(self.plate)vert = make_vertex(self.pnt)self.display.DisplayShape(vert, update=True)def radius(self, z):'''sets the height of the point constraining the plate, returnsthe radius at this point'''if isinstance(z, float):self.pnt.SetX(z)else:self.pnt.SetX(float(z[0]))self.build_surface()uv = uv_from_projected_point_on_face(self.plate, self.pnt)radius = radius_at_uv(self.plate, uv[0], uv[1])print('z: ', z, 'radius: ', radius)self.curr_radius = radiusreturn self.targetRadius - abs(radius)def solve(self):fsolve(self.radius, 1, maxfev=1000)return self.platedef solve_radius(event=None):if not HAS_SCIPY:print("sorry cannot run solve_radius, scipy was not found...")returndisplay.EraseAll()p1 = gp_Pnt(0, 0, 0)p2 = gp_Pnt(0, 10, 0)p3 = gp_Pnt(0, 10, 10)p4 = gp_Pnt(0, 0, 10)p5 = gp_Pnt(5, 5, 5)poly = make_closed_polygon([p1, p2, p3, p4])for i in (0.1, 0.5, 1.5, 2., 3., 0.2):rcs = RadiusConstrainedSurface(display, poly, p5, i)rcs.solve()print('Goal: %s radius: %s' % (i, rcs.curr_radius))time.sleep(0.1)def build_geom_plate(edges):bpSrf = GeomPlate_BuildPlateSurface(3, 9, 12)# add curve constraintsfor edg in edges:c = BRepAdaptor_HCurve()print('edge:', edg)c.ChangeCurve().Initialize(edg)constraint = BRepFill_CurveConstraint(c, 0)bpSrf.Add(constraint)# add point constrainttry:bpSrf.Perform()except RuntimeError:print('failed to build the geom plate surface ')srf = bpSrf.Surface()plate = GeomPlate_MakeApprox(srf, 0.01, 10, 5, 0.01, 0, GeomAbs_C0)uMin, uMax, vMin, vMax = srf.Bounds()face = make_face(plate.Surface(), uMin, uMax, vMin, vMax, 1e-6)return facedef build_curve_network(event=None):'''mimic the curve network surfacing command from rhino'''print('Importing IGES file...')iges_file = os.path.join('..', 'assets', 'models', 'curve_geom_plate.igs')iges = read_iges_file(iges_file)print('Building geomplate...')topo = TopologyExplorer(iges)edges_list = list(topo.edges())face = build_geom_plate(edges_list)print('done.')display.EraseAll()display.DisplayShape(edges_list)display.DisplayShape(face)display.FitAll()print('Cutting out of edges...')def exit(event=None):sys.exit()if __name__ == "__main__":add_menu('geom plate')add_function_to_menu('geom plate', geom_plate)add_function_to_menu('geom plate', solve_radius)add_function_to_menu('geom plate', build_curve_network)add_function_to_menu('geom plate', exit)build_curve_network()start_display()
运行效果:生成几何板条