source: cpp/frams/model/geometry/geometryutils.cpp @ 229

Last change on this file since 229 was 193, checked in by Maciej Komosinski, 11 years ago

Set svn:eol-style native for all textual files
  • Property svn:eol-style set to native
File size: 8.5 KB
Line 
1// This file is a part of the Framsticks GDK.
2// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
3// Refer to http://www.framsticks.com/ for further information.
4
5#include "geometryutils.h"
6
7#include <math.h>
8
9double GeometryUtils::pointPosition(const int pointIndex, const int numberOfPoints)
10{
11        return pointIndex / (numberOfPoints-1.0);
12}
13
14double GeometryUtils::pointOnAxis(const double scale, const double position)
15{
16        return (position-0.5) * scale;
17}
18
19double GeometryUtils::pointOnAxis(const double scale, const int pointIndex, const int numberOfPoints)
20{
21        return pointOnAxis(scale, pointPosition(pointIndex, numberOfPoints));
22}
23
24double GeometryUtils::combination(const double value1, const double value2, const double position)
25{
26        return value1 + position * (value2-value1);
27}
28
29double GeometryUtils::combination(const double value1, const double value2, const int pointIndex, const int numberOfPoints)
30{
31        return combination(value1, value2, pointPosition(pointIndex, numberOfPoints));
32}
33
34bool GeometryUtils::isPointInsideModelExcludingPart(const Pt3D &point, const Model *model, const int excludedPartIndex)
35{
36        for (int i = 0; i < excludedPartIndex; i++)
37        {
38                if (isPointInsidePart(point, model->getPart(i)))
39                {
40                        return true;
41                }
42        }
43       
44        for (int i = excludedPartIndex+1; i < model->getPartCount(); i++)
45        {
46                if (isPointStrictlyInsidePart(point, model->getPart(i)))
47                {
48                        return true;
49                }
50        }
51       
52        return false;
53}
54
55bool GeometryUtils::isPointInsideModel(const Pt3D &point, const Model &model)
56{
57        for (int i = 0; i < model.getPartCount(); i++)
58        {
59                if (isPointInsidePart(point, model.getPart(i)))
60                {
61                        return true;
62                }
63        }
64       
65        return false;
66}
67
68bool GeometryUtils::isPointInsidePart(const Pt3D &point, const Part *part)
69{
70        switch (part->shape)
71        {
72                case Part::SHAPE_ELLIPSOID:
73                        return isPointInsideEllipsoid(point, part);
74                        break;
75                       
76                case Part::SHAPE_CUBOID:
77                        return isPointInsideCuboid(point, part);
78                        break;
79                       
80                case Part::SHAPE_CYLINDER:
81                        return isPointInsideCylinder(point, part);
82                        break;
83        }
84}
85
86bool GeometryUtils::isPointStrictlyInsidePart(const Pt3D &point, const Part *part)
87{
88        switch (part->shape)
89        {
90                case Part::SHAPE_ELLIPSOID:
91                        return isPointStrictlyInsideEllipsoid(point, part);
92                        break;
93                       
94                case Part::SHAPE_CUBOID:
95                        return isPointStrictlyInsideCuboid(point, part);
96                        break;
97                       
98                case Part::SHAPE_CYLINDER:
99                        return isPointStrictlyInsideCylinder(point, part);
100                        break;
101        }
102}
103
104bool GeometryUtils::isPointInsideEllipsoid(const Pt3D &point, const Part *part)
105{
106        Pt3D moved = point - part->p;
107        Pt3D rotated;
108        part->o.revTransform(rotated, moved);
109       
110        double r
111                = (pow(rotated.x, 2.0) / pow(part->scale.x, 2.0))
112                + (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
113                + (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
114       
115        return r <= 1.0;
116}
117
118bool GeometryUtils::isPointStrictlyInsideEllipsoid(const Pt3D &point, const Part *part)
119{
120        Pt3D moved = point - part->p;
121        Pt3D rotated;
122        part->o.revTransform(rotated, moved);
123       
124        double r
125                = (pow(rotated.x, 2.0) / pow(part->scale.x, 2.0))
126                + (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
127                + (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
128       
129        return r < 1.0;
130}
131
132bool GeometryUtils::isPointInsideCuboid(const Pt3D &point, const Part *part)
133{
134        Pt3D moved = point - part->p;
135        Pt3D rotated;
136        part->o.revTransform(rotated, moved);
137       
138        return (fabs(rotated.x) <= part->scale.x)
139                && (fabs(rotated.y) <= part->scale.y)
140                && (fabs(rotated.z) <= part->scale.z);
141}
142
143bool GeometryUtils::isPointStrictlyInsideCuboid(const Pt3D &point, const Part *part)
144{
145        Pt3D moved = point - part->p;
146        Pt3D rotated;
147        part->o.revTransform(rotated, moved);
148       
149        return (fabs(rotated.x) < part->scale.x)
150                && (fabs(rotated.y) < part->scale.y)
151                && (fabs(rotated.z) < part->scale.z);
152}
153
154bool GeometryUtils::isPointInsideCylinder(const Pt3D &point, const Part *part)
155{
156        Pt3D moved = point - part->p;
157        Pt3D rotated;
158        part->o.revTransform(rotated, moved);
159       
160        double r
161                = (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
162                + (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
163       
164        return (fabs(rotated.x) <= part->scale.x) && (r <= 1.0);
165}
166
167bool GeometryUtils::isPointStrictlyInsideCylinder(const Pt3D &point, const Part *part)
168{
169        Pt3D moved = point - part->p;
170        Pt3D rotated;
171        part->o.revTransform(rotated, moved);
172       
173        double r
174                = (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
175                + (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
176       
177        return (fabs(rotated.x) < part->scale.x) && (r < 1.0);
178}
179
180void GeometryUtils::findSizesAndAxesOfPointsGroup(SListTempl<Pt3D> &points, Pt3D &sizes,
181        Orient &axes)
182{
183        findSizeAndAxisOfPointsGroup(points, sizes.x, axes.x);
184        orthographicProjectionToPlane(points, axes.x);
185        findSizeAndAxisOfPointsGroup(points, sizes.y, axes.y);
186        orthographicProjectionToPlane(points, axes.y);
187       
188        Pt3D min(points.get(0)), max(points.get(0));
189       
190        for (int i = 1; i < points.size(); i++)
191        {
192                min.getMin(points.get(i));
193                max.getMax(points.get(i));
194        }
195       
196        sizes.z = min.distanceTo(max);
197        axes.z.vectorProduct(axes.x, axes.y);
198}
199
200void GeometryUtils::findSizeAndAxisOfPointsGroup(const SListTempl<Pt3D> &points, double &size,
201        Pt3D &axis)
202{
203        int index1, index2;
204        size = findTwoFurthestPoints(points, index1, index2);
205        createAxisFromTwoPoints(axis, points.get(index1), points.get(index2));
206}
207
208double GeometryUtils::findTwoFurthestPoints(const SListTempl<Pt3D> &points, int &index1,
209        int &index2)
210{
211        double distance = 0;
212        index1 = index2 = 0;
213       
214        for (int i = 0; i < points.size()-1; i++)
215        {
216                Pt3D p1 = points.get(i);
217               
218                for (int j = i+1; j < points.size(); j++)
219                {
220                        Pt3D p2 = points.get(j);
221                        double d = p1.distanceTo(p2);
222                       
223                        if (d > distance)
224                        {
225                                distance = d;
226                                index1 = i;
227                                index2 = j;
228                        }
229                }
230        }
231       
232        return distance;
233}
234
235void GeometryUtils::createAxisFromTwoPoints(Pt3D &axis, const Pt3D &point1, const Pt3D &point2)
236{
237        Pt3D vector = point2 - point1;
238        vector.normalize();
239       
240        axis.x = vector.x;
241        axis.y = vector.y;
242        axis.z = vector.z;
243}
244
245void GeometryUtils::orthographicProjectionToPlane(SListTempl<Pt3D> &points,
246        const Pt3D &planeNormalVector)
247{
248        for (int i = 0; i < points.size(); i++)
249        {
250                Pt3D &point = points.get(i);
251               
252                double distance = pointDistanceToPlane(point, planeNormalVector);
253               
254                point.x -= planeNormalVector.x * distance;
255                point.y -= planeNormalVector.y * distance;
256                point.z -= planeNormalVector.z * distance;
257        }
258}
259
260double GeometryUtils::pointDistanceToPlane(const Pt3D &point, const Pt3D &planeNormalVector)
261{
262        return planeNormalVector.x*point.x + planeNormalVector.y*point.y + planeNormalVector.z*point.z;
263}
264
265void GeometryUtils::getRectangleApicesFromCuboid(const Part *part, const CuboidFaces::Face face, Pt3D &apex1, Pt3D &apex2, Pt3D &apex3, Pt3D &apex4)
266{
267        Pt3D temp1(part->scale), temp2(part->scale), temp3(part->scale), temp4(part->scale);
268       
269        if (CuboidFaces::isX(face))
270        {
271                temp2.z *= -1;
272                temp3.y *= -1;
273                temp4.z *= -1;
274                temp4.y *= -1;
275        }
276        else if (CuboidFaces::isY(face))
277        {
278                temp2.x *= -1;
279                temp3.z *= -1;
280                temp4.x *= -1;
281                temp4.z *= -1;
282        }
283        else if (CuboidFaces::isZ(face))
284        {
285                temp2.y *= -1;
286                temp3.x *= -1;
287                temp4.y *= -1;
288                temp4.x *= -1;
289        }
290       
291        if (CuboidFaces::isNegative(face))
292        {
293                temp1 *= -1;
294                temp2 *= -1;
295                temp3 *= -1;
296                temp4 *= -1;
297        }
298
299        part->o.transform(apex1, temp1);
300        part->o.transform(apex2, temp2);
301        part->o.transform(apex3, temp3);
302        part->o.transform(apex4, temp4);
303
304        apex1 += part->p;
305        apex2 += part->p;
306        apex3 += part->p;
307        apex4 += part->p;
308}
309
310void GeometryUtils::getRectangleApices(const double width, const double height, const Pt3D &position, const Orient &orient, Pt3D &apex1, Pt3D &apex2, Pt3D &apex3, Pt3D &apex4)
311{
312        Pt3D temp1(0.0, +width, +height);
313        Pt3D temp2(0.0, +width, -height);
314        Pt3D temp3(0.0, -width, +height);
315        Pt3D temp4(0.0, -width, -height);
316
317        orient.transform(apex1, temp1);
318        orient.transform(apex2, temp2);
319        orient.transform(apex3, temp3);
320        orient.transform(apex4, temp4);
321
322        apex1 += position;
323        apex2 += position;
324        apex3 += position;
325        apex4 += position;
326}
327
328void GeometryUtils::getNextEllipseSegmentationPoint(const double d, const double a, const double b, double &x, double &y)
329{
330        x += d / sqrt(1.0 + (b*b * x*x) / (a*a * (a*a - x*x)));
331        y = b * sqrt(1.0 - (x*x) / (a*a));
332}
333
334double GeometryUtils::ellipsoidArea(const Pt3D &sizes)
335{
336        return ellipsoidArea(sizes.x, sizes.y, sizes.z);
337}
338
339double GeometryUtils::ellipsoidArea(const double a, const double b, const double c)
340{
341        double p = 1.6075;
342        double ap = pow(a, p);
343        double bp = pow(b, p);
344        double cp = pow(c, p);
345        return 4*M_PI * pow((ap*bp + bp*cp + cp*ap) / 3.0, 1.0 / p);
346}
347
348double GeometryUtils::ellipsePerimeter(const double a, const double b)
349{
350        return M_PI * ((3 * (a+b)) - sqrt((3*a + b) * (a + 3*b)));
351}
Note: See TracBrowser for help on using the repository browser.