alignment api: spiral transition curves and cant for the PI method

Extends the PI method to horizontal alignments with clothoid spiral
transition curves and cant profiles for railway engineering, with the
geometry available as a pure, file-independent solver that any caller
(interactive editors, scripts, importers) can reuse.

solve_horizontal_alignment_by_pi_method() solves PI coordinates and radii
into a continuous sequence of HorizontalSegmentDefinition values whose
fields mirror IfcAlignmentHorizontalSegment. Each radii element is either
a circular curve radius R or a (R, Lin, Lout) sequence giving the entry
and exit spiral lengths. The circular curve is shifted inward so tangent
runs, spirals, and the circular curve remain continuous in position and
direction. Spiral end points are computed by Gauss-Legendre integration
of the clothoid position functions (compute_clothoid_end, exposed for
reuse), matching the exact clothoid to machine precision, and
compute_horizontal_segment_end() lets callers verify continuity of any
definition sequence without touching a file.

layout_horizontal_alignment_by_pi_method() is now a thin writer over the
solver: it accepts the same (R, Lin, Lout) elements and writes the
solved definitions to an IfcAlignmentHorizontal layout. When a cant
layout and cant values are given, cant segments are written one-for-one
with the horizontal segments, as expected by the cant segment lookup
used for Viennese bends: zero cant on tangent runs (CONSTANTCANT),
linearly varying cant over spirals (LINEARTRANSITION), and constant cant
over circular curves (CONSTANTCANT), applied to the rail on the outside
of the curve. Curves with a non-zero cant require both spirals so the
cant profile is continuous, which keeps the segmented reference curve
schema-valid. The parity is structural: every solved definition carries
its own cant values, so the two layouts cannot drift apart.

create_by_pi_method() passes radii elements through and accepts optional
cants and rail_head_distance; create() accepts rail_head_distance for
IfcAlignmentCant.RailHeadDistance. Plain-R input routes through the
verbatim legacy code path and produces unchanged output.

