go-hep.org/x/hep@v0.38.1/fastjet/internal/predicates/orientation.go (about)

     1  // Copyright ©2017 The go-hep Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package predicates
     6  
     7  import (
     8  	"fmt"
     9  	"math/big"
    10  
    11  	"gonum.org/v1/gonum/mat"
    12  )
    13  
    14  // OrientationKind indicates how three points are located in respect to each other.
    15  type OrientationKind int
    16  
    17  const (
    18  	// Counterclockwise
    19  	CCW OrientationKind = iota
    20  	// Clockwise
    21  	CW
    22  	Colinear
    23  	IndeterminateOrientation
    24  )
    25  
    26  func (o OrientationKind) String() string {
    27  	switch o {
    28  	case CCW:
    29  		return "Counterclockwise"
    30  	case CW:
    31  		return "Clockwise"
    32  	case Colinear:
    33  		return "Colinear"
    34  	case IndeterminateOrientation:
    35  		return "Indeterminate"
    36  	default:
    37  		panic(fmt.Errorf("predicates: unknown OrientationKind %d", int(o)))
    38  	}
    39  }
    40  
    41  // Orientation returns how the point (x,y) is oriented with respect to
    42  // the line defined by the points (x1,y1) and (x2,y2).
    43  func Orientation(x1, y1, x2, y2, x, y float64) OrientationKind {
    44  	o := simpleOrientation(x1, y1, x2, y2, x, y)
    45  	if o == IndeterminateOrientation {
    46  		// too close to 0 to give a definite answer.
    47  		// Therefore check with more expansive tests.
    48  		o = matOrientation(x1, y1, x2, y2, x, y)
    49  	}
    50  	return o
    51  }
    52  
    53  // simpleOrientation finds the orientation using the simple float64 type.
    54  // Its accuracy can't be guaranteed, therefore close decisions
    55  // return IndeterminateOrientation which signals Orientation that further
    56  // testing is necessary.
    57  //
    58  // It computes the determinant of the matrix and returns the orientation based
    59  // on the value of the determinant.
    60  //
    61  //	| x1 y1 1 |
    62  //	| x2 y2 1 |
    63  //	| x  y  1 |
    64  func simpleOrientation(x1, y1, x2, y2, x, y float64) OrientationKind {
    65  	if (x1 == x2 && x2 == x) || (y1 == y2 && y2 == y) {
    66  		// points are horizontally or vertically aligned
    67  		return Colinear
    68  	}
    69  	// Compute the determinant of the matrix
    70  	p := newFloat64Pred(x1).mulFloat64(y2).addFloat64Pred(newFloat64Pred(x2).mulFloat64(y)).
    71  		addFloat64Pred(newFloat64Pred(x).mulFloat64(y1)).subFloat64Pred(newFloat64Pred(x1).mulFloat64(y)).
    72  		subFloat64Pred(newFloat64Pred(x2).mulFloat64(y1)).subFloat64Pred(newFloat64Pred(x).mulFloat64(y2))
    73  	// det := x1*y2 + x2*y + x*y1 - x1*y - x2*y1 - x*y2
    74  	det := p.n
    75  	// e determines when the determinant in simpleOrientation is too close to 0 to rely on floating point operations.
    76  	e := p.e
    77  	if det > e {
    78  		return CCW
    79  	}
    80  	if det < -e {
    81  		return CW
    82  	}
    83  	return IndeterminateOrientation
    84  }
    85  
    86  // robustOrientation finds the orientation using the accurate big/Rat type.
    87  //
    88  // It computes the determinant of the matrix and returns the orientation based
    89  // on the value of the determinant.
    90  //
    91  //	| x1 y1 1 |
    92  //	| x2 y2 1 |
    93  //	| x  y  1 |
    94  func robustOrientation(x1, y1, x2, y2, x, y *big.Rat) OrientationKind {
    95  	// Compute the determinant of the matrix
    96  	// det := x1*y2 + x2*y + x*y1 - x1*y - x2*y1 - x*y2
    97  	det := bigSub(
    98  		bigAdd(bigAdd(bigMul(x1, y2), bigMul(x2, y)), bigMul(x, y1)),
    99  		bigAdd(bigAdd(bigMul(x1, y), bigMul(x2, y1)), bigMul(x, y2)),
   100  	)
   101  	sign := det.Sign()
   102  	switch sign {
   103  	case 1:
   104  		return CCW
   105  	case -1:
   106  		return CW
   107  	case 0:
   108  		return Colinear
   109  	}
   110  	return IndeterminateOrientation
   111  }
   112  
   113  // matOrientation determines the orientation using the mat package.
   114  //
   115  // It first computes the conditional number of the matrix. When the condition number
   116  // is higher than the Condition Tolerance, then we assume the matrix is singular and
   117  // the determinant is 0. If the determinant is not 0 the sign of the determinant is computed.
   118  //
   119  //	| x1 y1 1 |
   120  //	| x2 y2 1 |
   121  //	| x  y  1 |
   122  func matOrientation(x1, y1, x2, y2, x, y float64) OrientationKind {
   123  	if (x1 == x2 && x2 == x) || (y1 == y2 && y2 == y) {
   124  		// points are horizontally or vertically aligned
   125  		return Colinear
   126  	}
   127  	m := mat.NewDense(3, 3, []float64{x1, y1, 1, x2, y2, 1, x, y, 1})
   128  	var lu mat.LU
   129  	lu.Factorize(m)
   130  	cond := lu.Cond()
   131  	if cond > mat.ConditionTolerance {
   132  		return Colinear
   133  	}
   134  	// Since only the sign is needed LogDet achieves the result in faster time.
   135  	_, sign := lu.LogDet()
   136  	switch sign {
   137  	case 1:
   138  		return CCW
   139  	case -1:
   140  		return CW
   141  	}
   142  	return IndeterminateOrientation
   143  }