source: src/LinearAlgebra/Line.cpp@ fd4905

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable v1.0.3
Last change on this file since fd4905 was bcf653, checked in by Frederik Heber <heber@…>, 14 years ago

Added copyright note to each .cpp file and an extensive one to builder.cpp.

  • Property mode set to 100644
File size: 7.6 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010 University of Bonn. All rights reserved.
5 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
6 */
7
8/*
9 * Line.cpp
10 *
11 * Created on: Apr 30, 2010
12 * Author: crueger
13 */
14
15// include config.h
16#ifdef HAVE_CONFIG_H
17#include <config.h>
18#endif
19
20#include "Helpers/MemDebug.hpp"
21
22#include "LinearAlgebra/Line.hpp"
23
24#include <cmath>
25#include <iostream>
26
27#include "LinearAlgebra/Vector.hpp"
28#include "Helpers/Log.hpp"
29#include "Helpers/Verbose.hpp"
30#include "LinearAlgebra/gslmatrix.hpp"
31#include "Helpers/Info.hpp"
32#include "Exceptions/LinearDependenceException.hpp"
33#include "Exceptions/SkewException.hpp"
34#include "LinearAlgebra/Plane.hpp"
35
36using namespace std;
37
38Line::Line(const Vector &_origin, const Vector &_direction) :
39 direction(new Vector(_direction))
40{
41 direction->Normalize();
42 origin.reset(new Vector(_origin.partition(*direction).second));
43}
44
45Line::Line(const Line &src) :
46 origin(new Vector(*src.origin)),
47 direction(new Vector(*src.direction))
48{}
49
50Line::~Line()
51{}
52
53
54double Line::distance(const Vector &point) const{
55 // get any vector from line to point
56 Vector helper = point - *origin;
57 // partition this vector along direction
58 // the residue points from the line to the point
59 return helper.partition(*direction).second.Norm();
60}
61
62Vector Line::getClosestPoint(const Vector &point) const{
63 // get any vector from line to point
64 Vector helper = point - *origin;
65 // partition this vector along direction
66 // add only the part along the direction
67 return *origin + helper.partition(*direction).first;
68}
69
70Vector Line::getDirection() const{
71 return *direction;
72}
73
74Vector Line::getOrigin() const{
75 return *origin;
76}
77
78vector<Vector> Line::getPointsOnLine() const{
79 vector<Vector> res;
80 res.reserve(2);
81 res.push_back(*origin);
82 res.push_back(*origin+*direction);
83 return res;
84}
85
86/** Calculates the intersection of the two lines that are both on the same plane.
87 * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
88 * \param *out output stream for debugging
89 * \param *Line1a first vector of first line
90 * \param *Line1b second vector of first line
91 * \param *Line2a first vector of second line
92 * \param *Line2b second vector of second line
93 * \return true - \a this will contain the intersection on return, false - lines are parallel
94 */
95Vector Line::getIntersection(const Line& otherLine) const{
96 Info FunctionInfo(__func__);
97
98 pointset line1Points = getPointsOnLine();
99
100 Vector Line1a = line1Points[0];
101 Vector Line1b = line1Points[1];
102
103 pointset line2Points = otherLine.getPointsOnLine();
104
105 Vector Line2a = line2Points[0];
106 Vector Line2b = line2Points[1];
107
108 Vector res;
109
110 auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
111
112 M->SetAll(1.);
113 for (int i=0;i<3;i++) {
114 M->Set(0, i, Line1a[i]);
115 M->Set(1, i, Line1b[i]);
116 M->Set(2, i, Line2a[i]);
117 M->Set(3, i, Line2b[i]);
118 }
119
120 //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
121 //for (int i=0;i<4;i++) {
122 // for (int j=0;j<4;j++)
123 // cout << "\t" << M->Get(i,j);
124 // cout << endl;
125 //}
126 if (fabs(M->Determinant()) > MYEPSILON) {
127 Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
128 throw SkewException(__FILE__,__LINE__);
129 }
130
131 Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
132
133
134 // constuct a,b,c
135 Vector a = Line1b - Line1a;
136 Vector b = Line2b - Line2a;
137 Vector c = Line2a - Line1a;
138 Vector d = Line2b - Line1b;
139 Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
140 if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
141 res.Zero();
142 Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
143 throw LinearDependenceException(__FILE__,__LINE__);
144 }
145
146 // check for parallelity
147 Vector parallel;
148 double factor = 0.;
149 if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
150 parallel = Line1a - Line2a;
151 factor = parallel.ScalarProduct(a)/a.Norm();
152 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
153 res = Line2a;
154 Log() << Verbose(1) << "Lines conincide." << endl;
155 return res;
156 } else {
157 parallel = Line1a - Line2b;
158 factor = parallel.ScalarProduct(a)/a.Norm();
159 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
160 res = Line2b;
161 Log() << Verbose(1) << "Lines conincide." << endl;
162 return res;
163 }
164 }
165 Log() << Verbose(1) << "Lines are parallel." << endl;
166 res.Zero();
167 throw LinearDependenceException(__FILE__,__LINE__);
168 }
169
170 // obtain s
171 double s;
172 Vector temp1, temp2;
173 temp1 = c;
174 temp1.VectorProduct(b);
175 temp2 = a;
176 temp2.VectorProduct(b);
177 Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
178 if (fabs(temp2.NormSquared()) > MYEPSILON)
179 s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
180 else
181 s = 0.;
182 Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
183
184 // construct intersection
185 res = a;
186 res.Scale(s);
187 res += Line1a;
188 Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
189
190 return res;
191}
192
193/** Rotates the vector by an angle of \a alpha around this line.
194 * \param rhs Vector to rotate
195 * \param alpha rotation angle in radian
196 */
197Vector Line::rotateVector(const Vector &rhs, double alpha) const{
198 Vector helper = rhs;
199
200 // translate the coordinate system so that the line goes through (0,0,0)
201 helper -= *origin;
202
203 // partition the vector into a part that gets rotated and a part that lies along the line
204 pair<Vector,Vector> parts = helper.partition(*direction);
205
206 // we just keep anything that is along the axis
207 Vector res = parts.first;
208
209 // the rest has to be rotated
210 Vector a = parts.second;
211 // we only have to do the rest, if we actually could partition the vector
212 if(!a.IsZero()){
213 // construct a vector that is orthogonal to a and direction and has length |a|
214 Vector y = a;
215 // direction is normalized, so the result has length |a|
216 y.VectorProduct(*direction);
217
218 res += cos(alpha) * a + sin(alpha) * y;
219 }
220
221 // translate the coordinate system back
222 res += *origin;
223 return res;
224}
225
226Plane Line::getOrthogonalPlane(const Vector &origin) const{
227 return Plane(getDirection(),origin);
228}
229
230std::vector<Vector> Line::getSphereIntersections() const{
231 std::vector<Vector> res;
232
233 // line is kept in normalized form, so we can skip a lot of calculations
234 double discriminant = 1-origin->NormSquared();
235 // we might have 2, 1 or 0 solutions, depending on discriminant
236 if(discriminant>=0){
237 if(discriminant==0){
238 res.push_back(*origin);
239 }
240 else{
241 Vector helper = sqrt(discriminant)*(*direction);
242 res.push_back(*origin+helper);
243 res.push_back(*origin-helper);
244 }
245 }
246 return res;
247}
248
249Line makeLineThrough(const Vector &x1, const Vector &x2){
250 if(x1==x2){
251 throw LinearDependenceException(__FILE__,__LINE__);
252 }
253 return Line(x1,x1-x2);
254}
255
256ostream& operator<<(ostream& ost, const Line& m)
257{
258 const Vector origin = m.getOrigin();
259 const Vector direction = m.getDirection();
260 ost << "(";
261 for (int i=0;i<NDIM;i++) {
262 ost << origin[i];
263 if (i != 2)
264 ost << ",";
265 }
266 ost << ") -> (";
267 for (int i=0;i<NDIM;i++) {
268 ost << direction[i];
269 if (i != 2)
270 ost << ",";
271 }
272 ost << ")";
273 return ost;
274};
275
Note: See TracBrowser for help on using the repository browser.