EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2026-04-10 07:49:27
0001 /* 0002 * Copyright (c) 2019 Opticks Team. All Rights Reserved. 0003 * 0004 * This file is part of Opticks 0005 * (see https://bitbucket.org/simoncblyth/opticks). 0006 * 0007 * Licensed under the Apache License, Version 2.0 (the "License"); 0008 * you may not use this file except in compliance with the License. 0009 * You may obtain a copy of the License at 0010 * 0011 * http://www.apache.org/licenses/LICENSE-2.0 0012 * 0013 * Unless required by applicable law or agreed to in writing, software 0014 * distributed under the License is distributed on an "AS IS" BASIS, 0015 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 0016 * See the License for the specific language governing permissions and 0017 * limitations under the License. 0018 */ 0019 0020 #pragma once 0021 0022 0023 #include "csg.h" 0024 0025 /* 0026 0027 Robust Quadratic Roots : Numerically Stable Method for Solving Quadratic Equations 0028 ====================================================================================== 0029 0030 0031 0032 Reference 0033 ------------- 0034 0035 * https://people.csail.mit.edu/bkph/articles/Quadratics.pdf 0036 0037 Robustness in Geometric Computations 0038 Christoph M. Hoffmann 0039 Computer Science Purdue University 0040 0041 * http://www.cs.uoi.gr/~fudos/grad-material/robust4.pdf 0042 * ~/opticks_refs/robust_geometry_calc.pdf 0043 0044 0045 * :google:`numerically robust line plane intersection` 0046 0047 0048 Normal Quadratic 0049 ----------------- 0050 0051 d t^2 + 2b t + c = 0 0052 0053 0054 t = -2b +- sqrt((2b)^2 - 4dc ) 0055 ----------------------------- 0056 2d 0057 0058 -b +- sqrt( b^2 - d c ) 0059 ----------------------- 0060 d 0061 0062 0063 Alternative quadratic in 1/t 0064 -------------------------------- 0065 0066 0067 c (1/t)^2 + 2b (1/t) + d = 0 0068 0069 0070 1/t = -2b +- sqrt( (2b)^2 - 4dc ) 0071 ---------------------------- 0072 2c 0073 0074 1/t = -b +- sqrt( b^2 - d c ) 0075 ------------------------- 0076 c 0077 0078 c 0079 t = --------------------------- 0080 -b +- sqrt( b^2 - d c ) 0081 0082 0083 ---------------- 0084 0085 q = b + sign(b) sqrt( b^2 - d c ) 0086 0087 q = b + sqrt( b^2 - d c ) # b > 0 0088 q = b - sqrt( b^2 - d c ) # b < 0 0089 0090 0091 0092 */ 0093 0094 static csg_D 0095 void robust_quadratic_roots(float& t1, float &t2, float& disc, float& sdisc, const float d, const float b, const float c) 0096 { 0097 disc = b*b-d*c; 0098 sdisc = disc > 0.f ? sqrtf(disc) : 0.f ; // real roots for sdisc > 0.f 0099 0100 #ifdef NAIVE_QUADRATIC 0101 t1 = (-b - sdisc)/d ; 0102 t2 = (-b + sdisc)/d ; // t2 > t1 always, sdisc and d always +ve 0103 #else 0104 // picking robust quadratic roots that avoid catastrophic subtraction 0105 float q = b > 0.f ? -(b + sdisc) : -(b - sdisc) ; 0106 float root1 = q/d ; 0107 float root2 = c/q ; 0108 t1 = fminf( root1, root2 ); 0109 t2 = fmaxf( root1, root2 ); 0110 #endif 0111 } 0112 0113 0114 /** 0115 robust_quadratic_roots_disqualifying 0116 -------------------------------------- 0117 0118 This is a variant that avoids being forced to be used within a "disc block" 0119 by passing t_min as argument and disqualify roots by setting them 0120 to t_min when disc < 0.f 0121 0122 **/ 0123 0124 static csg_D 0125 void robust_quadratic_roots_disqualifying(const float t_min, float& t1, float &t2, float& disc, float& sdisc, const float d, const float b, const float c) 0126 { 0127 disc = b*b-d*c; 0128 sdisc = disc > 0.f ? sqrtf(disc) : 0.f ; // real roots for sdisc > 0.f 0129 float q = b > 0.f ? -(b + sdisc) : -(b - sdisc) ; // picking robust quadratic roots that avoid catastrophic subtraction 0130 float root1 = sdisc > 0.f ? q/d : t_min ; 0131 float root2 = sdisc > 0.f ? c/q : t_min ; 0132 t1 = fminf( root1, root2 ); 0133 t2 = fmaxf( root1, root2 ); 0134 } 0135 0136 0137 0138
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|