Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, * 20 // * subject to DEFCON 705 IPR conditions. 20 // * subject to DEFCON 705 IPR conditions. * 21 // * By using, copying, modifying or distri 21 // * By using, copying, modifying or distributing the software (or * 22 // * any work based on the software) you ag 22 // * any work based on the software) you agree to acknowledge its * 23 // * use in resulting scientific publicati 23 // * use in resulting scientific publications, and indicate your * 24 // * acceptance of all terms of the Geant4 Sof 24 // * acceptance of all terms of the Geant4 Software license. * 25 // ******************************************* 25 // ******************************************************************** 26 // 26 // 27 // G4TessellatedGeometryAlgorithms implementat << 27 // >> 28 // >> 29 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> 30 // >> 31 // CHANGE HISTORY >> 32 // -------------- 28 // 33 // 29 // 07 August 2007, P R Truscott, QinetiQ Ltd, 34 // 07 August 2007, P R Truscott, QinetiQ Ltd, UK - Created, with member 30 // functions based on the work 35 // functions based on the work of Rickard Holmberg. >> 36 // >> 37 // 26 September 2007 >> 38 // P R Truscott, qinetiQ Ltd, UK >> 39 // Updated to assign values of location array, not update >> 40 // just the pointer. >> 41 // 31 // 12 October 2012, M Gayer, CERN, - Reviewed 42 // 12 October 2012, M Gayer, CERN, - Reviewed optimized implementation. 32 // ------------------------------------------- << 43 // >> 44 /////////////////////////////////////////////////////////////////////////////// 33 45 34 #include "G4TessellatedGeometryAlgorithms.hh" 46 #include "G4TessellatedGeometryAlgorithms.hh" 35 47 36 #include <cfloat> 48 #include <cfloat> 37 49 38 ////////////////////////////////////////////// 50 /////////////////////////////////////////////////////////////////////////////// 39 // 51 // 40 // IntersectLineAndTriangle2D 52 // IntersectLineAndTriangle2D 41 // 53 // 42 // Determines whether there is an intersection 54 // Determines whether there is an intersection between a line defined 43 // by r = p + s.v and a triangle defined by ve 55 // by r = p + s.v and a triangle defined by verticies p0, p0+e0 and p0+e1. 44 // 56 // 45 // Here: 57 // Here: 46 // p = 2D vector 58 // p = 2D vector 47 // s = scaler on [0,infinity) 59 // s = scaler on [0,infinity) 48 // v = 2D vector 60 // v = 2D vector 49 // p0, e0 and e1 are 2D vectors 61 // p0, e0 and e1 are 2D vectors 50 // Information about where the intersection oc 62 // Information about where the intersection occurs is returned in the 51 // variable location. 63 // variable location. 52 // 64 // 53 // This is based on the work of Rickard Holmbe 65 // This is based on the work of Rickard Holmberg. 54 // 66 // 55 G4bool G4TessellatedGeometryAlgorithms::Inters 67 G4bool G4TessellatedGeometryAlgorithms::IntersectLineAndTriangle2D ( 56 const G4TwoVector& p, const G4TwoVector& v, << 68 const G4TwoVector &p, const G4TwoVector &v, 57 const G4TwoVector& p0, const G4TwoVector& e0 << 69 const G4TwoVector &p0, const G4TwoVector &e0, const G4TwoVector &e1, 58 G4TwoVector location[2]) 70 G4TwoVector location[2]) 59 { 71 { 60 G4TwoVector loc0[2]; 72 G4TwoVector loc0[2]; 61 G4int e0i = IntersectLineAndLineSegment2D (p 73 G4int e0i = IntersectLineAndLineSegment2D (p,v,p0,e0,loc0); 62 if (e0i == 2) 74 if (e0i == 2) 63 { 75 { 64 location[0] = loc0[0]; 76 location[0] = loc0[0]; 65 location[1] = loc0[1]; 77 location[1] = loc0[1]; 66 return true; 78 return true; 67 } 79 } 68 80 69 G4TwoVector loc1[2]; 81 G4TwoVector loc1[2]; 70 G4int e1i = IntersectLineAndLineSegment2D (p 82 G4int e1i = IntersectLineAndLineSegment2D (p,v,p0,e1,loc1); 71 if (e1i == 2) 83 if (e1i == 2) 72 { 84 { 73 location[0] = loc1[0]; 85 location[0] = loc1[0]; 74 location[1] = loc1[1]; 86 location[1] = loc1[1]; 75 return true; 87 return true; 76 } 88 } 77 89 78 if ((e0i == 1) && (e1i == 1)) 90 if ((e0i == 1) && (e1i == 1)) 79 { 91 { 80 if ((loc0[0]-p).mag2() < (loc1[0]-p).mag2( 92 if ((loc0[0]-p).mag2() < (loc1[0]-p).mag2()) 81 { 93 { 82 location[0] = loc0[0]; 94 location[0] = loc0[0]; 83 location[1] = loc1[0]; 95 location[1] = loc1[0]; 84 } 96 } 85 else 97 else 86 { 98 { 87 location[0] = loc1[0]; 99 location[0] = loc1[0]; 88 location[1] = loc0[0]; 100 location[1] = loc0[0]; 89 } 101 } 90 return true; 102 return true; 91 } 103 } 92 104 93 G4TwoVector p1 = p0 + e0; 105 G4TwoVector p1 = p0 + e0; 94 G4TwoVector DE = e1 - e0; 106 G4TwoVector DE = e1 - e0; 95 G4TwoVector loc2[2]; 107 G4TwoVector loc2[2]; 96 G4int e2i = IntersectLineAndLineSegment2D (p 108 G4int e2i = IntersectLineAndLineSegment2D (p,v,p1,DE,loc2); 97 if (e2i == 2) 109 if (e2i == 2) 98 { 110 { 99 location[0] = loc2[0]; 111 location[0] = loc2[0]; 100 location[1] = loc2[1]; 112 location[1] = loc2[1]; 101 return true; 113 return true; 102 } 114 } 103 115 104 if ((e0i == 0) && (e1i == 0) && (e2i == 0)) 116 if ((e0i == 0) && (e1i == 0) && (e2i == 0)) return false; 105 117 106 if ((e0i == 1) && (e2i == 1)) 118 if ((e0i == 1) && (e2i == 1)) 107 { 119 { 108 if ((loc0[0]-p).mag2() < (loc2[0]-p).mag2( 120 if ((loc0[0]-p).mag2() < (loc2[0]-p).mag2()) 109 { 121 { 110 location[0] = loc0[0]; 122 location[0] = loc0[0]; 111 location[1] = loc2[0]; 123 location[1] = loc2[0]; 112 } 124 } 113 else 125 else 114 { 126 { 115 location[0] = loc2[0]; 127 location[0] = loc2[0]; 116 location[1] = loc0[0]; 128 location[1] = loc0[0]; 117 } 129 } 118 return true; 130 return true; 119 } 131 } 120 132 121 if ((e1i == 1) && (e2i == 1)) 133 if ((e1i == 1) && (e2i == 1)) 122 { 134 { 123 if ((loc1[0]-p).mag2() < (loc2[0]-p).mag2( 135 if ((loc1[0]-p).mag2() < (loc2[0]-p).mag2()) 124 { 136 { 125 location[0] = loc1[0]; 137 location[0] = loc1[0]; 126 location[1] = loc2[0]; 138 location[1] = loc2[0]; 127 } 139 } 128 else 140 else 129 { 141 { 130 location[0] = loc2[0]; 142 location[0] = loc2[0]; 131 location[1] = loc1[0]; 143 location[1] = loc1[0]; 132 } 144 } 133 return true; 145 return true; 134 } 146 } 135 147 136 return false; 148 return false; 137 } 149 } 138 150 139 ////////////////////////////////////////////// 151 /////////////////////////////////////////////////////////////////////////////// 140 // 152 // 141 // IntersectLineAndLineSegment2D 153 // IntersectLineAndLineSegment2D 142 // 154 // 143 // Determines whether there is an intersection 155 // Determines whether there is an intersection between a line defined 144 // by r = p0 + s.d0 and a line-segment with en 156 // by r = p0 + s.d0 and a line-segment with endpoints p1 and p1+d1. 145 // Here: 157 // Here: 146 // p0 = 2D vector 158 // p0 = 2D vector 147 // s = scaler on [0,infinity) 159 // s = scaler on [0,infinity) 148 // d0 = 2D vector 160 // d0 = 2D vector 149 // p1 and d1 are 2D vectors 161 // p1 and d1 are 2D vectors 150 // 162 // 151 // This function returns: 163 // This function returns: 152 // 0 - if there is no intersection; 164 // 0 - if there is no intersection; 153 // 1 - if there is a unique intersection; 165 // 1 - if there is a unique intersection; 154 // 2 - if the line and line-segments overlap, 166 // 2 - if the line and line-segments overlap, and the intersection is a 155 // segment itself. 167 // segment itself. 156 // Information about where the intersection oc 168 // Information about where the intersection occurs is returned in the 157 // as ??. 169 // as ??. 158 // 170 // 159 // This is based on the work of Rickard Holmbe 171 // This is based on the work of Rickard Holmberg as well as material published 160 // by Philip J Schneider and David H Eberly, " 172 // by Philip J Schneider and David H Eberly, "Geometric Tools for Computer 161 // Graphics," ISBN 1-55860-694-0, pp 244-245, 173 // Graphics," ISBN 1-55860-694-0, pp 244-245, 2003. 162 // 174 // 163 G4int G4TessellatedGeometryAlgorithms::Interse 175 G4int G4TessellatedGeometryAlgorithms::IntersectLineAndLineSegment2D ( 164 const G4TwoVector& p0, const G4TwoVector& d0 << 176 const G4TwoVector &p0, const G4TwoVector &d0, 165 const G4TwoVector& p1, const G4TwoVector& d1 << 177 const G4TwoVector &p1, const G4TwoVector &d1, >> 178 G4TwoVector location[2]) 166 { 179 { 167 G4TwoVector e = p1 - p0; 180 G4TwoVector e = p1 - p0; 168 G4double kross = cross(d0,d1); 181 G4double kross = cross(d0,d1); 169 G4double sqrKross = kross * kross; 182 G4double sqrKross = kross * kross; 170 G4double sqrLen0 = d0.mag2(); 183 G4double sqrLen0 = d0.mag2(); 171 G4double sqrLen1 = d1.mag2(); 184 G4double sqrLen1 = d1.mag2(); 172 location[0] = G4TwoVector(0.0,0.0); 185 location[0] = G4TwoVector(0.0,0.0); 173 location[1] = G4TwoVector(0.0,0.0); 186 location[1] = G4TwoVector(0.0,0.0); 174 187 175 if (sqrKross > DBL_EPSILON * DBL_EPSILON * s 188 if (sqrKross > DBL_EPSILON * DBL_EPSILON * sqrLen0 * sqrLen1) 176 { 189 { 177 // 190 // 178 // The line and line segment are not paral 191 // The line and line segment are not parallel. Determine if the intersection 179 // is in positive s where r=p0 + s*d0, and 192 // is in positive s where r=p0 + s*d0, and for 0<=t<=1 where r=p1 + t*d1. 180 // 193 // 181 G4double ss = cross(e,d1)/kross; 194 G4double ss = cross(e,d1)/kross; 182 if (ss < 0) return 0; // Intersect 195 if (ss < 0) return 0; // Intersection does not occur for positive ss 183 G4double t = cross(e,d0)/kross; 196 G4double t = cross(e,d0)/kross; 184 if (t < 0 || t > 1) return 0; // Intersect 197 if (t < 0 || t > 1) return 0; // Intersection does not occur on line-segment 185 // 198 // 186 // Intersection of lines is a single point 199 // Intersection of lines is a single point on the forward-propagating line 187 // defined by r=p0 + ss*d0, and the line s 200 // defined by r=p0 + ss*d0, and the line segment defined by r=p1 + t*d1. 188 // 201 // 189 location[0] = p0 + ss*d0; 202 location[0] = p0 + ss*d0; 190 return 1; 203 return 1; 191 } 204 } 192 // 205 // 193 // Line and line segment are parallel. Deter 206 // Line and line segment are parallel. Determine whether they overlap or not. 194 // 207 // 195 G4double sqrLenE = e.mag2(); 208 G4double sqrLenE = e.mag2(); 196 kross = cross(e,d0); 209 kross = cross(e,d0); 197 sqrKross = kross * kross; 210 sqrKross = kross * kross; 198 if (sqrKross > DBL_EPSILON * DBL_EPSILON * s 211 if (sqrKross > DBL_EPSILON * DBL_EPSILON * sqrLen0 * sqrLenE) 199 { 212 { 200 return 0; //Lines are different. 213 return 0; //Lines are different. 201 } 214 } 202 // 215 // 203 // Lines are the same. Test for overlap. 216 // Lines are the same. Test for overlap. 204 // 217 // 205 G4double s0 = d0.dot(e)/sqrLen0; 218 G4double s0 = d0.dot(e)/sqrLen0; 206 G4double s1 = s0 + d0.dot(d1)/sqrLen0; 219 G4double s1 = s0 + d0.dot(d1)/sqrLen0; 207 G4double smin = 0.0; 220 G4double smin = 0.0; 208 G4double smax = 0.0; 221 G4double smax = 0.0; 209 222 210 if (s0 < s1) {smin = s0; smax = s1;} 223 if (s0 < s1) {smin = s0; smax = s1;} 211 else {smin = s1; smax = s0;} 224 else {smin = s1; smax = s0;} 212 225 213 if (smax < 0.0) return 0; 226 if (smax < 0.0) return 0; 214 else if (smin < 0.0) 227 else if (smin < 0.0) 215 { 228 { 216 location[0] = p0; 229 location[0] = p0; 217 location[1] = p0 + smax*d0; 230 location[1] = p0 + smax*d0; 218 return 2; 231 return 2; 219 } 232 } 220 else 233 else 221 { 234 { 222 location[0] = p0 + smin*d0; 235 location[0] = p0 + smin*d0; 223 location[1] = p0 + smax*d0; 236 location[1] = p0 + smax*d0; 224 return 2; 237 return 2; 225 } 238 } 226 } 239 } 227 240 228 ////////////////////////////////////////////// 241 /////////////////////////////////////////////////////////////////////////////// 229 // 242 // 230 // CrossProduct 243 // CrossProduct 231 // 244 // 232 // This is just a ficticious "cross-product" f 245 // This is just a ficticious "cross-product" function for two 2D vectors... 233 // "ficticious" because such an operation is n 246 // "ficticious" because such an operation is not relevant to 2D space compared 234 // with 3D space. 247 // with 3D space. 235 // 248 // 236 G4double G4TessellatedGeometryAlgorithms::cros << 249 G4double G4TessellatedGeometryAlgorithms::cross(const G4TwoVector &v1, 237 << 250 const G4TwoVector &v2) 238 { 251 { 239 return v1.x()*v2.y() - v1.y()*v2.x(); 252 return v1.x()*v2.y() - v1.y()*v2.x(); 240 } 253 } 241 254