source: src/ellipsoid.cpp@ c38826

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
Last change on this file since c38826 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: 17.4 KB
RevLine 
[bcf653]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
[6ac7ee]8/*
9 * ellipsoid.cpp
10 *
[042f82]11 * Created on: Jan 20, 2009
12 * Author: heber
[6ac7ee]13 */
14
[bf3817]15// include config.h
16#ifdef HAVE_CONFIG_H
17#include <config.h>
18#endif
19
[112b09]20#include "Helpers/MemDebug.hpp"
21
[357fba]22#include <gsl/gsl_multimin.h>
23#include <gsl/gsl_vector.h>
24
[f66195]25#include <iomanip>
26
27#include <set>
28
[d74077]29#include "BoundaryPointSet.hpp"
[357fba]30#include "boundary.hpp"
[6ac7ee]31#include "ellipsoid.hpp"
[f66195]32#include "linkedcell.hpp"
[952f38]33#include "Helpers/Log.hpp"
[f66195]34#include "tesselation.hpp"
[57f243]35#include "LinearAlgebra/Vector.hpp"
36#include "LinearAlgebra/Matrix.hpp"
[952f38]37#include "Helpers/Verbose.hpp"
[6ac7ee]38
39/** Determines squared distance for a given point \a x to surface of ellipsoid.
40 * \param x given point
41 * \param EllipsoidCenter center of ellipsoid
42 * \param EllipsoidLength[3] three lengths of half axis of ellipsoid
43 * \param EllipsoidAngle[3] three rotation angles of ellipsoid
44 * \return squared distance from point to surface
45 */
46double SquaredDistanceToEllipsoid(Vector &x, Vector &EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle)
47{
[042f82]48 Vector helper, RefPoint;
49 double distance = -1.;
[c94eeb]50 Matrix Matrix;
[042f82]51 double InverseLength[3];
52 double psi,theta,phi; // euler angles in ZX'Z'' convention
53
[e138de]54 //Log() << Verbose(3) << "Begin of SquaredDistanceToEllipsoid" << endl;
[042f82]55
56 for(int i=0;i<3;i++)
57 InverseLength[i] = 1./EllipsoidLength[i];
58
59 // 1. translate coordinate system so that ellipsoid center is in origin
[273382]60 RefPoint = helper = x - EllipsoidCenter;
[e138de]61 //Log() << Verbose(4) << "Translated given point is at " << RefPoint << "." << endl;
[042f82]62
63 // 2. transform coordinate system by inverse of rotation matrix and of diagonal matrix
64 psi = EllipsoidAngle[0];
65 theta = EllipsoidAngle[1];
66 phi = EllipsoidAngle[2];
[a679d1]67 Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
68 Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
69 Matrix.set(2,0, sin(psi)*sin(theta));
70 Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
71 Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
72 Matrix.set(2,1, -cos(psi)*sin(theta));
73 Matrix.set(0,2, sin(theta)*sin(phi));
74 Matrix.set(1,2, sin(theta)*cos(phi));
75 Matrix.set(2,2, cos(theta));
[5108e1]76 helper *= Matrix;
[1bd79e]77 helper.ScaleAll(InverseLength);
[e138de]78 //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
[042f82]79
80 // 3. construct intersection point with unit sphere and ray between origin and x
81 helper.Normalize(); // is simply normalizes vector in distance direction
[e138de]82 //Log() << Verbose(4) << "Transformed intersection is at " << helper << "." << endl;
[042f82]83
84 // 4. transform back the constructed intersection point
85 psi = -EllipsoidAngle[0];
86 theta = -EllipsoidAngle[1];
87 phi = -EllipsoidAngle[2];
[1bd79e]88 helper.ScaleAll(EllipsoidLength);
[a679d1]89 Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
90 Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
91 Matrix.set(2,0, sin(psi)*sin(theta));
92 Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
93 Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
94 Matrix.set(2,1, -cos(psi)*sin(theta));
95 Matrix.set(0,2, sin(theta)*sin(phi));
96 Matrix.set(1,2, sin(theta)*cos(phi));
97 Matrix.set(2,2, cos(theta));
[5108e1]98 helper *= Matrix;
[e138de]99 //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;
[042f82]100
101 // 5. determine distance between backtransformed point and x
[273382]102 distance = RefPoint.DistanceSquared(helper);
[e138de]103 //Log() << Verbose(4) << "Squared distance between intersection and RefPoint is " << distance << "." << endl;
[042f82]104
105 return distance;
[e138de]106 //Log() << Verbose(3) << "End of SquaredDistanceToEllipsoid" << endl;
[6ac7ee]107};
108
109/** structure for ellipsoid minimisation containing points to fit to.
110 */
111struct EllipsoidMinimisation {
[042f82]112 int N; //!< dimension of vector set
113 Vector *x; //!< array of vectors
[6ac7ee]114};
115
116/** Sum of squared distance to ellipsoid to be minimised.
117 * \param *x parameters for the ellipsoid
118 * \param *params EllipsoidMinimisation with set of data points to minimise distance to and dimension
119 * \return sum of squared distance, \sa SquaredDistanceToEllipsoid()
120 */
121double SumSquaredDistance (const gsl_vector * x, void * params)
122{
[042f82]123 Vector *set= ((struct EllipsoidMinimisation *)params)->x;
124 int N = ((struct EllipsoidMinimisation *)params)->N;
125 double SumDistance = 0.;
126 double distance;
127 Vector Center;
128 double EllipsoidLength[3], EllipsoidAngle[3];
129
130 // put parameters into suitable ellipsoid form
131 for (int i=0;i<3;i++) {
[0a4f7f]132 Center[i] = gsl_vector_get(x, i+0);
[042f82]133 EllipsoidLength[i] = gsl_vector_get(x, i+3);
134 EllipsoidAngle[i] = gsl_vector_get(x, i+6);
135 }
136
137 // go through all points and sum distance
138 for (int i=0;i<N;i++) {
139 distance = SquaredDistanceToEllipsoid(set[i], Center, EllipsoidLength, EllipsoidAngle);
140 if (!isnan(distance)) {
141 SumDistance += distance;
142 } else {
143 SumDistance = GSL_NAN;
144 break;
145 }
146 }
147
[e138de]148 //Log() << Verbose(0) << "Current summed distance is " << SumDistance << "." << endl;
[042f82]149 return SumDistance;
[6ac7ee]150};
151
152/** Finds best fitting ellipsoid parameter set in Least square sense for a given point set.
153 * \param *out output stream for debugging
154 * \param *set given point set
155 * \param N number of points in set
156 * \param EllipsoidParamter[3] three parameters in ellipsoid equation
157 * \return true - fit successful, false - fit impossible
158 */
[e138de]159bool FitPointSetToEllipsoid(Vector *set, int N, Vector *EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle)
[6ac7ee]160{
[042f82]161 int status = GSL_SUCCESS;
[a67d19]162 DoLog(2) && (Log() << Verbose(2) << "Begin of FitPointSetToEllipsoid " << endl);
[042f82]163 if (N >= 3) { // check that enough points are given (9 d.o.f.)
164 struct EllipsoidMinimisation par;
165 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
166 gsl_multimin_fminimizer *s = NULL;
167 gsl_vector *ss, *x;
168 gsl_multimin_function minex_func;
169
170 size_t iter = 0;
171 double size;
172
173 /* Starting point */
174 x = gsl_vector_alloc (9);
175 for (int i=0;i<3;i++) {
[0a4f7f]176 gsl_vector_set (x, i+0, EllipsoidCenter->at(i));
[042f82]177 gsl_vector_set (x, i+3, EllipsoidLength[i]);
178 gsl_vector_set (x, i+6, EllipsoidAngle[i]);
179 }
180 par.x = set;
181 par.N = N;
182
183 /* Set initial step sizes */
184 ss = gsl_vector_alloc (9);
185 for (int i=0;i<3;i++) {
186 gsl_vector_set (ss, i+0, 0.1);
187 gsl_vector_set (ss, i+3, 1.0);
188 gsl_vector_set (ss, i+6, M_PI/20.);
189 }
190
191 /* Initialize method and iterate */
192 minex_func.n = 9;
193 minex_func.f = &SumSquaredDistance;
194 minex_func.params = (void *)&par;
195
196 s = gsl_multimin_fminimizer_alloc (T, 9);
197 gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
198
199 do {
200 iter++;
201 status = gsl_multimin_fminimizer_iterate(s);
202
203 if (status)
204 break;
205
206 size = gsl_multimin_fminimizer_size (s);
207 status = gsl_multimin_test_size (size, 1e-2);
208
209 if (status == GSL_SUCCESS) {
210 for (int i=0;i<3;i++) {
[0a4f7f]211 EllipsoidCenter->at(i) = gsl_vector_get (s->x,i+0);
[042f82]212 EllipsoidLength[i] = gsl_vector_get (s->x, i+3);
213 EllipsoidAngle[i] = gsl_vector_get (s->x, i+6);
214 }
[a67d19]215 DoLog(4) && (Log() << Verbose(4) << setprecision(3) << "Converged fit at: " << *EllipsoidCenter << ", lengths " << EllipsoidLength[0] << ", " << EllipsoidLength[1] << ", " << EllipsoidLength[2] << ", angles " << EllipsoidAngle[0] << ", " << EllipsoidAngle[1] << ", " << EllipsoidAngle[2] << " with summed distance " << s->fval << "." << endl);
[042f82]216 }
217
218 } while (status == GSL_CONTINUE && iter < 1000);
219
220 gsl_vector_free(x);
221 gsl_vector_free(ss);
222 gsl_multimin_fminimizer_free (s);
223
224 } else {
[a67d19]225 DoLog(3) && (Log() << Verbose(3) << "Not enough points provided for fit to ellipsoid." << endl);
[042f82]226 return false;
227 }
[a67d19]228 DoLog(2) && (Log() << Verbose(2) << "End of FitPointSetToEllipsoid" << endl);
[042f82]229 if (status == GSL_SUCCESS)
230 return true;
231 else
232 return false;
[6ac7ee]233};
234
235/** Picks a number of random points from a LC neighbourhood as a fitting set.
236 * \param *out output stream for debugging
237 * \param *T Tesselation containing boundary points
238 * \param *LC linked cell list of all atoms
239 * \param *&x random point set on return (not allocated!)
240 * \param PointsToPick number of points in set to pick
241 */
[e138de]242void PickRandomNeighbouredPointSet(class Tesselation *T, class LinkedCell *LC, Vector *&x, size_t PointsToPick)
[6ac7ee]243{
[70c333f]244 size_t PointsLeft = 0;
245 size_t PointsPicked = 0;
[042f82]246 int Nlower[NDIM], Nupper[NDIM];
247 set<int> PickedAtomNrs; // ordered list of picked atoms
248 set<int>::iterator current;
249 int index;
[357fba]250 TesselPoint *Candidate = NULL;
[a67d19]251 DoLog(2) && (Log() << Verbose(2) << "Begin of PickRandomPointSet" << endl);
[042f82]252
253 // allocate array
254 if (x == NULL) {
255 x = new Vector[PointsToPick];
256 } else {
[58ed4a]257 DoeLog(2) && (eLog()<< Verbose(2) << "Given pointer to vector array seems already allocated." << endl);
[042f82]258 }
259
260 do {
261 for(int i=0;i<NDIM;i++) // pick three random indices
262 LC->n[i] = (rand() % LC->N[i]);
[a67d19]263 DoLog(2) && (Log() << Verbose(2) << "INFO: Center cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " ... ");
[042f82]264 // get random cell
[734816]265 const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
[042f82]266 if (List == NULL) { // set index to it
267 continue;
268 }
[a67d19]269 DoLog(2) && (Log() << Verbose(2) << "with No. " << LC->index << "." << endl);
[042f82]270
[a67d19]271 DoLog(2) && (Log() << Verbose(2) << "LC Intervals:");
[042f82]272 for (int i=0;i<NDIM;i++) {
273 Nlower[i] = ((LC->n[i]-1) >= 0) ? LC->n[i]-1 : 0;
274 Nupper[i] = ((LC->n[i]+1) < LC->N[i]) ? LC->n[i]+1 : LC->N[i]-1;
[a67d19]275 DoLog(0) && (Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ");
[042f82]276 }
[a67d19]277 DoLog(0) && (Log() << Verbose(0) << endl);
[042f82]278
279 // count whether there are sufficient atoms in this cell+neighbors
280 PointsLeft=0;
281 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
282 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
283 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
[734816]284 const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
[042f82]285 PointsLeft += List->size();
286 }
[a67d19]287 DoLog(2) && (Log() << Verbose(2) << "There are " << PointsLeft << " atoms in this neighbourhood." << endl);
[042f82]288 if (PointsLeft < PointsToPick) { // ensure that we can pick enough points in its neighbourhood at all.
289 continue;
290 }
291
292 // pre-pick a fixed number of atoms
293 PickedAtomNrs.clear();
294 do {
295 index = (rand() % PointsLeft);
296 current = PickedAtomNrs.find(index); // not present?
297 if (current == PickedAtomNrs.end()) {
[e138de]298 //Log() << Verbose(2) << "Picking atom nr. " << index << "." << endl;
[042f82]299 PickedAtomNrs.insert(index);
300 }
301 } while (PickedAtomNrs.size() < PointsToPick);
302
303 index = 0; // now go through all and pick those whose from PickedAtomsNr
304 PointsPicked=0;
305 current = PickedAtomNrs.begin();
306 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
307 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
308 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
[734816]309 const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
[e138de]310// Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
[042f82]311 if (List != NULL) {
312// if (List->begin() != List->end())
[e138de]313// Log() << Verbose(2) << "Going through candidates ... " << endl;
[042f82]314// else
[e138de]315// Log() << Verbose(2) << "Cell is empty ... " << endl;
[734816]316 for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
[042f82]317 if ((current != PickedAtomNrs.end()) && (*current == index)) {
318 Candidate = (*Runner);
[d74077]319 DoLog(2) && (Log() << Verbose(2) << "Current picked node is " << (*Runner)->getName() << " with index " << index << "." << endl);
320 x[PointsPicked++] = Candidate->getPosition(); // we have one more atom picked
[042f82]321 current++; // next pre-picked atom
322 }
323 index++; // next atom nr.
324 }
325// } else {
[e138de]326// Log() << Verbose(2) << "List for this index not allocated!" << endl;
[042f82]327 }
328 }
[a67d19]329 DoLog(2) && (Log() << Verbose(2) << "The following points were picked: " << endl);
[042f82]330 for (size_t i=0;i<PointsPicked;i++)
[a67d19]331 DoLog(2) && (Log() << Verbose(2) << x[i] << endl);
[042f82]332 if (PointsPicked == PointsToPick) // break out of loop if we have all
333 break;
334 } while(1);
335
[a67d19]336 DoLog(2) && (Log() << Verbose(2) << "End of PickRandomPointSet" << endl);
[6ac7ee]337};
338
339/** Picks a number of random points from a set of boundary points as a fitting set.
340 * \param *out output stream for debugging
341 * \param *T Tesselation containing boundary points
342 * \param *&x random point set on return (not allocated!)
343 * \param PointsToPick number of points in set to pick
344 */
[e138de]345void PickRandomPointSet(class Tesselation *T, Vector *&x, size_t PointsToPick)
[6ac7ee]346{
[70c333f]347 size_t PointsLeft = (size_t) T->PointsOnBoundaryCount;
348 size_t PointsPicked = 0;
[042f82]349 double value, threshold;
350 PointMap *List = &T->PointsOnBoundary;
[a67d19]351 DoLog(2) && (Log() << Verbose(2) << "Begin of PickRandomPointSet" << endl);
[042f82]352
353 // allocate array
354 if (x == NULL) {
355 x = new Vector[PointsToPick];
356 } else {
[58ed4a]357 DoeLog(2) && (eLog()<< Verbose(2) << "Given pointer to vector array seems already allocated." << endl);
[042f82]358 }
359
360 if (List != NULL)
361 for (PointMap::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
362 threshold = 1. - (double)(PointsToPick - PointsPicked)/(double)PointsLeft;
363 value = (double)rand()/(double)RAND_MAX;
[e138de]364 //Log() << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": ";
[042f82]365 if (value > threshold) {
[d74077]366 x[PointsPicked] = (Runner->second->node->getPosition());
[042f82]367 PointsPicked++;
[e138de]368 //Log() << Verbose(0) << "IN." << endl;
[042f82]369 } else {
[e138de]370 //Log() << Verbose(0) << "OUT." << endl;
[042f82]371 }
372 PointsLeft--;
373 }
[a67d19]374 DoLog(2) && (Log() << Verbose(2) << "The following points were picked: " << endl);
[042f82]375 for (size_t i=0;i<PointsPicked;i++)
[a67d19]376 DoLog(3) && (Log() << Verbose(3) << x[i] << endl);
[042f82]377
[a67d19]378 DoLog(2) && (Log() << Verbose(2) << "End of PickRandomPointSet" << endl);
[6ac7ee]379};
380
381/** Finds best fitting ellipsoid parameter set in least square sense for a given point set.
382 * \param *out output stream for debugging
383 * \param *T Tesselation containing boundary points
384 * \param *LCList linked cell list of all atoms
385 * \param N number of unique points in ellipsoid fit, must be greater equal 6
386 * \param number of fits (i.e. parameter sets in output file)
387 * \param *filename name for output file
388 */
[e138de]389void FindDistributionOfEllipsoids(class Tesselation *T, class LinkedCell *LCList, int N, int number, const char *filename)
[6ac7ee]390{
[042f82]391 ofstream output;
392 Vector *x = NULL;
393 Vector Center;
394 Vector EllipsoidCenter;
395 double EllipsoidLength[3];
396 double EllipsoidAngle[3];
397 double distance, MaxDistance, MinDistance;
[a67d19]398 DoLog(0) && (Log() << Verbose(0) << "Begin of FindDistributionOfEllipsoids" << endl);
[042f82]399
400 // construct center of gravity of boundary point set for initial ellipsoid center
401 Center.Zero();
402 for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++)
[d74077]403 Center += (Runner->second->node->getPosition());
[042f82]404 Center.Scale(1./T->PointsOnBoundaryCount);
[a67d19]405 DoLog(1) && (Log() << Verbose(1) << "Center is at " << Center << "." << endl);
[042f82]406
407 // Output header
408 output.open(filename, ios::trunc);
409 output << "# Nr.\tCenterX\tCenterY\tCenterZ\ta\tb\tc\tpsi\ttheta\tphi" << endl;
410
411 // loop over desired number of parameter sets
412 for (;number >0;number--) {
[a67d19]413 DoLog(1) && (Log() << Verbose(1) << "Determining data set " << number << " ... " << endl);
[042f82]414 // pick the point set
415 x = NULL;
[e138de]416 //PickRandomPointSet(T, LCList, x, N);
417 PickRandomNeighbouredPointSet(T, LCList, x, N);
[042f82]418
419 // calculate some sensible starting values for parameter fit
420 MaxDistance = 0.;
[273382]421 MinDistance = x[0].ScalarProduct(x[0]);
[042f82]422 for (int i=0;i<N;i++) {
[273382]423 distance = x[i].ScalarProduct(x[i]);
[042f82]424 if (distance > MaxDistance)
425 MaxDistance = distance;
426 if (distance < MinDistance)
427 MinDistance = distance;
428 }
[e138de]429 //Log() << Verbose(2) << "MinDistance " << MinDistance << ", MaxDistance " << MaxDistance << "." << endl;
[273382]430 EllipsoidCenter = Center; // use Center of Gravity as initial center of ellipsoid
[042f82]431 for (int i=0;i<3;i++)
432 EllipsoidAngle[i] = 0.;
433 EllipsoidLength[0] = sqrt(MaxDistance);
434 EllipsoidLength[1] = sqrt((MaxDistance+MinDistance)/2.);
435 EllipsoidLength[2] = sqrt(MinDistance);
436
437 // fit the parameters
[e138de]438 if (FitPointSetToEllipsoid(x, N, &EllipsoidCenter, &EllipsoidLength[0], &EllipsoidAngle[0])) {
[a67d19]439 DoLog(1) && (Log() << Verbose(1) << "Picking succeeded!" << endl);
[042f82]440 // output obtained parameter set
441 output << number << "\t";
442 for (int i=0;i<3;i++)
[0a4f7f]443 output << setprecision(9) << EllipsoidCenter[i] << "\t";
[042f82]444 for (int i=0;i<3;i++)
445 output << setprecision(9) << EllipsoidLength[i] << "\t";
446 for (int i=0;i<3;i++)
447 output << setprecision(9) << EllipsoidAngle[i] << "\t";
448 output << endl;
449 } else { // increase N to pick one more
[a67d19]450 DoLog(1) && (Log() << Verbose(1) << "Picking failed!" << endl);
[042f82]451 number++;
452 }
453 delete[](x); // free allocated memory for point set
454 }
455 // close output and finish
456 output.close();
457
[a67d19]458 DoLog(0) && (Log() << Verbose(0) << "End of FindDistributionOfEllipsoids" << endl);
[6ac7ee]459};
Note: See TracBrowser for help on using the repository browser.