TBC Macros and Extensions

 View Only
Back to discussions
Expand all | Collapse all

Precision Challenges with Perpendicular Dot Product (perpDot) in TBC

  • 1.  Precision Challenges with Perpendicular Dot Product (perpDot) in TBC

    Posted 12-16-2024 12:03
    Edited by Vitalii Vovk 12-16-2024 19:09

    For the past couple of weeks, I've been exploring the dot product and perpendicular dot product (perpDot) practical features .   The perpDot product, in particular, has some unique applications, such as:

    • Determining if two lines are parallel.
    • Identifying which side of the line a point lies on.
    • Finding the intersection point between two line segments.

    While figuring out the math behind these operations was relatively straightforward, I've hit a wall regarding achieving consistent precision for the perpDot calculations. This issue becomes especially noticeable when comparing small segments (less than 0.1') with large segments (more than 1000').

    I have a questions to the  developers and others experienced with TBC's coding:
    What is the best approach to calculate the perpDot product when precision is critical?

    1. Should I rely on rounding or truncating numbers to improve precision?
    2. Are there any built-in TBC tools, settings, or practices that could help ensure accuracy in such calculations?

    I've tried different ways  to make the perpDot calculation more accurate and consistent:

    • Using native TBC tools to calculate perpDot.
    • Truncating input values (e.g., limiting decimal places for point coordinates).
    • Rounding numbers to different levels of precision.

    These attempts have improved the results to some extent, but the precision still isn't ideal for my needs.

    Any pointers and help will be appreciated. Also, I share some of the code that was given the most desirable results among all of the other options (code might not be efficient or redundant in some places)

    
import clr
clr.AddReference('IronPython.Wpf')
import wpf
from System import Math, Action, Func, MidpointRounding
from System.Text import Encoding
from System.Collections.Generic import List
from System.Windows import Visibility
from System.Windows.Controls import StackPanel
from System.Windows.Input import Keyboard, AccessKeyManager
from System import Type, Tuple, Array, Double
from System.IO import StreamReader, MemoryStream
from System.Windows.Threading import Dispatcher
clr.AddReference("System.Drawing")
from System.Drawing import Bitmap
clr.AddReference("System.Windows.Forms")
from System.Windows.Forms import Control, Keys
 
try:
    from Trimble.Vce.ViewModel.CommandLine.CommandLineCommands import OffsetLineCalculator
except:
    clr.AddReference("Trimble.Sdk.UI.Commands")
    from Trimble.Sdk.UI.Commands.CommandLine import OffsetLineCalculator
 
 
try:
    clr.AddReference("Trimble.Sdk") # In version 5.50, model assemblies are "pre-referenced"
except:
    clr.AddReference("Trimble.Vce.Core")
    clr.AddReference("Trimble.Vce.Alignment")
    clr.AddReference("Trimble.Vce.ForeignCad")
    clr.AddReference("Trimble.Vce.Data.Construction")
    clr.AddReference("Trimble.Vce.Geometry")
    clr.AddReference("Trimble.Vce.Units")
    clr.AddReference("Trimble.Vce.Interfaces")
    clr.AddReference("Trimble.DiskIO")
    clr.AddReference("Trimble.Vce.Gem")
 
from Trimble.Vce.Core import TransactMethodCall, ModelEvents
from Trimble.Vce.Core.Components import WorldView, Project, Layer, TextStyle
from Trimble.Vce.SurveyCAD import ProjectLineStyle
from Trimble.DiskIO import SearchPath, OptionsManager
 
from Trimble.Vce.Alignment.Linestring import Linestring,ElementFactory, IStraightSegment, IXYZLocation,ITangentArcSegment, ArcTangentType
from Trimble.Vce.Alignment.Parameters import Bearing, Elevation
 
from Trimble.Vce.Data.Construction.Settings import ConstructionCommandsSettings
 
from Trimble.Vce.Geometry import Point3D, Side, Vector3D
 
from Trimble.Vce.Interfaces import Client
clr.AddReference("Trimble.Vce.UI.BaseCommands")
from Trimble.Vce.UI.BaseCommands import ViewHelper
clr.AddReference("Trimble.Vce.UI.UIManager")
from Trimble.Vce.UI.UIManager import UIEvents, TrimbleOffice
clr.AddReference("Trimble.Vce.UI.Controls")
# later version of TBC moved these to a different assembly. If not found in new location, look at old
 
#This part of the script was replaced to be able to start it in TBC 5.9 and up
try:
    from Trimble.Sdk.UI import InputMethod, MousePosition, CursorStyle, UIEventArgs
    from Trimble.Sdk.Interfaces.UI import ControlBoolean
except:
    try:
        from Trimble.Sdk.Interfaces.UI import InputMethod, MousePosition, CursorStyle, ControlBoolean, UIEventArgs
    except:
        try:
            from Trimble.Vce.UI.UIManager import UIEventArgs,TrimbleOffice
            # version 5.40
            from Trimble.CustomControl.Interfaces.Enums import InputMethod
            clr.AddReference("Trimble.Vce.ViewModel")    
            from Trimble.CustomControl.Interfaces import MousePosition, ControlBoolean
            from Trimble.CustomControl.Interfaces.Enums import CursorStyle
        except:
            from Trimble.Vce.UI.Controls import MousePosition, CursorStyle, ControlBoolean,VceMessageBox
    
            
 
 
# my added imports 
 
from Trimble.Vce.Geometry.PolySeg.Segment import Segment, Line
from Trimble.Vce.Geometry import Point3D, Line2D, Side, Vector3D, Intersections
from Trimble.Vce.ForeignCad import TextEntity, Leader, MText, Point as CadPoint, TextUtilities, AttachmentPoint
from System.Text.RegularExpressions import Regex, Match, RegexOptions
from System import Decimal
from operator import itemgetter
from Trimble.Vce.Interfaces.SnapIn import IPolyseg
#from Trimble.Vce.UI.ConstructionCommands import EditTextCmd
from Trimble.Vce.Interfaces.Units import LinearType
from Trimble.Vce.UI.UIManager import TrimbleOffice, UIEvents
from Trimble.Vce.UI.Controls import VceMessageBox
 
 
from Trimble.Vce.UI.Controls.Wpf import OffsetEdit, SelectEntity, VerticalAngleEdit, TemplateSlope, DistanceEdit, ElevationEdit
from System import String
from System.Windows.Documents import Hyperlink
from System.Diagnostics import Process
 
from System.IO import FileStream, FileMode, FileAccess, Path
from System.Windows.Markup import XamlReader
from System.Windows.Controls import Label
from System.Windows.Input import Key
import time
import threading
from Trimble.Vce.Geometry.PolySeg.Node import Node
 
from Trimble.Vce.Alignment.Linestring.Linestring import ExtendType
 
from Trimble.Vce.Interfaces.SnapIn import ISnapIn
 
import os
import sys
parent_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.insert(0, parent_dir)  # Add the parent directory to sys.path
from TeichertUtilities import VectorsAngle, JoinPadLines, Units, Utilities, Validator
 
 
 
 
from Trimble.Vce.Coordinates import Point as CoordPoint
 
from Trimble.Vce.Core.Components import WorldView 
from System import UInt32
from Trimble.Vce.Core.Components.TextStyle import TextHeightMode, TextJustification
from Trimble.Vce.Interfaces.SnapIn import IPolyseg, IViewMember, ISettableName, IName, IViewSite, ISnapInHelpers
 
 
 
import wpf
from System.Collections.ObjectModel import ObservableCollection
from System.ComponentModel import INotifyPropertyChanged, PropertyChangedEventArgs
from System.Windows.Controls import StackPanel, Button
from System.Windows.Media.Imaging import BitmapImage
from System.Windows import Application, Window
from System.IO import StreamReader
from System import Array, Exception
 
# Data Model for grid view
from System.Windows.Controls import DataGrid, DataGridTextColumn, DataGridTemplateColumn, Button
from System.Windows import DataTemplate, Style
from System.Windows.Controls.Primitives import ButtonBase
from System.Windows import FrameworkElementFactory
from System.Windows.Data import Binding
from System.Windows import RoutedEventArgs  # Correct import for RoutedEventArgs
from System.Collections.ObjectModel import ObservableCollection
from System.Windows import Setter
from System.Windows import DependencyProperty
 
from System.Windows import RoutedEventHandler
from System.Windows.Controls import DataGridLength
 
from System import Uri, UriKind
from System.Windows.Controls import Image
from System.Windows import UIElement
from System.Windows.Media import Stretch 
from System.Windows import Thickness
from System.Windows.Data import RelativeSource
from System.Windows import HorizontalAlignment, VerticalAlignment
 
from Trimble.Vce.Geometry.PolySeg import PolySeg
from decimal import Decimal as pDecimal, getcontext #important to change name to Decimal not get confused with .net decimals
 
 
 
# addin this vector that will be inherited by custom vector class for the purpose of the instance check
class CustomVector:
    pass
 
class SegmentsSettings:
    def __init__(self, **kwargs):
        # Validate and initialize attributes
        self._container = kwargs.get("container", None)
        self._segment = kwargs.get("segment", None)
        self._segmentA = kwargs.get("segmentA", None)
        self._segmentB = kwargs.get("segmentB", None)
        self._point = kwargs.get("point", None)
        
 
        # Validate container at initialization
        if self._container is not None and not isinstance(self._container, WorldView):
            raise TypeError("Passed container is empty or incorrect")
 
    @property
    def container(self):
        return self._container
 
    @container.setter
    def container(self, value):
        # Validate container type
        if value is None or not isinstance(value, WorldView):
            raise TypeError("Passed container is empty or incorrect")
        self._container = value
 
    
    @property
    def vPoly(self):
        return self._vPoly
 
    @vPoly.setter
    def vPoly(self, value ):
        test = clr.GetClrType(PolySeg)
        if not isinstance(value, PolySeg):
            TypeError("object should be type of vertical PolySeg")
        self._vPoly = value
            
 
 
    @property
    def point(self):
        return self._point
 
    @point.setter
    def point(self, value):
        if value is not None and not isinstance(value, Point3D):
            raise TypeError("point must be type of Point3D")
        self._point = value
 
    @property
    def segment(self):
        return self._segment
 
    @segment.setter
    def segment(self, value):
        # Validate type using isinstance
        if value is not None and not isinstance(value, Segment):
            raise TypeError("segment must be of type Segment.")
        self._segment = value
 
    @property
    def segmentA(self):
        return self._segmentA
 
    @segmentA.setter
    def segmentA(self, value):
        if value is not None and not isinstance(value, Segment):
            raise TypeError("segmentA must be of type Segment.")
        self._segmentA = value
 
    @property
    def segmentB(self):
        return self._segmentB
 
    @segmentB.setter
    def segmentB(self, value):
        if value is not None and not isinstance(value, Segment):
            raise TypeError("segmentB must be of type Segment.")
        self._segmentB = value
 
 
class CustomSegment(SegmentsSettings):
    def __init__(self, **kwargs):
        
        super(CustomSegment, self).__init__(**kwargs)
 
    def overlappedSegments(self, segA, segB):
        """
        This method calculates the overlap properties between two segments (Segment A and Segment B).
 
        Segment A is considered the primary segment, and all calculations are based on its start point.
 
        **Key Outputs:**
        - `overlap`: Boolean, True if any part of Segment B overlaps with Segment A.
        - `tB0`, `tB1`: t-values of Segment B's start and end points projected onto Segment A AKA  ParametricLocation.
        - `p_B0_overlap`: Boolean, True if Segment B's start point is within the bounds of Segment A.
        - `p_B1_overlap`: Boolean, True if Segment B's end point is within the bounds of Segment A.
        - `stationB0`: The station on Segment A corresponding to Segment B's start point if overlapping.
        - `stationB1`: The station on Segment A corresponding to Segment B's end point if overlapping.
        - `segBinSegA`: Boolean, True if Segment B is fully inside or equal to Segment A.
        - `segAinSegB`: Boolean, True if Segment A is fully inside of  Segment B.
 
        **Notes:**
        - If Segment A is fully inside Segment B, `stationB0`, `stationB1`, `p_B0_overlap`, and `p_B1_overlap` will remain None.
        - Calculations use canonical line equations {x=mt+x0; y=mt+y0} or the equivalent gA(s) = A0 + t * (A1 - A0).
        - Parallelism is determined using `PerpDot` (Perpendicular Dot Product). If segments are parallel, they are checked for collinearity.
 
  
        """
 
        if segA is None or segB is None:
            raise ValueError("segmentA and segmentB should be instances of Segment in overlappedSegments method.")
 
        # Assign segments
        segA = self.segment = segA
        segB = self.segment = segB
 
 
        if segA.IsZeroLength2D():
            raise ValueError("Segment A Length 2D is 0")
        if segB.IsZeroLength2D():
            raise ValueError("Segment A Length 2D is 0")
 
 
        #!!!!! need to check if segments are not 0 length !!!!!!!!!!
 
        # Initialize 2D vectors for calculations (to ensure compatibility)
 
        # Direction vector for Segment A and B
        vA = self.__calculateVectorPrecise(segA.BeginPoint, segA.EndPoint,  imitateVector = True)
        #vAt = segA.EndPoint - segA.BeginPoint # for testing purposes only
       
 
        vB = self.__calculateVectorPrecise(segB.BeginPoint, segB.EndPoint, imitateVector = True)
        #vBt = segB.EndPoint - segB.BeginPoint # for testing purposes only
 
 
        #Vector from start of Segment A to Start of Segment B
        vU_A0_B0 = self.__calculateVectorPrecise(segA.BeginPoint, segB.BeginPoint,  imitateVector = True )
         #Vector from start of Segment A to End of Segment B
        vU_A0_B1 = self.__calculateVectorPrecise(segA.BeginPoint, segB.EndPoint,  imitateVector = True )
 
 
        # Initialize return values
        overlap = False  # Flag indicating if any overlap exists
        tB0 = None       # t-value for Segment B's start point on Segment A
        tB1 = None       # t-value for Segment B's end point on Segment A
        p_B0_overlap = False  # True if Segment B's start is within Segment A
        p_B1_overlap = False  # True if Segment B's end is within Segment A
        stationB0 = None      # Station of Segment B's start point on Segment A
        stationB1 = None      # Station of Segment B's end point on Segment A
        segAinSegB = False    # True if Segment A is fully inside Segment B
        segBinSegA = False    # True if Segment B is fully inside Segment A
 
        # Default return values
        defaultValues = (overlap, tB0, tB1, p_B0_overlap, p_B1_overlap, stationB0, stationB1, segBinSegA, segAinSegB)
 
        # Check if the segments are collinear (parallel)
 
   
 
        if not self.__calculatePerpDotPrecise(vA, vB, normalize = True, autothreshHold=True) == 0: # is 15 16 digits prcision for doubles in c# ?
        #if not self.simpleCollinear(vA, vB) : # is 15 16 digits prcision for doubles in c# ?
            return defaultValues  # Segments are not parallel
 
        # Check if the segments are collinear, if true it also means they are located on the same line 
        if not self.__calculatePerpDotPrecise(vA, vU_A0_B0, normalize=True, autothreshHold=True) == 0:
        #if not self.simpleCollinear(vA, vU_A0_B0):
            return defaultValues  # Segments do not coincide on the same line
 
        # Calculate t-values for Segment B's start and end points projected onto Segment A
        tB0 = self.__parametricValuePrecise(vU_A0_B0,vA, roundValueDigits = 9)
        tB1 = self.__parametricValuePrecise(vU_A0_B1,vA , roundValueDigits = 9)
 
       
 
        # Check if Segment B's start point is on Segment A
        if 0 <= tB0 <= 1:
            p_B0_overlap = True  # Overlap exists
            stationB0 = float(pDecimal(str(segA.BeginStation)) + (pDecimal(str(segA.EndStation)) *  pDecimal(str(tB0)))) # the same as  segA.Station(tB0)      
      
 
        # Check if Segment B's end point is on Segment A
        if 0 <= tB1 <= 1:
            p_B1_overlap = True  # Overlap exists
            try:
               stationB1 = float(pDecimal(str(segA.BeginStation)) + pDecimal(str(segA.EndStation)) * pDecimal(str(tB1)))
 # the same as segA.Station(tB1)
            except Exception as err:
            # Handle errors and display an error message.
                VceMessageBox.Show(TrimbleOffice.TheOffice.MainWindow.Form, "Error: {}".format(err))
 
       # Determine overlap and containment relationships
        if p_B0_overlap and p_B1_overlap:
            # Segment B is fully within Segment A
            if (tB0 in {0.0, 1.0}) and (tB1 in {0.0, 1.0}): # checking if segment end points are indentical
                  segAinSegB = True
            segBinSegA = True
            overlap = True
        elif not p_B0_overlap and not p_B1_overlap:
            if tB0 is not None and tB1 is not None:
                # Segment A is fully within Segment B
 
                if tB0 * tB1 < 0: #Only if to ends with oposite sign it meats the segment A is fully inside of the segment B
                    # NEEDDS TO BE FIXED ONE SIDE CAN EQUALS 0
                    segAinSegB = True
                    overlap = True
            # clear Parametric value for Segment A because end points are outside
            tB0 = None
            tB1 = None
        else:
            # needs to be fixed added if 0 and t should give you 
            if (tB0 in {0.0, 1.0}) or (tB1 in {0.0, 1.0}):
                    if tB0 == 0.0 and (tB0 + tB1) > 0:
                        segAinSegB = True
                    elif tB1 == 0.0 and (tB0 + tB1) > 0:
                        segAinSegB = True
                    elif tB0 == 1.0 and (tB0 * tB1) < 0:
                        segAinSegB = True
                    elif tB1 == 1.0 and (tB0 * tB1) < 0:
                        segAinSegB = True
 
            # Partial overlap: At least one endpoint of Segment B overlaps with Segment A
            overlap = True
            tB0 = tB0 if p_B0_overlap else None
            tB1 = tB1 if p_B1_overlap else None
 
        # Update return values if overlap exists
        if overlap:
            defaultValues = (overlap, tB0, tB1, p_B0_overlap, p_B1_overlap, stationB0, stationB1, segBinSegA, segAinSegB)
        return defaultValues
 
    def simpleCollinear(self, vA, vB, truncateValue=6):
        
        # Case 1: Both vectors are vertical 
        if vA.X == 0 and vB.X == 0:
            return True
        # Case 2: Both vectors are horizontal (slope = 0)
        if vA.Y == 0 and vB.Y == 0:
            return True
        # Case 3: Neither vector is vertical, compare tangents
        if vA.X != 0 and vB.X != 0:
            getcontext().prec = 200
            tanA = Math.Atan(pDecimal(str(vA.Y)) / pDecimal(str(vA.X)))  # Tangent for vA
            tanB = Math.Atan(pDecimal(str(vB.Y)) / pDecimal(str(vB.X)))  # Tangent for vB
            # Truncate values to the specified precision
            tanA = self.truncate_to(tanA, digits=truncateValue)
            tanB = self.truncate_to(tanB, digits=truncateValue)
            # Compare tangents for collinearity
            if tanA == tanB:
                return True
        # Case 4: Non-matching cases
        return False
          
 
 
        
 
        
        
 
 
    
 
    
 
    def __calculateVectorPrecise(self, StartV, EndV, imitateVector = False):
        """
        Parameters:
        StartV, EndV: Objects with attributes X, Y, and Z (representing 3D points).
 
        Returns:
        VectorImitation: A new object representing the vector difference.
        """
        getcontext().prec = 200  # Set precision for Decimal calculations.
 
        
        class VectorImitationPrecise(CustomVector):
            def __init__(self, X, Y, Z, length2D=None):
                """
                Initializes the VectorImitationPrecise instance.
 
                Parameters:
                X, Y, Z (Decimal or compatible): Components of the vector.
                length2D (Decimal, optional): Precomputed 2D length. If None, it is calculated.
                """
                if not all(isinstance(value, (Decimal, int, float, pDecimal)) for value in [X, Y, Z] if value is not None):
                    raise ValueError("X, Y, and Z must be Decimal, int, or float.")
                self.X = pDecimal(X)
                self.Y = pDecimal(Y)
                self.Z = pDecimal(Z) if Z is not None else None
                self.Length2D = self.calcLength2D() if length2D is None else length2D
                self.Length = self.calcLength()
                self.Angle = Math.Atan2(self.Y, self.X)
                self.Azimuth = self.azimuth()
 
            def azimuth(self):
                if self.X != 0:
                    angle = abs(Math.Atan(self.Y/self.X))
                    if self.X >= 0  and self.Y >= 0: #first quadrant
                        return (Math.PI/2) - angle
                    elif self.X < 0 and self.Y >= 0: # second quadrant
                        return (3*Math.PI/2) + angle
                    elif self.X < 0 and self.Y < 0: #third quadrant
                        return (3*Math.PI/2) - angle 
                    elif self.X >= 0 and self.Y < 0:
                        return (Math.PI/2) + angle
 
 
                    
 
 
 
 
            def calcLength(self):
                """
                Calculates the 3D length of the vector with high precision.
                Returns:
                Decimal: The 3D length of the vector
                """
                if self.Z is None:
                    return self.sqrt_precise(self.sum_of_squares(self.X, self.Y))
 
                return self.sqrt_precise(self.sum_of_squares(self.X, self.Y, self.Z))
 
            def Clone(self):
               """
               Creates a copy of the current instance.
               
               Returns:
               VectorImitationPrecise: A new instance with the same attribute values.
               """
               return VectorImitationPrecise(self.X, self.Y, self.Z, length2D=self.Length2D)
 
 
            def calcLength2D(self):
                """
                Calculates the 2D length of the vector (ignoring Z) with high precision.
                Returns:
                Decimal: The 2D length of the vector.
                """
                return self.sqrt_precise(self.sum_of_squares(self.X, self.Y))
 
            def sum_of_squares(self, *values):
                """
                Computes the sum of squares for multiple components with high precision.
 
                Parameters:
                components (Decimal): Components for which to calculate the sum of squares.
 
                Returns:
                Decimal: The sum of squares of the components.
                """
                result = sum(value * value for value in values)
                return result
 
            def sqrt_precise(self, value, precision=200):
                if value < 0:
                    raise ValueError("Cannot compute the square root of a negative number.")
                if value == 0:
                    return pDecimal(0)
 
                getcontext().prec = precision
                guess = value / 2
                epsilon = pDecimal('1e-{0}'.format(precision // 2))
 
                while True:
                    next_guess = (guess + value / guess) / 2
                    if abs(next_guess - guess) < epsilon:
                        return next_guess
                    guess = next_guess
 
            def __add__(self, object):
                """
                Adds a Point3D object to the current Point3D instance.
 
                Parameters:
                object (Point3D): The Point3D object to add.
 
                Returns:
                Point3D: A new Point3D object representing the sum of the two points.
                """
                if isinstance(object, clr.GetClrType(Point3D)):
                    getcontext().prec = 200
                    x_Value = self.X + pDecimal(str(object.X))
                    y_Value = self.Y + pDecimal(str(object.Y))
                    z_Value = (
                        self.Z + pDecimal(str(object.Z))
                        if not Double.IsNaN(object.Z) and self.Z is not None
                        else Double.NaN
                    )
                    return Point3D(x_Value, y_Value, z_Value)
 
            def __mul__(self, value):
                """
                Multiplies the current Point3D instance by a scalar value.
            
                Parameters:
                value (Double or Decimal): The scalar value to multiply by.
            
                Returns:
                VectorImitation: A new VectorImitation object representing the scaled point.
                """
                if isinstance(value, (Double, Decimal)):
                    getcontext().prec = 200
                    scalar = pDecimal(str(value))
                    x_Value = scalar * self.X
                    y_Value = scalar * self.Y
                    z_Value = scalar * self.Z if self.Z is not None else None
                    return VectorImitationPrecise (x_Value, y_Value, z_Value)
                else:
                    raise TypeError("Unsupported operand for ImitateVectorPrecise")
 
            def __rmul__(self, scalar):
                """
                Multiplies the vector by a scalar value (scalar * vector).
 
                Parameters:
                scalar (float or Decimal): The scalar value.
 
                Returns:
                VectorImitationPrecise: A new scaled vector.
                """
                return self.__mul__(scalar)
                   
                    
                    
 
 
        try:
 
            # Calculate precise differences for X, Y, and Z components.
            x_diff = pDecimal(str(EndV.X)) - pDecimal(str(StartV.X))
            y_diff = pDecimal(str(EndV.Y)) - pDecimal(str(StartV.Y))
            
            if Double.IsNaN(StartV.Z) or Double.IsNaN(EndV.Z):
                z_diff = None
            else:
                z_diff = pDecimal(str(EndV.Z)) - pDecimal(str(StartV.Z))
 
            #test = float(VectorImitationPrecise(pDecimal(54.706244987), pDecimal(84.603718944), None).length2D)
 
            # Create and return the result vector.
            if imitateVector:
                return VectorImitationPrecise(x_diff, y_diff, z_diff)
            return Vector3D(float(x_diff), float(y_diff), Double.NaN if z_diff is None else float(z_diff))
 
        except Exception as err:
            # Handle errors and display an error message.
            VceMessageBox.Show(TrimbleOffice.TheOffice.MainWindow.Form, "Error: {}".format(err))
            return None
 
    def __calculatePerpDotPrecise(self, vA, vB, autothreshHold=False, normalize = False):
        """
        Calculates the perpendicular dot product (PerpDot) of two vectors with high precision.
        Parameters:
        vA (Vector3D): The first vector.
        vB (Vector3D): The second vector.
        precision (int): The precision to use for the pDecimal calculations.
        Returns:
        pDecimal: The PerpDot product of the vectors with high precision.
        """
        getcontext().prec = 200
        if normalize:
            vA = self.__normalizedPrecice(vA)
            vB = self.__normalizedPrecice(vB)
            
        if autothreshHold:
            minL = vA.Length2D if vA.Length2D < vB.Length2D else vB.Length2D
            treshHold = self.__setThreshHold(minL)
 
        perpDot  = (pDecimal(str(vA.X)) *  pDecimal(str(vB.Y))) - (pDecimal(str(vA.Y)) * pDecimal(str(vB.X)))
        if autothreshHold:
            if abs(perpDot) <= treshHold:
               return pDecimal(0)
 
        return perpDot
 
 
    def __setThreshHold(self, value):
        """
        Determines the appropriate threshold based on the input value.
        Iterates through thresholds in ascending order of keys.
 
        Parameters:
        value (float): The input value to determine the threshold.
 
        Returns:
        float: The appropriate threshold for the given value.
        """
        thresholds = {
            0.005: 1E-5,
            0.01: 1E-6,
            0.1: 1E-7,
            1: 1E-8,
            10: 1E-9,
            100: 1E-10,
            10000: 1E-11,
            100000: 1E-12,
            1000000: 1E-13,
        }
 
        # Iterate over keys in ascending order
        for key in sorted(thresholds.keys()):
            if value <= key:
                return thresholds[key]
        return 1E-14
 
 
 
 
    
    def __normalizedPrecice(self, v, tresHold=pDecimal("1E-8")):
        """
        Normalizes the vector `v` with high precision, setting its 2D length to 1.
 
        Parameters:
        v (Vector2D): The vector to normalize.
        tresHold (Decimal): The threshold below which the vector is considered too small to normalize.
 
        Returns:
        Vector2D: The normalized vector with its length set to 1.
        """
        v = v.Clone()
 
        # Check if the vector is effectively zero
        if isinstance(v, CustomVector):
            if pDecimal(str(v.Length2D)) < tresHold:
                v.Length2D = pDecimal(0)
                v.X = pDecimal(0)
                v.Y = pDecimal(0)
                if v.Z:
                   v.Z = pDecimal(0)
                return v
 
            # Set precision for Decimal calculations
            getcontext().prec = 200
 
            # Normalize the vector
            v.X = pDecimal(str(v.X)) / pDecimal(str(v.Length2D))
            v.Y = pDecimal(str(v.Y)) / pDecimal(str(v.Length2D))
            v.Length2D = pDecimal(1)  # Set normalized length only works for imitate vector
            return v
        else:
            return v.Normalized()
 
 
 
 
 
    def __parametricValuePrecise(self, firstV, secondV, precision=200, roundValue=True, roundValueDigits = 15, returnFloat=True):
        """
        Calculates precise parametric values using high precision with the Decimal module.
        Parameters:
        firstV (Vector3D): The vector representing the numerator.
        secondV (Vector3D): The vector representing the denominator.
        Returns:
        float or Decimal: The parametric value, rounded to the specified precision.
        """
        test1 = float(firstV.X)
        test2 = float(secondV.X)
        result = test1/test2
 
        # Set the precision for Decimal calculations
        getcontext().prec = precision
 
        try:
            # Calculate parametric value using X components
            t = pDecimal(str(firstV.X)) / pDecimal(str(secondV.X))
        except (ZeroDivisionError, InvalidOperation):
            try:
                # Fallback to Y components if X components fail
                t = pDecimal(str(firstV.Y)) / pDecimal(str(secondV.Y))
            except (ZeroDivisionError, InvalidOperation):
                # Should not occur if the vectors are valid
                raise ZeroDivisionError("Division by zero encountered for both X and Y components of the second vector.")
 
        # Round the result to the specified number of decimal places
        if roundValue:
            rounding_format = '1e-' + str(roundValueDigits)
            t = t.quantize(pDecimal(rounding_format))
        if returnFloat:
            t = float(t)
 
        # Return as float or Decimal based on user preference
        return t
 
 
    def pointSide(self, segment, point, thresHold = float("1E-7")):
        """
        it returns the point is to the segment
        """
        seg = self.segment = segment
        point = self.point = point
 
        v = seg.EndPoint - seg.BeginPoint
        vP = point - seg.EndPoint
 
        vt = self.__calculateVectorPrecise(seg.BeginPoint, seg.EndPoint)
        vPt = self.__calculateVectorPrecise(seg.EndPoint, point)
        perpDot = Math.Round(Vector3D.PerpDot(v, vP), 10)
        test = self.__calculatePerpDotPrecise(vt, vPt)
        test = float(test)
 
        # if may not catch if point directly on the line 
        if abs(perpDot) < thresHold:
            return Side.On
        else:
            return Side.Right if perpDot < 0 else Side.Left
    
 
 
 
    def segmentInersection(self, segA, segB,  elevationMode = 1):
        """
        Base for finding intersection is canonical formula for a line 
        in a form {x= mt + x0 y= nt + y0} that can be writen as in a matrix form
        A  = at + A0 if two line have the solution that is the point of itersection
        as+A0 = bt+B0 and multiply one of the vectors on perpDot of itself that would give us 0
 
         **Key Outputs:**
        - point that represent point of intesection
        - s parametric value (scalar for vector a to point of intersection) for segemenA
        - t parametric value for sefment B 
 
        **Notes:**
        - if parametic valuw is from 0 to 1 it means point of intersection on that segment or
          both of the segments if both parametric values are in the range from 0 to 1
        - This method will work well horizontal nodes. If line has horizontal Nodes and VPI
          you woudl need to convert VPI to horizontal nodes to make this method give correct ele
          vation
 
        **Elevation Mdodes**
        - 1 - return intersection point with grade projection based on Segment A
        - 2 - return intersection point with grade projection based on Segment B
        - 3 - return intersection point with grade that would mean elevation between
              segment A and segment B
        - 4 - returns minimul elevation between segment A and Segment B
        - 5 - returns maximum elevation between Segment A and Segment B
 
        

        """
 
 
 
        segA = self.segment = segA #run segment throught the setter
        segB = self.segment = segB  #run segment throught the setter
 
 
        vA = self.__calculateVectorPrecise(segA.BeginPoint, segA.EndPoint) # create a vector based on segment 
        vB = self.__calculateVectorPrecise(segB.BeginPoint, segB.EndPoint) # create a vector based on segment 
        isParallel = self.simpleCollinear(vA, vB,  truncateValue=6) # make sure this segmens are not parallel
        if  isParallel:
            return False
 
        vU_A0_B0 = self.__calculateVectorPrecise(segA.BeginPoint, segB.BeginPoint) #addtiional vector to make formula work
        vU_B0_A0 = self.__calculateVectorPrecise(segB.BeginPoint, segA.BeginPoint)
 
 
 
        #computation for segment A 
        s = self. __calculatePerpDotPrecise(vB,vU_A0_B0 ) / self. __calculatePerpDotPrecise(vB,vA ) # parametric Value for segment A 
        s = float(s)
 
        # have to flip in order to use my custom __add__ (it meters on which side '+' is )
        #it works faster if you dont use imitateVectorPrecice (CustomVector)
 
        if isinstance(vA, CustomVector):
            pA =  (s * vA) + segA.BeginPoint #point with elevation based on Segment A 
        else:
            pA =  segA.BeginPoint + (s * vA)
 
        
 
        #computation for serment B 
        t = self. __calculatePerpDotPrecise(vA,vU_B0_A0 ) / self. __calculatePerpDotPrecise(vA, vB ) # parametric Value parametrivale for Segment B
        t = float(t)
        # have to flip in order to use my custom __add__ (it meters on which side '+' is )
        #it works faster if you dont use imitateVectorPrecice (CustomVector)
        if isinstance(vB, CustomVector):
            pB = (t * vB)  + segB.BeginPoint #point with elevation based on Segment B
        else:
             pB = segB.BeginPoint + (t * vB)
 
 
        point = None
 
        # Adjust resulted point for the selected elevation mode
        if elevationMode == 1:
            point = pA
        elif elevationMode == 2:
            point =  pB
        elif elevationMode == 3:
            meanEl = (pA.Z + pB.Z)/2
            point =  Point3D(pA.X, pA.Y, meanEl )
        elif elevationMode == 4:
            minEl = pA.Z if pA.Z < pB.Z else pB.Z
            point = Point3D(pA.X, pA.Y,minEl )
        elif elevationMode == 5:
            maxEl = pA.Z if pA.Z > pB.Z else pB.Z
            point = Point3D(pA.X, pA.Y, maxEl)
 
        roundedPoint = self.__roundPointValues(point) #round up poitn x,y, and z value
 
        return (roundedPoint, s, t )
 
    def __roundPointValues(self, point, digits=6):
        if isinstance(point, Point3D):
            x = Math.Round(point.X, digits)
            y = Math.Round(point.Y, digits)
            z  = Math.Round(point.Z, digits) if point.Z is not Double.NaN else Double.NaN
            return Point3D(x,y,z)
        return None
 
 
 
 
 
 
 
 
    def truncate_to(self, value, digits=15):
        """
        Truncate a floating-point number to the specified number of decimal places without rounding.
 
        Parameters:
        value (float): The number to be truncated.
        digits (int): The number of decimal places to keep.
 
        Returns:
        float: The truncated number.
        """
        value_str = "%.20f" % value  # Use classic string formatting for high precision
        if '.' in value_str:
            integer_part, decimal_part = value_str.split('.')
            truncated_str = integer_part + '.' + decimal_part[:digits]
            return float(truncated_str)
        return float(value_str)
 
 
    def truncate_entity(self, entity, digits=15):
        """
        Truncate all components of a Vector3D or Point3D to a fixed number of decimal places.
 
        Parameters:
        entity (Vector3D or Point3D): The vector or point to be truncated.
        digits (int): The number of decimal places to keep.
        type (str): The type of object to return ("Vector" or "Point").
 
        Returns:
        Vector3D or Point3D: A new object with truncated components.
        """
        truncated_x = self.truncate_to(entity.X, digits)
        truncated_y = self.truncate_to(entity.Y, digits)
        truncated_z = self.truncate_to(entity.Z, digits)
        
        if isinstance(entity, clr.GetClrType(Vector3D)):
            return Vector3D(truncated_x, truncated_y, truncated_z)
        elif isinstance(entity, clr.GetClrType(Point3D)):
            return Point3D(truncated_x, truncated_y, truncated_z)
        else:
            raise ValueError("Invalid type specified. Use 'Vector' or 'Point'.")



    ------------------------------
    Vitalii Vovk
    ------------------------------



  • 2.  RE: Precision Challenges with Perpendicular Dot Product (perpDot) in TBC

    Posted 12-17-2024 02:24
    Edited by Ronny Schneider 12-17-2024 02:30

    Hi Vitalii,

    even without looking at your code I can see that there will always be an issue to make a 100 % accurate statement whether an extremely short segment is parallel to a very long one.

    The problem I see is the source of the coordinates for those segments. If I assume that they come from some survey activity the best you can hope for is probably 1/10 mm.

    That means the definition of your segments will be accurate to 1/10 over 1 inch and 1/10 over 1000 foot. The spatial definition of the longer segment will be much more accurate than the shorter one.

    I'd probably compute the offset of the shorter segment nodes to the longer segment. And when the difference of the perpendicular offset is smaller than my expected coordinate precision I'd declare the segments parallel.


    I haven't tested your code or run some experiments on my own yet.

    Do you see inconsistent offset results if you declare a long segment (i.e. 0 to 100000) and a short one (i.e. +0.0001 to +0.0002)? Do you get as answer that the short segment may be on the left side of the line?


    ------------------------------
    Ronny Schneider
    ------------------------------

    Original Message


  • 3.  RE: Precision Challenges with Perpendicular Dot Product (perpDot) in TBC

    Posted 12-17-2024 09:08

    I was analyzing CAD line segments and noticed that the input points used to build vectors in TBC might be inherently incorrect, probably due to floating-point precision issues. I was wondering how much I can trust the precision of the input from the Segment Points in TBC.

    You can tell a lot of things by utilizing the perpDot Product for segments.

    The approach that I took on overlapping segments:

    1. I built vectors that correspond to the given segments.

    2. If segments are parallel, the perpDot of their vectors should be 0.

    3. Now we test if the segments are on the same line by comparing the vector of Segment A to a new vector built from any point of Segment A to any point of Segment B. If they are on the same line, the perpDot will be 0.

    4. You have to choose a vector and point to which all results will be referenced. I chose Segment A and its start point.

    5. Now you can build additional vectors from the start point of Segment A: one to the start point of Segment B and another to the end point of Segment B.

    6. By finding out the scalar between the vector of Segment A and, for example, the vector from the start point of Segment A to the start point of Segment B (simply dividing one of the coordinates by the other), you will get the parametric value for the start point of Segment B in relation to the start point of Segment A. If it is between 0 and 1, the start point of Segment B is within Segment A. If it is 0, it is the same as the start point of Segment A. If it equals 1, it is the same point as the end of Segment A. If it is less than 0, it is before the start point of Segment A. If it is more than 1, it is past the limits of Segment A in the direction of the vector from Segment A.



    ------------------------------
    Vitalii Vovk
    ------------------------------

    Original Message


  • 4.  RE: Precision Challenges with Perpendicular Dot Product (perpDot) in TBC

    Posted 01-13-2025 19:58

    Where exactly do those CAD line segments originate from? Another CAD package, survey data, or based on some coordinates and then rotated/offset/shifted all the time?

    Arguing about the 7th or 8th decimal is completely dependent on the sources origin/accuracy.

    You say that the point information you use to define a Vector3D might be incorrect? How do you come to this conclusion? Did you compare the Point3D information, that you've retrieved from a polyseg/segment???, and that doesn't match the data you see in another package.

    Please upload some sample data?



    ------------------------------
    Ronny Schneider
    ------------------------------

    Original Message


  • 5.  RE: Precision Challenges with Perpendicular Dot Product (perpDot) in TBC

    Posted 01-16-2025 20:56

    My source data comes from simply drawing two lines (each with only two points: start and stop) in TBC at 90 degrees to each other.

    If I convert those lines into vectors, the dot product of the vectors should theoretically be zero. However, when the lines are based around large coordinates, such as 2,000,000 and 6,000,000 (plus some 15-17 decimal points ), I was never able to achieve a dot product result of exactly zero.

    I understand the floating-point precision issues  and I tried everything I could think of to maintain consistent precision in the calculations. At that point, I wanted to determine when the dot product would be zero and when it wouldn't.

    I did have some success figuring out overlapping segments in about 80% of cases, but the solution was not 100% reliable, which is frustrating. The reliability also seemed to depend on the length of the vectors, especially in cases where one vector was very large (e.g., 1 million meters) and the other was very small (e.g., less than 1 cm) (I did tey to normalize vectors as well).  I had better success finding the intersection point between two segments as you have a little more room for accuracy ther

    I concluded (though I might be wrong) that the error might originate in the coordinates within TBC itself. For instance, if a coordinate like 656389634.3456353534534553 is being used, it's likely not accurate beyond the 6th or 7th decimal place when using a float type.

    I'm curious: if a dot product method exists in the TBC SDK, what approach do developers use to ensure high accuracy in such calculations when precision is must?



    ------------------------------
    Vitalii Vovk
    ------------------------------

    Original Message


  • 6.  RE: Precision Challenges with Perpendicular Dot Product (perpDot) in TBC

    Posted 01-16-2025 22:20

    Have you tried the Vector3D.DotProduct ?

    The following is really stretching accuracy limits. I can't see a geospatial use case where this becomes critical.

    Here is probably some inaccuracy of the Sin/Cos computation for the rotation kicking in.

    If I shorten v2 after I rotate it that arbitrary amount, I get a DotProduct of 0.

    For geospatial use I'd probably just declare every deviation after the 10th decimal as unimportant.



    ------------------------------
    Ronny Schneider
    ------------------------------

    Original Message