The tests check the quadrature against an independent dense Simpson
reference, pure-solver continuity (max position gap 1e-12), continuity
of the written geometry on the C++ evaluator (max position gap 8.5e-7,
the engine's clothoid series truncation), the cant profile and its
segment parity, error conditions, and an end-to-end
tangent-clothoid-arc-clothoid-tangent authoring example that validates
cleanly against the schema and express rules.

Part of #6890.

Generated with the assistance of an AI coding tool.
This commit is contained in:
Petru Conduraru
2026-07-21 17:49:40 +03:00
parent e52e5e2e58
commit 1a81292eb4
6 changed files with 794 additions and 108 deletions
@@ -32,7 +32,10 @@ This API does not determine alignment parameters based on rules, such as minimum
This API is under development and subject to code breaking changes in the future.
Presently, this API supports:
1. Creating alignments, both horizontal and vertical, using the PI method. Alignment definition can be read from a CSV file.
1. Creating alignments, both horizontal and vertical, using the PI method, including clothoid transition spirals and cant.
The horizontal PI solve is also available as a pure geometric computation (solve_horizontal_alignment_by_pi_method)
for callers that need segment parameters without writing to a file, such as interactive editors.
Alignment definition can be read from a CSV file.
2. Creating alignments segment by segment.
3. Automatic creation of geometric definitions (IfcCompositeCurve, IfcGradientCurve, IfcSegmentedReferenceCurve)
4. Automatic definition of stationing
@@ -40,11 +43,10 @@ Presently, this API supports:
6. Utility functions for printing business logical and geometric representations, as well as minimal geometry evaluations
Future versions of this API may support:
1. Defining alignments using the PI method, including transition spirals
2. Updating horizontal curve definitions by revising transition spiral parameters and circular curve radii
3. Updating vertical curve definitions by revising horizontal length of curves
4. Removing a segment at any location along a curve
5. Adding a segment at any location along a curve
1. Updating horizontal curve definitions by revising transition spiral parameters and circular curve radii
2. Updating vertical curve definitions by revising horizontal length of curves
3. Removing a segment at any location along a curve
4. Adding a segment at any location along a curve
"""
from ._get_segment_start_point_label import register_referent_name_callback
@@ -89,6 +91,12 @@ from .layout_vertical_alignment_by_pi_method import (
layout_vertical_alignment_by_pi_method,
)
from .name_segments import name_segments
from .solve_horizontal_alignment_by_pi_method import (
HorizontalSegmentDefinition,
compute_clothoid_end,
compute_horizontal_segment_end,
solve_horizontal_alignment_by_pi_method,
)
from .update_end_point import update_end_point
from .update_fallback_position import update_fallback_position
from .util import *
@@ -127,10 +135,14 @@ __all__ = [
"get_referent_nest",
"get_vertical_layout",
"has_zero_length_segment",
"HorizontalSegmentDefinition",
"compute_clothoid_end",
"compute_horizontal_segment_end",
"layout_horizontal_alignment_by_pi_method",
"layout_vertical_alignment_by_pi_method",
"name_segments",
"register_referent_name_callback",
"solve_horizontal_alignment_by_pi_method",
"update_end_point",
"update_fallback_position",
"get_mapped_segments",
@@ -35,6 +35,7 @@ def create(
include_cant: bool = False,
include_geometry: bool = True,
start_station: float = 0.0,
rail_head_distance: float = 1.0,
) -> entity_instance:
"""
Creates a new alignment with a horizontal layout. Optionally, vertical and cant layouts can be created as well.
@@ -58,6 +59,7 @@ def create(
:param include_cant: If True, IfcAlignmentCant is created. IfcSegmentedReferenceCurve is created if include_geometry is True
:param include_geometry: If True, the geometric representations are added
:param start_station: station value at the start of the alignment.
:param rail_head_distance: value assigned to IfcAlignmentCant.RailHeadDistance when include_cant is True
:return: Returns an IfcAlignment
"""
alignment = file.createIfcAlignment(
@@ -79,7 +81,9 @@ def create(
alignment_layouts.append(file.createIfcAlignmentVertical(GlobalId=ifcopenshell.guid.new()))
if include_cant:
alignment_layouts.append(file.createIfcAlignmentCant(GlobalId=ifcopenshell.guid.new(), RailHeadDistance=1.0))
alignment_layouts.append(
file.createIfcAlignmentCant(GlobalId=ifcopenshell.guid.new(), RailHeadDistance=rail_head_distance)
)
ifcopenshell.api.nest.assign_object(file, related_objects=alignment_layouts, relating_object=alignment)
@@ -17,6 +17,7 @@
# along with IfcOpenShell. If not, see <http://www.gnu.org/licenses/>.
from collections.abc import Sequence
from typing import Optional, Union
import ifcopenshell
import ifcopenshell.api.alignment
@@ -27,28 +28,52 @@ def create_by_pi_method(
file: ifcopenshell.file,
name: str,
hpoints: Sequence[Sequence[float]],
radii: Sequence[float],
radii: Sequence[Union[float, Sequence[float]]],
vpoints: Sequence[Sequence[float]] = None,
lengths: Sequence[float] = None,
start_station: float = 0.0,
cants: Optional[Sequence[float]] = None,
rail_head_distance: float = 1.0,
) -> entity_instance:
"""
Create an alignment using the PI layout method for both horizontal and vertical alignments.
If vpoints and lengths are omitted, only a horizontal alignment is created.
Each element of radii is either a circular curve radius R or a (R, Lin, Lout) sequence where
Lin and Lout are the lengths of clothoid spiral transition curves ahead of and following the
circular curve (see layout_horizontal_alignment_by_pi_method).
If cants is provided, a cant layout is created as well. Cant segments are created one-for-one
with the horizontal segments: zero cant on tangent runs, linearly varying cant over spiral
transitions, and constant cant over circular curves, applied to the rail on the outside of the
curve. A vertical alignment is required when cants is provided.
:param name: value for Name attribute
:param points: (X,Y) pairs denoting the location of the horizontal PIs, including start and end
:param radii: radii values to use for transition
:param radii: radii values to use for transition, optionally with spiral transition lengths
:param vpoints: (distance_along, Z_height) pairs denoting the location of the vertical PIs, including start and end.
:param lengths: parabolic vertical curve horizontal length values to use for transition
:param cants: cant values, one per PI curve, applied to the outer rail
:param rail_head_distance: value assigned to IfcAlignmentCant.RailHeadDistance
:return: Returns an IfcAlignment
"""
include_vertical = True if vpoints and lengths else False
include_cant = cants is not None
if include_cant and not include_vertical:
raise ValueError("a vertical alignment is required when cants is provided; supply vpoints and lengths")
alignment = ifcopenshell.api.alignment.create(
file, name, include_vertical=include_vertical, start_station=start_station
file,
name,
include_vertical=include_vertical,
include_cant=include_cant,
start_station=start_station,
rail_head_distance=rail_head_distance,
)
cant_layout = ifcopenshell.api.alignment.get_cant_layout(alignment) if include_cant else None
horizontal_layout = ifcopenshell.api.alignment.get_horizontal_layout(alignment)
ifcopenshell.api.alignment.layout_horizontal_alignment_by_pi_method(file, horizontal_layout, hpoints, radii)
ifcopenshell.api.alignment.layout_horizontal_alignment_by_pi_method(
file, horizontal_layout, hpoints, radii, cant_layout=cant_layout, cants=cants
)
if include_vertical:
vertical_layout = ifcopenshell.api.alignment.get_vertical_layout(alignment)
ifcopenshell.api.alignment.layout_vertical_alignment_by_pi_method(file, vertical_layout, vpoints, lengths)
@@ -16,128 +16,103 @@
# You should have received a copy of the GNU Lesser General Public License
# along with IfcOpenShell. If not, see <http://www.gnu.org/licenses/>.
import math
from collections.abc import Sequence
from typing import Optional, Union
import ifcopenshell
import ifcopenshell.api.alignment
import ifcopenshell.util.unit
from ifcopenshell import entity_instance
from ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method import (
HorizontalSegmentDefinition,
solve_horizontal_alignment_by_pi_method,
)
def _create_cant_segment(
file: ifcopenshell.file,
cant_layout: entity_instance,
segment: HorizontalSegmentDefinition,
) -> None:
"""
Appends the cant segment corresponding to one horizontal segment definition. Cant is applied
to a single rail (the rail on the outside of the curve). Constant cant is modeled with
CONSTANTCANT and varying cant (over a transition curve) with LINEARTRANSITION.
"""
is_transition = segment.start_cant != segment.end_cant
if segment.raise_left_rail:
start_left, start_right = segment.start_cant, 0.0
end_left, end_right = segment.end_cant, 0.0
else:
start_left, start_right = 0.0, segment.start_cant
end_left, end_right = 0.0, segment.end_cant
design_parameters = file.createIfcAlignmentCantSegment(
StartTag=None,
EndTag=None,
StartDistAlong=segment.start_dist_along,
HorizontalLength=segment.segment_length,
StartCantLeft=start_left,
EndCantLeft=end_left if is_transition else None,
StartCantRight=start_right,
EndCantRight=end_right if is_transition else None,
PredefinedType="LINEARTRANSITION" if is_transition else "CONSTANTCANT",
)
ifcopenshell.api.alignment.create_layout_segment(file, cant_layout, design_parameters)
def layout_horizontal_alignment_by_pi_method(
file: ifcopenshell.file, layout: entity_instance, hpoints: Sequence[Sequence[float]], radii: Sequence[float]
file: ifcopenshell.file,
layout: entity_instance,
hpoints: Sequence[Sequence[float]],
radii: Sequence[Union[float, Sequence[float]]],
cant_layout: Optional[entity_instance] = None,
cants: Optional[Sequence[float]] = None,
) -> None:
"""
Appends IfcAlignmentHorizontalSegment to a previously defined IfcAlignmentHorizontal using the PI layout method.
The zero length segment is updated.
The geometry is computed by solve_horizontal_alignment_by_pi_method; see that function for the
meaning of hpoints, radii, and cants. This function writes the resulting segment definitions to
the layout.
Optionally, a cant profile can be created alongside the horizontal layout. Cant segments are
created one-for-one with the horizontal segments: zero cant on tangent runs (CONSTANTCANT),
linearly varying cant over spiral transitions (LINEARTRANSITION), and constant cant over
circular curves (CONSTANTCANT). The cant is applied to the rail on the outside of the curve.
Cant values are expressed in the project length unit. IfcAlignmentCant.RailHeadDistance is
taken from the cant_layout. Curves with a non-zero cant require entry and exit spiral
transition curves so the cant profile is continuous.
:param file: file
:param layout: An IfcAlignmentHorizontal layout
:param hpoints: (X, Y) pairs denoting the location of the horizontal PIs, including start (POB) and end (POE).
:param radii: radius values to use for transition
:param radii: radius values to use for transition, optionally with spiral transition lengths
:param cant_layout: An IfcAlignmentCant layout to receive the cant segments. Required when cants is provided.
:param cants: cant values, one per PI curve, applied to the outer rail. Required when cant_layout is provided.
:return: None
"""
if not (len(hpoints) - 2 == len(radii)):
raise ValueError("radii should have two fewer elements that hpoints")
if (cant_layout is None) != (cants is None):
raise ValueError("cant_layout and cants must be provided together")
angle_unit_scale = ifcopenshell.util.unit.calculate_unit_scale(file, "PLANEANGLEUNIT")
xBT, yBT = hpoints[0]
xPI, yPI = hpoints[1]
i = 1
for radius in radii:
# back tangent
dxBT = xPI - xBT
dyBT = yPI - yBT
angleBT = math.atan2(dyBT, dxBT)
lengthBT = math.sqrt(dxBT * dxBT + dyBT * dyBT)
# forward tangent
i += 1
xFT, yFT = hpoints[i]
dxFT = xFT - xPI
dyFT = yFT - yPI
angleFT = math.atan2(dyFT, dxFT)
delta = angleFT - angleBT
tangent = abs(radius * math.tan(delta / 2))
lc = abs(radius * delta)
radius *= delta / abs(delta)
xPC = xPI - tangent * math.cos(angleBT)
yPC = yPI - tangent * math.sin(angleBT)
xPT = xPI + tangent * math.cos(angleFT)
yPT = yPI + tangent * math.sin(angleFT)
tangent_run = lengthBT - tangent
# create back tangent run
if 1.0e-03 < tangent_run:
pt = file.createIfcCartesianPoint(
Coordinates=(xBT, yBT),
)
design_parameters = file.createIfcAlignmentHorizontalSegment(
StartTag=None,
EndTag=None,
StartPoint=pt,
StartDirection=angleBT / angle_unit_scale,
StartRadiusOfCurvature=0.0,
EndRadiusOfCurvature=0.0,
SegmentLength=tangent_run,
GravityCenterLineHeight=None,
PredefinedType="LINE",
)
ifcopenshell.api.alignment.create_layout_segment(file, layout, design_parameters)
# create circular curve
if radius != 0.0:
pc = file.createIfcCartesianPoint(
Coordinates=(xPC, yPC),
)
design_parameters = file.createIfcAlignmentHorizontalSegment(
StartTag=None,
EndTag=None,
StartPoint=pc,
StartDirection=angleBT / angle_unit_scale,
StartRadiusOfCurvature=float(radius),
EndRadiusOfCurvature=float(radius),
SegmentLength=lc,
GravityCenterLineHeight=None,
PredefinedType="CIRCULARARC",
)
ifcopenshell.api.alignment.create_layout_segment(file, layout, design_parameters)
xBT = xPT
yBT = yPT
xPI = xFT
yPI = yFT
# done processing radii
# create last tangent run
dx = xPI - xBT
dy = yPI - yBT
angleBT = math.atan2(dy, dx)
tangent_run = math.sqrt(dx * dx + dy * dy)
if 1.0e-03 < tangent_run:
pt = file.createIfcCartesianPoint(Coordinates=(xBT, yBT))
for segment in solve_horizontal_alignment_by_pi_method(hpoints, radii, cants):
if cant_layout is not None:
_create_cant_segment(file, cant_layout, segment)
start_point = file.createIfcCartesianPoint(
Coordinates=segment.start_point,
)
design_parameters = file.createIfcAlignmentHorizontalSegment(
StartTag=None,
EndTag=None,
StartPoint=pt,
StartDirection=angleBT / angle_unit_scale,
StartRadiusOfCurvature=0.0,
EndRadiusOfCurvature=0.0,
SegmentLength=tangent_run,
StartPoint=start_point,
StartDirection=segment.start_direction / angle_unit_scale,
StartRadiusOfCurvature=segment.start_radius_of_curvature,
EndRadiusOfCurvature=segment.end_radius_of_curvature,
SegmentLength=segment.segment_length,
GravityCenterLineHeight=None,
PredefinedType="LINE",
PredefinedType=segment.predefined_type,
)
ifcopenshell.api.alignment.create_layout_segment(file, layout, design_parameters)
@@ -0,0 +1,423 @@
# IfcOpenShell - IFC toolkit and geometry engine
# Copyright (C) 2025 Thomas Krijnen <thomas@aecgeeks.com>
#
# This file is part of IfcOpenShell.
#
# IfcOpenShell 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.
#
# IfcOpenShell 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 IfcOpenShell. If not, see <http://www.gnu.org/licenses/>.
# This file was generated with the assistance of an AI coding tool.
import math
from collections.abc import Sequence
from typing import NamedTuple, Optional, Union
import numpy as np
# Gauss-Legendre quadrature nodes and weights used to integrate the clothoid position functions
_gauss_legendre_points = np.polynomial.legendre.leggauss(32)
class HorizontalSegmentDefinition(NamedTuple):
"""
Parameters of one horizontal alignment segment produced by the PI method solver.
The fields mirror IfcAlignmentHorizontalSegment so a definition can be written to a file
without further computation, but the definition itself is independent of any file. Directions
are in radians, lengths and coordinates in the caller's length unit.
"""
start_point: tuple[float, float]
"""(X, Y) of the segment start"""
start_direction: float
"""direction of the tangent at the segment start, in radians"""
start_radius_of_curvature: float
"""radius at the segment start; 0.0 for straight, positive curving left, negative curving right"""
end_radius_of_curvature: float
"""radius at the segment end, with the same sign convention as start_radius_of_curvature"""
segment_length: float
"""length of the segment along the curve"""
predefined_type: str
"""IfcAlignmentHorizontalSegmentTypeEnum value: LINE, CLOTHOID, or CIRCULARARC"""
start_dist_along: float = 0.0
"""distance along the alignment at the segment start"""
start_cant: float = 0.0
"""cant at the segment start, applied to the rail on the outside of the curve"""
end_cant: float = 0.0
"""cant at the segment end, applied to the rail on the outside of the curve"""
raise_left_rail: bool = False
"""True when the outside of the curve is the left rail (a curve to the right)"""
def compute_clothoid_end(length: float, start_curvature: float, end_curvature: float) -> tuple[float, float, float]:
"""
Computes the end point of a clothoid transition whose curvature varies linearly from
start_curvature to end_curvature over length.
The result (dx, dy, dtheta) is relative to the start of the transition, with the x-axis in the
direction of the tangent at the start. Curvatures are signed: positive curving left, negative
curving right. The position is computed with 32 point Gauss-Legendre quadrature of the clothoid
integrals.
:param length: length of the transition, measured along the curve
:param start_curvature: curvature at the start (1/R, 0.0 for a straight)
:param end_curvature: curvature at the end (1/R, 0.0 for a straight)
:return: (dx, dy, dtheta) displacement and change in tangent direction over the transition
"""
u, w = _gauss_legendre_points
l = 0.5 * length * (u + 1.0) # map quadrature points from (-1,1) onto (0,length)
theta = start_curvature * l + (end_curvature - start_curvature) * l * l / (2.0 * length)
dx = 0.5 * length * float(np.sum(w * np.cos(theta)))
dy = 0.5 * length * float(np.sum(w * np.sin(theta)))
dtheta = 0.5 * (start_curvature + end_curvature) * length
return dx, dy, dtheta
def compute_horizontal_segment_end(segment: HorizontalSegmentDefinition) -> tuple[float, float, float]:
"""
Computes the end point and end direction of a horizontal segment definition.
Useful for checking position and direction continuity between consecutive segments: the result
for one segment should match the start_point and start_direction of the next.
:param segment: the segment definition
:return: (x, y, direction) at the end of the segment, direction in radians
"""
x, y = segment.start_point
direction = segment.start_direction
length = segment.segment_length
if segment.predefined_type == "LINE":
return (x + length * math.cos(direction), y + length * math.sin(direction), direction)
start_curvature = 1.0 / segment.start_radius_of_curvature if segment.start_radius_of_curvature != 0.0 else 0.0
end_curvature = 1.0 / segment.end_radius_of_curvature if segment.end_radius_of_curvature != 0.0 else 0.0
if segment.predefined_type == "CIRCULARARC":
dtheta = start_curvature * length
dx = math.sin(dtheta) / start_curvature
dy = (1.0 - math.cos(dtheta)) / start_curvature
elif segment.predefined_type == "CLOTHOID":
dx, dy, dtheta = compute_clothoid_end(length, start_curvature, end_curvature)
else:
raise NotImplementedError(f"unsupported predefined type '{segment.predefined_type}'")
return (
x + dx * math.cos(direction) - dy * math.sin(direction),
y + dx * math.sin(direction) + dy * math.cos(direction),
direction + dtheta,
)
def solve_horizontal_alignment_by_pi_method(
hpoints: Sequence[Sequence[float]],
radii: Sequence[Union[float, Sequence[float]]],
cants: Optional[Sequence[float]] = None,
) -> list[HorizontalSegmentDefinition]:
"""
Solves a horizontal alignment defined by the PI layout method into a continuous sequence of
segment definitions.
This is a pure geometric computation: no file is read or written. Use
layout_horizontal_alignment_by_pi_method to write the solution to an IfcAlignmentHorizontal
layout, or consume the returned definitions directly, for example to preview an alignment in
an interactive editor before committing it to a file.
Each element of radii defines the transition at the corresponding PI and is either:
R - radius of a circular curve (tangent runs connect directly to the circular curve), or
(R, Lin, Lout) - radius of a circular curve with clothoid spiral transition curves of length
Lin ahead of the curve and Lout following the curve. When spiral transitions are used the
circular curve shifts inward relative to the tangent runs so the tangent runs, spirals, and
circular curve are continuous in position and direction. Lin and Lout can be 0.0 for a
spiral-less connection on that end of the curve.
If cants is provided, each definition also carries the cant at the segment start and end,
applied to the rail on the outside of the curve: zero cant on tangent runs, linearly varying
cant over spiral transitions, and constant cant over circular curves. Because every horizontal
segment carries its own cant values, a cant layout built from the definitions is one-for-one
with the horizontal layout. Curves with a non-zero cant require entry and exit spiral
transition curves so the cant profile is continuous.
:param hpoints: (X, Y) pairs denoting the location of the horizontal PIs, including start (POB) and end (POE).
:param radii: radius values to use for transition, optionally with spiral transition lengths
:param cants: cant values, one per PI curve, applied to the outer rail
:return: list of segment definitions, in order, continuous in position and direction
"""
if not (len(hpoints) - 2 == len(radii)):
raise ValueError("radii should have two fewer elements that hpoints")
if cants is not None and not (len(cants) == len(radii)):
raise ValueError("cants should have the same number of elements as radii")
segments: list[HorizontalSegmentDefinition] = []
xBT, yBT = hpoints[0]
xPI, yPI = hpoints[1]
i = 1
dist_along = 0.0 # distance along the horizontal alignment at the start of the next segment
for curve_index, curve in enumerate(radii):
if isinstance(curve, (int, float)):
radius = float(curve)
entry_length = 0.0
exit_length = 0.0
else:
if len(curve) != 3:
raise ValueError("each radii element should be a radius R or a (R, Lin, Lout) sequence")
radius, entry_length, exit_length = (float(v) for v in curve)
if radius == 0.0 and (entry_length != 0.0 or exit_length != 0.0):
raise ValueError("spiral transition lengths require a non-zero radius")
cant = float(cants[curve_index]) if cants is not None else 0.0
if cant != 0.0 and (entry_length == 0.0 or exit_length == 0.0):
raise ValueError(
"curves with a non-zero cant require entry and exit spiral transition curves; "
"otherwise the cant profile is discontinuous"
)
# back tangent
dxBT = xPI - xBT
dyBT = yPI - yBT
angleBT = math.atan2(dyBT, dxBT)
lengthBT = math.sqrt(dxBT * dxBT + dyBT * dyBT)
# forward tangent
i += 1
xFT, yFT = hpoints[i]
dxFT = xFT - xPI
dyFT = yFT - yPI
angleFT = math.atan2(dyFT, dxFT)
delta = angleFT - angleBT
if entry_length == 0.0 and exit_length == 0.0:
# tangent runs connect directly to the circular curve
tangent = abs(radius * math.tan(delta / 2))
lc = abs(radius * delta)
radius *= delta / abs(delta)
xPC = xPI - tangent * math.cos(angleBT)
yPC = yPI - tangent * math.sin(angleBT)
xPT = xPI + tangent * math.cos(angleFT)
yPT = yPI + tangent * math.sin(angleFT)
tangent_run = lengthBT - tangent
# back tangent run
if 1.0e-03 < tangent_run:
segments.append(
HorizontalSegmentDefinition(
start_point=(xBT, yBT),
start_direction=angleBT,
start_radius_of_curvature=0.0,
end_radius_of_curvature=0.0,
segment_length=tangent_run,
predefined_type="LINE",
start_dist_along=dist_along,
raise_left_rail=delta < 0.0,
)
)
dist_along += tangent_run
# circular curve
if radius != 0.0:
segments.append(
HorizontalSegmentDefinition(
start_point=(xPC, yPC),
start_direction=angleBT,
start_radius_of_curvature=float(radius),
end_radius_of_curvature=float(radius),
segment_length=lc,
predefined_type="CIRCULARARC",
start_dist_along=dist_along,
start_cant=cant,
end_cant=cant,
raise_left_rail=delta < 0.0,
)
)
dist_along += lc
else:
# tangent runs connect to the circular curve with clothoid spiral transition curves.
# normalize the deflection angle onto (-pi, pi)
delta = math.atan2(math.sin(delta), math.cos(delta))
if delta == 0.0:
raise ValueError("PI deflection angle is zero; spiral transitions cannot be created")
R = abs(radius)
s = 1.0 if 0.0 < delta else -1.0 # +1 curve to the left, -1 curve to the right
theta1 = entry_length / (2.0 * R) # deflection of the entry spiral
theta2 = exit_length / (2.0 * R) # deflection of the exit spiral
theta_c = abs(delta) - theta1 - theta2 # deflection of the circular curve
if theta_c < 0.0:
raise ValueError(
"spiral transition curves are too long; their combined deflection exceeds the PI deflection angle"
)
lc = R * theta_c
# compose the displacement from the start of the entry spiral (TS) to the end of the
# exit spiral (ST), in a frame with the x-axis along the back tangent.
# pieces are computed for a curve to the left and mirrored by s.
pieces = []
if 0.0 < entry_length:
pieces.append(compute_clothoid_end(entry_length, 0.0, 1.0 / R))
pieces.append((R * math.sin(theta_c), R * (1.0 - math.cos(theta_c)), theta_c))
if 0.0 < exit_length:
pieces.append(compute_clothoid_end(exit_length, 1.0 / R, 0.0))
x = 0.0
y = 0.0
direction = 0.0
for dx_, dy_, dtheta_ in pieces:
x += dx_ * math.cos(direction) - s * dy_ * math.sin(direction)
y += dx_ * math.sin(direction) + s * dy_ * math.cos(direction)
direction += s * dtheta_
# locate TS on the back tangent and ST on the forward tangent so that the curve ends on
# the forward tangent. this accounts for the inward shift of the circular curve.
ts_to_pi = x - y / math.tan(delta) # distance from TS to the PI, along the back tangent
pi_to_st = y / math.sin(delta) # distance from the PI to ST, along the forward tangent
tangent_run = lengthBT - ts_to_pi
# back tangent run
if 1.0e-03 < tangent_run:
segments.append(
HorizontalSegmentDefinition(
start_point=(xBT, yBT),
start_direction=angleBT,
start_radius_of_curvature=0.0,
end_radius_of_curvature=0.0,
segment_length=tangent_run,
predefined_type="LINE",
start_dist_along=dist_along,
raise_left_rail=delta < 0.0,
)
)
dist_along += tangent_run
signed_radius = s * R
cur_x = xPI - ts_to_pi * math.cos(angleBT)
cur_y = yPI - ts_to_pi * math.sin(angleBT)
cur_direction = angleBT
# entry spiral
if 0.0 < entry_length:
segments.append(
HorizontalSegmentDefinition(
start_point=(cur_x, cur_y),
start_direction=cur_direction,
start_radius_of_curvature=0.0,
end_radius_of_curvature=signed_radius,
segment_length=entry_length,
predefined_type="CLOTHOID",
start_dist_along=dist_along,
start_cant=0.0,
end_cant=cant,
raise_left_rail=delta < 0.0,
)
)
dist_along += entry_length
dx_, dy_, dtheta_ = compute_clothoid_end(entry_length, 0.0, 1.0 / R)
cur_x += dx_ * math.cos(cur_direction) - s * dy_ * math.sin(cur_direction)
cur_y += dx_ * math.sin(cur_direction) + s * dy_ * math.cos(cur_direction)
cur_direction += s * dtheta_
# circular curve
if 1.0e-03 < lc:
segments.append(
HorizontalSegmentDefinition(
start_point=(cur_x, cur_y),
start_direction=cur_direction,
start_radius_of_curvature=signed_radius,
end_radius_of_curvature=signed_radius,
segment_length=lc,
predefined_type="CIRCULARARC",
start_dist_along=dist_along,
start_cant=cant,
end_cant=cant,
raise_left_rail=delta < 0.0,
)
)
dist_along += lc
cur_x += R * math.sin(theta_c) * math.cos(cur_direction) - s * R * (1.0 - math.cos(theta_c)) * math.sin(
cur_direction
)
cur_y += R * math.sin(theta_c) * math.sin(cur_direction) + s * R * (1.0 - math.cos(theta_c)) * math.cos(
cur_direction
)
cur_direction += s * theta_c
# exit spiral
if 0.0 < exit_length:
segments.append(
HorizontalSegmentDefinition(
start_point=(cur_x, cur_y),
start_direction=cur_direction,
start_radius_of_curvature=signed_radius,
end_radius_of_curvature=0.0,
segment_length=exit_length,
predefined_type="CLOTHOID",
start_dist_along=dist_along,
start_cant=cant,
end_cant=0.0,
raise_left_rail=delta < 0.0,
)
)
dist_along += exit_length
xPT = xPI + pi_to_st * math.cos(angleFT)
yPT = yPI + pi_to_st * math.sin(angleFT)
xBT = xPT
yBT = yPT
xPI = xFT
yPI = yFT
# done processing radii
# last tangent run
dx = xPI - xBT
dy = yPI - yBT
angleBT = math.atan2(dy, dx)
tangent_run = math.sqrt(dx * dx + dy * dy)
if 1.0e-03 < tangent_run:
segments.append(
HorizontalSegmentDefinition(
start_point=(xBT, yBT),
start_direction=angleBT,
start_radius_of_curvature=0.0,
end_radius_of_curvature=0.0,
segment_length=tangent_run,
predefined_type="LINE",
start_dist_along=dist_along,
)
)
return segments
@@ -0,0 +1,247 @@
# IfcOpenShell - IFC toolkit and geometry engine
# Copyright (C) 2025 Thomas Krijnen <thomas@aecgeeks.com>
#
# This file is part of IfcOpenShell.
#
# IfcOpenShell 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.
#
# IfcOpenShell 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 IfcOpenShell. If not, see <http://www.gnu.org/licenses/>.
# This file was generated with the assistance of an AI coding tool.
import logging
import math
import numpy as np
import pytest
import ifcopenshell.api.alignment
import ifcopenshell.api.context
import ifcopenshell.api.unit
import ifcopenshell.geom
import ifcopenshell.validate
from ifcopenshell import ifcopenshell_wrapper
def _create_file():
file = ifcopenshell.file(schema="IFC4X3")
project = file.createIfcProject(GlobalId=ifcopenshell.guid.new(), Name="Test")
length = ifcopenshell.api.unit.add_si_unit(file, unit_type="LENGTHUNIT")
ifcopenshell.api.unit.assign_unit(file, units=[length])
geometric_representation_context = ifcopenshell.api.context.add_context(file, context_type="Model")
axis_model_representation_subcontext = ifcopenshell.api.context.add_context(
file,
context_type="Model",
context_identifier="Axis",
target_view="MODEL_VIEW",
parent=geometric_representation_context,
)
return file
def _reference_clothoid_end(length, start_curvature, end_curvature, steps=20000):
"""Composite Simpson integration of the clothoid position functions, as an independent check."""
l = np.linspace(0.0, length, 2 * steps + 1)
theta = start_curvature * l + (end_curvature - start_curvature) * l * l / (2.0 * length)
h = length / (2.0 * steps)
weights = np.ones(2 * steps + 1)
weights[1:-1:2] = 4.0
weights[2:-1:2] = 2.0
dx = h / 3.0 * float(np.sum(weights * np.cos(theta)))
dy = h / 3.0 * float(np.sum(weights * np.sin(theta)))
return dx, dy
def test_compute_clothoid_end():
for length, k1, k2 in [(200.0, 0.0, 1.0 / 1000.0), (150.0, 1.0 / 1000.0, 0.0), (120.0, -1.0 / 800.0, 1.0 / 500.0)]:
dx, dy, dtheta = ifcopenshell.api.alignment.compute_clothoid_end(length, k1, k2)
ref_dx, ref_dy = _reference_clothoid_end(length, k1, k2)
assert dx == pytest.approx(ref_dx, abs=1.0e-12)
assert dy == pytest.approx(ref_dy, abs=1.0e-12)
assert dtheta == pytest.approx(0.5 * (k1 + k2) * length)
# signed curvatures mirror the unsigned result
dx, dy, dtheta = ifcopenshell.api.alignment.compute_clothoid_end(200.0, 0.0, 1.0 / 1000.0)
mx, my, mtheta = ifcopenshell.api.alignment.compute_clothoid_end(200.0, 0.0, -1.0 / 1000.0)
assert mx == pytest.approx(dx)
assert my == pytest.approx(-dy)
assert mtheta == pytest.approx(-dtheta)
def test_solve_produces_continuous_segments():
hpoints = [(500.0, 2500.0), (3340.0, 660.0), (4340.0, 5000.0), (7600.0, 4560.0), (8480.0, 2010.0)]
radii = [(1000.0, 200.0, 150.0), (1250.0, 180.0, 180.0), (950.0, 0.0, 120.0)]
segments = ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, radii)
expected_types = [
"LINE",
"CLOTHOID",
"CIRCULARARC",
"CLOTHOID",
"LINE",
"CLOTHOID",
"CIRCULARARC",
"CLOTHOID",
"LINE",
"CIRCULARARC",
"CLOTHOID",
"LINE",
]
assert [s.predefined_type for s in segments] == expected_types
# the solution starts at the POB, in the direction of the first PI
assert segments[0].start_point == pytest.approx((500.0, 2500.0))
assert segments[0].start_direction == pytest.approx(math.atan2(660.0 - 2500.0, 3340.0 - 500.0))
# spirals run from zero curvature to the curve radius and vice versa
entry_spiral = segments[1]
assert entry_spiral.start_radius_of_curvature == 0.0
assert entry_spiral.end_radius_of_curvature == pytest.approx(1000.0) # positive, curve to the left
assert entry_spiral.segment_length == pytest.approx(200.0)
exit_spiral = segments[3]
assert exit_spiral.start_radius_of_curvature == pytest.approx(1000.0)
assert exit_spiral.end_radius_of_curvature == 0.0
assert exit_spiral.segment_length == pytest.approx(150.0)
assert segments[5].end_radius_of_curvature == pytest.approx(-1250.0) # curve to the right
# each segment ends exactly where the next one starts, in position and direction
dist_along = 0.0
for segment, next_segment in zip(segments[:-1], segments[1:]):
assert segment.start_dist_along == pytest.approx(dist_along)
end_x, end_y, end_direction = ifcopenshell.api.alignment.compute_horizontal_segment_end(segment)
assert end_x == pytest.approx(next_segment.start_point[0], abs=1.0e-9)
assert end_y == pytest.approx(next_segment.start_point[1], abs=1.0e-9)
direction_gap = end_direction - next_segment.start_direction
assert math.atan2(math.sin(direction_gap), math.cos(direction_gap)) == pytest.approx(0.0, abs=1.0e-12)
dist_along += segment.segment_length
def test_solve_plain_radius_matches_spiral_free_tuple():
hpoints = [(500.0, 2500.0), (3340.0, 660.0), (4340.0, 5000.0), (8480.0, 2010.0)]
segments1 = ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, [1000.0, 1250.0])
segments2 = ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(
hpoints, [(1000.0, 0.0, 0.0), (1250.0, 0.0, 0.0)]
)
assert segments1 == segments2
def test_solve_cant_profile():
hpoints = [(500.0, 2500.0), (3340.0, 660.0), (4340.0, 5000.0), (7600.0, 4560.0), (8480.0, 2010.0)]
radii = [(1000.0, 200.0, 150.0), (1250.0, 180.0, 180.0), (950.0, 140.0, 120.0)]
cants = [0.15, 0.12, 0.1]
segments = ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, radii, cants)
# cant varies linearly over spirals, is constant over circular curves, and is zero on tangents
for segment in segments:
if segment.predefined_type == "LINE":
assert segment.start_cant == 0.0 and segment.end_cant == 0.0
elif segment.predefined_type == "CIRCULARARC":
assert segment.start_cant == segment.end_cant != 0.0
entry_spiral = segments[1]
assert (entry_spiral.start_cant, entry_spiral.end_cant) == (0.0, 0.15)
exit_spiral = segments[3]
assert (exit_spiral.start_cant, exit_spiral.end_cant) == (0.15, 0.0)
# the first curve is to the left, so the outer rail is the right rail
assert entry_spiral.raise_left_rail is False
# the second curve is to the right, so the outer rail is the left rail
assert segments[5].raise_left_rail is True
# the cant profile is continuous across segment boundaries
for segment, next_segment in zip(segments[:-1], segments[1:]):
assert segment.end_cant == pytest.approx(next_segment.start_cant)
def test_solve_errors():
hpoints = [(0.0, 0.0), (1000.0, 0.0), (2000.0, 1000.0)]
with pytest.raises(ValueError): # radii count mismatch
ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, [500.0, 500.0])
with pytest.raises(ValueError): # cants count mismatch
ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, [(500.0, 50.0, 50.0)], [0.1, 0.1])
with pytest.raises(ValueError): # malformed radii element
ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, [(500.0, 50.0)])
with pytest.raises(ValueError): # spiral lengths without a radius
ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, [(0.0, 50.0, 50.0)])
with pytest.raises(ValueError): # cant without spiral transitions is discontinuous
ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, [(500.0, 0.0, 0.0)], [0.1])
with pytest.raises(ValueError): # spirals deflect more than the PI deflection angle
ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(hpoints, [(500.0, 5000.0, 5000.0)])
with pytest.raises(ValueError): # zero deflection angle
ifcopenshell.api.alignment.solve_horizontal_alignment_by_pi_method(
[(0.0, 0.0), (1000.0, 0.0), (2000.0, 0.0)], [(500.0, 50.0, 50.0)]
)
def test_author_transition_curve_alignment():
"""
End-to-end example: author a tangent -> clothoid -> circular arc -> clothoid -> tangent
alignment with cant, then check the written geometry for continuity with the geometry engine
and validate the file against the schema and express rules.
"""
file = _create_file()
alignment = ifcopenshell.api.alignment.create_by_pi_method(
file,
"TestAlignment",
[(500.0, 2500.0), (3340.0, 660.0), (4340.0, 5000.0)],
[(1000.0, 200.0, 150.0)],
[(0.0, 100.0), (2000.0, 135.0), (4000.0, 105.0)],
[1600.0],
cants=[0.15],
rail_head_distance=1.5,
)
horizontal_layout = ifcopenshell.api.alignment.get_horizontal_layout(alignment)
segment_nest = ifcopenshell.api.alignment.get_alignment_segment_nest(horizontal_layout)
expected_types = ["LINE", "CLOTHOID", "CIRCULARARC", "CLOTHOID", "LINE", "LINE"] # last is the zero length segment
assert [s.DesignParameters.PredefinedType for s in segment_nest.RelatedObjects] == expected_types
cant_layout = ifcopenshell.api.alignment.get_cant_layout(alignment)
assert cant_layout.RailHeadDistance == pytest.approx(1.5)
cant_nest = ifcopenshell.api.alignment.get_alignment_segment_nest(cant_layout)
# cant segments correspond one-for-one with the horizontal segments
assert len(cant_nest.RelatedObjects) == len(segment_nest.RelatedObjects)
# verify continuity of position and direction between consecutive segments of the
# geometric representation
curve = ifcopenshell.api.alignment.get_layout_curve(horizontal_layout)
settings = ifcopenshell.geom.settings()
for segment, next_segment in zip(curve.Segments[:-1], curve.Segments[1:]):
fn = ifcopenshell_wrapper.map_shape(settings, segment.wrapped_data)
evaluator = ifcopenshell_wrapper.function_item_evaluator(settings, fn)
end = np.array(evaluator.evaluate(fn.end()))
end_position = end[0:2, 3]
end_direction = math.atan2(end[1, 0], end[0, 0])
start_position = next_segment.Placement.Location.Coordinates
d = next_segment.Placement.RefDirection.DirectionRatios
start_direction = math.atan2(d[1], d[0])
assert end_position[0] == pytest.approx(start_position[0], abs=1.0e-5)
assert end_position[1] == pytest.approx(start_position[1], abs=1.0e-5)
direction_gap = math.atan2(math.sin(end_direction - start_direction), math.cos(end_direction - start_direction))
assert direction_gap == pytest.approx(0.0, abs=1.0e-9)
# the file is schema and express rule valid
logger = ifcopenshell.validate.json_logger()
ifcopenshell.validate.validate(file, logger, express_rules=True)
assert [entry for entry in logger.statements if entry["level"] == logging.ERROR] == []
test_compute_clothoid_end()
test_solve_produces_continuous_segments()
test_solve_plain_radius_matches_spiral_free_tuple()
test_solve_cant_profile()
test_solve_errors()
test_author_transition_curve_alignment()