1 /*
2  Copyright (C) 2010-2014 Kristian Duske
3 
4  This file is part of TrenchBroom.
5 
6  TrenchBroom is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  TrenchBroom is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with TrenchBroom. If not, see <http://www.gnu.org/licenses/>.
18  */
19 
20 #include "PlanePointFinder.h"
21 
22 namespace TrenchBroom {
23     namespace Model {
24         class GridSearchCursor {
25         private:
26             static const size_t Center = 4;
27 
28             static const Vec2 MoveOffsets[];
29 
30             const Plane3& m_plane;
31             const FloatType m_frequency;
32 
33             Vec2 m_position;
34             FloatType m_errors[9];
35         public:
GridSearchCursor(const Plane3 & plane,const FloatType frequency)36             GridSearchCursor(const Plane3& plane, const FloatType frequency) :
37             m_plane(plane),
38             m_frequency(frequency) {
39                 for (size_t i = 0; i < 9; ++i)
40                     m_errors[i] = 0.0;
41             }
42 
findMinimum(const Vec3 & initialPosition)43             Vec3 findMinimum(const Vec3& initialPosition) {
44                 for (size_t i = 0; i < 2; ++i)
45                     m_position[i] = Math::round(initialPosition[i]);
46 
47                 findLocalMinimum();
48                 const Vec2 localMinPos = m_position;
49                 const FloatType localMinErr = m_errors[Center];
50 
51                 Vec2 globalMinPos = localMinPos;
52                 FloatType globalMinErr = localMinErr;
53 
54                 if (globalMinErr > 0.0) {
55                     // To escape local minima, let's search some adjacent quadrants
56                     // The number of extra quadrants should depend on the frequency: The higher the frequency, the
57                     // more quadrants should be searched.
58                     const size_t numQuadrants = static_cast<size_t>(std::ceil(m_frequency * m_frequency * 3.0));
59                     for (size_t i = 0; i < numQuadrants && globalMinErr > 0.0; ++i) {
60                         if (i != Center) {
61                             m_position = localMinPos + i * 3.0 * MoveOffsets[i];
62                             findLocalMinimum();
63                             const FloatType newError = m_errors[Center];
64                             if (newError < globalMinErr) {
65                                 globalMinPos = m_position;
66                                 globalMinErr = newError;
67                             }
68                         }
69                     }
70                 }
71 
72                 return Vec3(globalMinPos.x(),
73                             globalMinPos.y(),
74                             Math::round(m_plane.zAt(globalMinPos)));
75             }
76         private:
findLocalMinimum()77             void findLocalMinimum() {
78                 updateErrors();
79 
80                 size_t smallestError = findSmallestError();
81                 while (smallestError != Center)
82                     smallestError = moveCursor(smallestError);
83             }
84 
moveCursor(const size_t direction)85             size_t moveCursor(const size_t direction) {
86                 m_position += MoveOffsets[direction];
87                 updateErrors();
88                 return findSmallestError();
89             }
90 
updateErrors()91             void updateErrors() {
92                 for (size_t i = 0; i < 9; ++i)
93                     m_errors[i] = computeError(i);
94             }
95 
computeError(const size_t location) const96             FloatType computeError(const size_t location) const {
97                 const FloatType z = m_plane.zAt(m_position + MoveOffsets[location]);
98                 return std::abs(z - Math::round(z));
99             }
100 
findSmallestError()101             size_t findSmallestError() {
102                 size_t smallest = Center;
103                 for (size_t i = 0; i < 9; ++i) {
104                     if (m_errors[i] < m_errors[smallest])
105                         smallest = i;
106                 }
107                 return smallest;
108             }
109         };
110 
111         const Vec2 GridSearchCursor::MoveOffsets[] = {
112             Vec2(-1.0,  1.0), Vec2( 0.0,  1.0), Vec2( 1.0,  1.0),
113             Vec2(-1.0,  0.0), Vec2( 0.0,  0.0), Vec2( 1.0,  0.0),
114             Vec2(-1.0, -1.0), Vec2( 0.0, -1.0), Vec2( 1.0, -1.0)
115         };
116 
117         FloatType computePlaneFrequency(const Plane3& plane);
118         void setDefaultPlanePoints(const Plane3& plane, BrushFace::Points& points);
119 
computePlaneFrequency(const Plane3 & plane)120         FloatType computePlaneFrequency(const Plane3& plane) {
121             static const FloatType c = 1.0 - std::sin(Math::C::pi() / 4.0);
122 
123             const Vec3& axis = plane.normal.firstAxis();
124             const FloatType d = plane.normal.dot(axis);
125             assert(d != 0.0);
126 
127             return (1.0 - d) / c;
128         }
129 
setDefaultPlanePoints(const Plane3 & plane,BrushFace::Points & points)130         void setDefaultPlanePoints(const Plane3& plane, BrushFace::Points& points) {
131             points[0] = plane.anchor().rounded();
132             const Math::Axis::Type axis = plane.normal.firstComponent();
133             switch (axis) {
134                 case Math::Axis::AX:
135                     if (plane.normal.x() > 0.0) {
136                         points[1] = points[0] + 64.0 * Vec3::PosZ;
137                         points[2] = points[0] + 64.0 * Vec3::PosY;
138                     } else {
139                         points[1] = points[0] + 64.0 * Vec3::PosY;
140                         points[2] = points[0] + 64.0 * Vec3::PosZ;
141                     }
142                     break;
143                 case Math::Axis::AY:
144                     if (plane.normal.y() > 0.0) {
145                         points[1] = points[0] + 64.0 * Vec3::PosX;
146                         points[2] = points[0] + 64.0 * Vec3::PosZ;
147                     } else {
148                         points[1] = points[0] + 64.0 * Vec3::PosZ;
149                         points[2] = points[0] + 64.0 * Vec3::PosX;
150                     }
151                     break;
152                 default:
153                     if  (plane.normal.z() > 0.0) {
154                         points[1] = points[0] + 64.0 * Vec3::PosY;
155                         points[2] = points[0] + 64.0 * Vec3::PosX;
156                     } else {
157                         points[1] = points[0] + 64.0 * Vec3::PosX;
158                         points[2] = points[0] + 64.0 * Vec3::PosY;
159                     }
160                     break;
161             }
162         }
163 
findPoints(const Plane3 & plane,BrushFace::Points & points,const size_t numPoints)164         void PlanePointFinder::findPoints(const Plane3& plane, BrushFace::Points& points, const size_t numPoints) {
165             using std::swap;
166 
167             assert(numPoints <= 3);
168 
169             if (numPoints == 3 && points[0].isInteger() && points[1].isInteger() && points[2].isInteger())
170                 return;
171 
172             const FloatType frequency = computePlaneFrequency(plane);
173             if (Math::zero(frequency, 1.0 / 7084.0)) {
174                 setDefaultPlanePoints(plane, points);
175                 return;
176             }
177 
178             const Math::Axis::Type axis = plane.normal.firstComponent();
179             const Plane3 swizzledPlane(plane.distance, swizzle(plane.normal, axis));
180             for (size_t i = 0; i < 3; ++i)
181                 points[i] = swizzle(points[i], axis);
182 
183             const FloatType waveLength = 1.0 / frequency;
184             const FloatType pointDistance = std::min(64.0, waveLength);
185 
186             FloatType multiplier = 10.0;
187             GridSearchCursor cursor(swizzledPlane, frequency);
188             if (numPoints == 0)
189                 points[0] = cursor.findMinimum(swizzledPlane.anchor());
190             else if (!points[0].isInteger())
191                 points[0] = cursor.findMinimum(points[0]);
192 
193             Vec3 v1, v2;
194             FloatType cos;
195             size_t count = 0;
196             do {
197                 if (numPoints < 2 || !points[1].isInteger())
198                     points[1] = cursor.findMinimum(points[0] + 0.33 * multiplier * pointDistance * Vec3::PosX);
199                 points[2] = cursor.findMinimum(points[0] + multiplier * (pointDistance * Vec3::PosY - pointDistance / 2.0 * Vec3::PosX));
200                 v1 = points[2] - points[0];
201                 v2 = points[1] - points[0];
202                 v1.normalize();
203                 v2.normalize();
204                 cos = v1.dot(v2);
205                 multiplier *= 1.5f;
206                 ++count;
207             } while (Math::isnan(cos) || std::abs(cos) > 0.9);
208 
209             cross(v1, v2);
210             if ((v1.z() > 0.0) != (swizzledPlane.normal.z() > 0.0))
211                 swap(points[0], points[2]);
212 
213             for (size_t i = 0; i < 3; ++i)
214                 points[i] = unswizzle(points[i], axis);
215         }
216     }
217 }
218