17046961068

Committed 18 Aug 2025 04:51PM UTC coverage: 17.141% (+1.5%) from 15.675%

Build # 17046961068

Build Type

push

github

Committed by

web-flow

Commit Message

update to v0.18.4

Run Details

277 of 2402 new or added lines in 25 files covered. (11.53%)

13 existing lines in 6 files now uncovered.

3231 of 18850 relevant lines covered (17.14%)

0.51 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

14.43

/src/sdynpy/signal_processing/sdynpy_geometry_fitting.py

# -*- coding: utf-8 -*-
"""
Functions to fit geometry to point data

This module defines a Geometry object as well as all of the subcomponents of
a geometry object: nodes, elements, tracelines and coordinate system.  Geometry
plotting is also handled in this module.

Copyright 2022 National Technology & Engineering Solutions of Sandia,
LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
Government retains certain rights in this software.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""


import numpy as np
import scipy.optimize as opt


def cone_fit(points, origin, direction, angle):
    a = origin
    p = points
    n = direction
    height = np.einsum('...i,i->...', (a - p), n)
    distance = np.linalg.norm((a-p) - (height[:, np.newaxis])*n, axis=1)
    expected_distance = -np.tan(angle)*height
    return np.linalg.norm((distance-expected_distance))


def create_cone(ncircum, naxial, height, angle, origin, direction):
    circum_angles = np.linspace(0, 2*np.pi, ncircum, endpoint=False)
    heights = np.linspace(0, height, naxial+1)[1:]
    points = []
    for i, cangle in enumerate(circum_angles):
        for j, h in enumerate(heights):
            radius = np.tan(angle)*h
            points.append([np.sin(cangle)*radius, np.cos(cangle)*radius, h])
    points = np.array(points)
    rotation_matrix_z = direction/np.linalg.norm(direction)
    rotation_matrix_y = np.cross(rotation_matrix_z, [1.0, 0.0, 0.0])
    if np.linalg.norm(rotation_matrix_y) < 1e-10:
        rotation_matrix_x = np.cross([0.0, 1.0, 0.0], rotation_matrix_z)
        rotation_matrix_x /= np.linalg.norm(rotation_matrix_x)
        rotation_matrix_y = np.cross(rotation_matrix_z, rotation_matrix_x)
        rotation_matrix_y /= np.linalg.norm(rotation_matrix_y)
    else:
        rotation_matrix_y /= np.linalg.norm(rotation_matrix_y)
        rotation_matrix_x = np.cross(rotation_matrix_y, rotation_matrix_z)
        rotation_matrix_x /= np.linalg.norm(rotation_matrix_x)
    R = np.array((rotation_matrix_x, rotation_matrix_y, rotation_matrix_z)).T
    new_points = ((R@points.T).T+origin)
    return new_points


def cone_error_fn_fixed_angle(points, angle):
    return lambda x: cone_fit(points, x[:3], x[3:], angle)


def cone_error_fn_free_angle(points):
    return lambda x: cone_fit(points, x[:3], x[3:6], x[6])


def fit_cone_fixed_angle(points, angle):

    def constraint_eqn(x):
        return np.linalg.norm(x[3:])

    opt_fn = cone_error_fn_fixed_angle(points, angle)
    origin_guess = np.mean(points, axis=0)
    u, s, vh = np.linalg.svd(points-origin_guess)
    direction_guess = vh[0]
    x0_1 = np.concatenate((origin_guess, direction_guess))
    x0_2 = np.concatenate((origin_guess, -direction_guess))
    # Constrain to a unit vector
    constraint = opt.NonlinearConstraint(constraint_eqn, 1.0, 1.0)
    result_1 = opt.minimize(opt_fn, x0_1, method='trust-constr', constraints=constraint)
    result_2 = opt.minimize(opt_fn, x0_2, method='trust-constr', constraints=constraint)
    # Figure out which is the better one
    result = result_1 if result_1.fun < result_2.fun else result_2
    cone_origin = result.x[:3]
    cone_direction = result.x[3:]
    return result.fun, cone_origin, cone_direction


def distance_point_line(points, origin, direction):
    return np.linalg.norm(np.cross((points-origin), direction), axis=-1)/np.linalg.norm(direction)


def distance_point_plane(point, plane_point, plane_direction):
    return abs(np.einsum('...i,...i->...', point-plane_point, plane_direction))/np.linalg.norm(plane_direction, axis=-1)


def cylinder_fit(points, origin, direction, radius):
    return np.linalg.norm(distance_point_line(points, origin, direction) - radius)


def fit_cylinder(points, origin_guess=None, direction_guess=None, radius_guess=None):
    def opt_fn(x):
        return cylinder_fit(points, x[:3], x[3:6], x[6:7])

    def constraint_eqn(x):
        return np.linalg.norm(x[3:6])

    if origin_guess is None:
        origin_guess = np.mean(points, axis=0)
    if direction_guess is None or direction_guess == 'largest':
        u, s, vh = np.linalg.svd(points-origin_guess)
        direction_guess = vh[0]
    elif direction_guess == 'smallest':
        u, s, vh = np.linalg.svd(points-origin_guess)
        direction_guess = vh[-1]
    if radius_guess is None:
        u, s, vh = np.linalg.svd(points-origin_guess)
        radius_guess = s[0]/2
    x0 = np.concatenate((origin_guess, direction_guess, [radius_guess]))
    # Constrain to a unit vector
    constraint = opt.NonlinearConstraint(constraint_eqn, 1.0, 1.0)
    result = opt.minimize(opt_fn, x0, method='trust-constr', constraints=constraint)
    origin = result.x[:3]
    direction = result.x[3:6]
    radius = result.x[6]
    return result.fun, origin, direction, radius


def fit_line_point_cloud(points):
    center = np.mean(points, axis=0)
    shifted_points = points-center
    cov = points.T@shifted_points
    evals, evects = np.linalg.eigh(cov)
    direction = evects[:, np.argmax(evals)]
    direction /= np.linalg.norm(direction)
    return lambda t: center+t*direction

def intersection_point_multiple_lines(P0,P1):
    """Computes the intersection point of multiple lines in a least squares sense
    
    P0 and P1 are points on a line.  Each are NxD where N is the number of lines
    and D is the dimensionality of the space (D=2 corresponds to lines on a plane,
    D=3 is points in 3D space).
    """
    n = (P1-P0)/np.linalg.norm(P1-P0,axis=1)[:,np.newaxis]
    projectors = np.eye(n.shape[1]) - n[:,:,np.newaxis]*n[:,np.newaxis]  # I - n*n.T
    R = projectors.sum(axis=0)
    q = (projectors @ P0[:,:,np.newaxis]).sum(axis=0)
    p = np.linalg.lstsq(R,q,rcond=None)[0]

    return p

1	# -- coding: utf-8 --
2	"""	1✔
3	Functions to fit geometry to point data
4
5	This module defines a Geometry object as well as all of the subcomponents of
6	a geometry object: nodes, elements, tracelines and coordinate system. Geometry
7	plotting is also handled in this module.
8
9	Copyright 2022 National Technology & Engineering Solutions of Sandia,
10	LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
11	Government retains certain rights in this software.
12
13	This program is free software: you can redistribute it and/or modify
14	it under the terms of the GNU General Public License as published by
15	the Free Software Foundation, either version 3 of the License, or
16	(at your option) any later version.
17
18	This program is distributed in the hope that it will be useful,
19	but WITHOUT ANY WARRANTY; without even the implied warranty of
20	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21	GNU General Public License for more details.
22
23	You should have received a copy of the GNU General Public License
24	along with this program. If not, see <https://www.gnu.org/licenses/>.
25	"""
26
27
28	import numpy as np	3✔
29	import scipy.optimize as opt	3✔
30
31
32	def cone_fit(points, origin, direction, angle):	3✔
33	a = origin	×
34	p = points	×
35	n = direction	×
36	height = np.einsum('...i,i->...', (a - p), n)	×
37	distance = np.linalg.norm((a-p) - (height[:, np.newaxis])*n, axis=1)	×
38	expected_distance = -np.tan(angle)*height	×
39	return np.linalg.norm((distance-expected_distance))	×
40
41
42	def create_cone(ncircum, naxial, height, angle, origin, direction):	3✔
43	circum_angles = np.linspace(0, 2*np.pi, ncircum, endpoint=False)	×
44	heights = np.linspace(0, height, naxial+1)[1:]	×
45	points = []	×
46	for i, cangle in enumerate(circum_angles):	×
47	for j, h in enumerate(heights):	×
48	radius = np.tan(angle)*h	×
49	points.append([np.sin(cangle)radius, np.cos(cangle)radius, h])	×
50	points = np.array(points)	×
51	rotation_matrix_z = direction/np.linalg.norm(direction)	×
52	rotation_matrix_y = np.cross(rotation_matrix_z, [1.0, 0.0, 0.0])	×
53	if np.linalg.norm(rotation_matrix_y) < 1e-10:	×
54	rotation_matrix_x = np.cross([0.0, 1.0, 0.0], rotation_matrix_z)	×
55	rotation_matrix_x /= np.linalg.norm(rotation_matrix_x)	×
56	rotation_matrix_y = np.cross(rotation_matrix_z, rotation_matrix_x)	×
57	rotation_matrix_y /= np.linalg.norm(rotation_matrix_y)	×
58	else:
59	rotation_matrix_y /= np.linalg.norm(rotation_matrix_y)	×
60	rotation_matrix_x = np.cross(rotation_matrix_y, rotation_matrix_z)	×
61	rotation_matrix_x /= np.linalg.norm(rotation_matrix_x)	×
62	R = np.array((rotation_matrix_x, rotation_matrix_y, rotation_matrix_z)).T	×
63	new_points = ((R@points.T).T+origin)	×
64	return new_points	×
65
66
67	def cone_error_fn_fixed_angle(points, angle):	3✔
68	return lambda x: cone_fit(points, x[:3], x[3:], angle)	×
69
70
71	def cone_error_fn_free_angle(points):	3✔
72	return lambda x: cone_fit(points, x[:3], x[3:6], x[6])	×
73
74
75	def fit_cone_fixed_angle(points, angle):	3✔
76
77	def constraint_eqn(x):	×
78	return np.linalg.norm(x[3:])	×
79
80	opt_fn = cone_error_fn_fixed_angle(points, angle)	×
81	origin_guess = np.mean(points, axis=0)	×
82	u, s, vh = np.linalg.svd(points-origin_guess)	×
83	direction_guess = vh[0]	×
84	x0_1 = np.concatenate((origin_guess, direction_guess))	×
85	x0_2 = np.concatenate((origin_guess, -direction_guess))	×
86	# Constrain to a unit vector
87	constraint = opt.NonlinearConstraint(constraint_eqn, 1.0, 1.0)	×
88	result_1 = opt.minimize(opt_fn, x0_1, method='trust-constr', constraints=constraint)	×
89	result_2 = opt.minimize(opt_fn, x0_2, method='trust-constr', constraints=constraint)	×
90	# Figure out which is the better one
91	result = result_1 if result_1.fun < result_2.fun else result_2	×
92	cone_origin = result.x[:3]	×
93	cone_direction = result.x[3:]	×
94	return result.fun, cone_origin, cone_direction	×
95
96
97	def distance_point_line(points, origin, direction):	3✔
98	return np.linalg.norm(np.cross((points-origin), direction), axis=-1)/np.linalg.norm(direction)	×
99
100
101	def distance_point_plane(point, plane_point, plane_direction):	3✔
102	return abs(np.einsum('...i,...i->...', point-plane_point, plane_direction))/np.linalg.norm(plane_direction, axis=-1)	×
103
104
105	def cylinder_fit(points, origin, direction, radius):	3✔
106	return np.linalg.norm(distance_point_line(points, origin, direction) - radius)	×
107
108
109	def fit_cylinder(points, origin_guess=None, direction_guess=None, radius_guess=None):	3✔
110	def opt_fn(x):	×
111	return cylinder_fit(points, x[:3], x[3:6], x[6:7])	×
112
113	def constraint_eqn(x):	×
114	return np.linalg.norm(x[3:6])	×
115
116	if origin_guess is None:	×
117	origin_guess = np.mean(points, axis=0)	×
118	if direction_guess is None or direction_guess == 'largest':	×
119	u, s, vh = np.linalg.svd(points-origin_guess)	×
120	direction_guess = vh[0]	×
121	elif direction_guess == 'smallest':	×
122	u, s, vh = np.linalg.svd(points-origin_guess)	×
123	direction_guess = vh[-1]	×
124	if radius_guess is None:	×
125	u, s, vh = np.linalg.svd(points-origin_guess)	×
126	radius_guess = s[0]/2	×
127	x0 = np.concatenate((origin_guess, direction_guess, [radius_guess]))	×
128	# Constrain to a unit vector
129	constraint = opt.NonlinearConstraint(constraint_eqn, 1.0, 1.0)	×
130	result = opt.minimize(opt_fn, x0, method='trust-constr', constraints=constraint)	×
131	origin = result.x[:3]	×
132	direction = result.x[3:6]	×
133	radius = result.x[6]	×
134	return result.fun, origin, direction, radius	×
135
136
137	def fit_line_point_cloud(points):	3✔
138	center = np.mean(points, axis=0)	×
139	shifted_points = points-center	×
140	cov = points.T@shifted_points	×
141	evals, evects = np.linalg.eigh(cov)	×
142	direction = evects[:, np.argmax(evals)]	×
143	direction /= np.linalg.norm(direction)	×
144	return lambda t: center+t*direction	×
145
146	def intersection_point_multiple_lines(P0,P1):	3✔
147	"""Computes the intersection point of multiple lines in a least squares sense
148
149	P0 and P1 are points on a line. Each are NxD where N is the number of lines
150	and D is the dimensionality of the space (D=2 corresponds to lines on a plane,
151	D=3 is points in 3D space).
152	"""
NEW 153	n = (P1-P0)/np.linalg.norm(P1-P0,axis=1)[:,np.newaxis]	×
NEW 154	projectors = np.eye(n.shape[1]) - n[:,:,np.newaxis]n[:,np.newaxis] # I - nn.T	×
NEW 155	R = projectors.sum(axis=0)	×
NEW 156	q = (projectors @ P0[:,:,np.newaxis]).sum(axis=0)	×
NEW 157	p = np.linalg.lstsq(R,q,rcond=None)[0]	×
158
NEW 159	return p	×

sandialabs / sdynpy / 17046961068

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous