Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Microstation DGN Access Library
4 : * Purpose: Code to stroke Arcs/Ellipses into polylines.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2001, Avenza Systems Inc, http://www.avenza.com/
9 : *
10 : * Permission is hereby granted, free of charge, to any person obtaining a
11 : * copy of this software and associated documentation files (the "Software"),
12 : * to deal in the Software without restriction, including without limitation
13 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 : * and/or sell copies of the Software, and to permit persons to whom the
15 : * Software is furnished to do so, subject to the following conditions:
16 : *
17 : * The above copyright notice and this permission notice shall be included
18 : * in all copies or substantial portions of the Software.
19 : *
20 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 : * DEALINGS IN THE SOFTWARE.
27 : ****************************************************************************/
28 :
29 : #include "dgnlibp.h"
30 : #include <cmath>
31 :
32 : constexpr double DEG_TO_RAD = M_PI / 180.0;
33 :
34 : /************************************************************************/
35 : /* ComputePointOnArc() */
36 : /************************************************************************/
37 :
38 438 : static void ComputePointOnArc2D(double dfPrimary, double dfSecondary,
39 : double dfAxisRotation, double dfAngle,
40 : double *pdfX, double *pdfY)
41 :
42 : {
43 : // dfAxisRotation and dfAngle are supposed to be in Radians
44 438 : const double dfCosRotation = cos(dfAxisRotation);
45 438 : const double dfSinRotation = sin(dfAxisRotation);
46 438 : const double dfEllipseX = dfPrimary * cos(dfAngle);
47 438 : const double dfEllipseY = dfSecondary * sin(dfAngle);
48 :
49 438 : *pdfX = dfEllipseX * dfCosRotation - dfEllipseY * dfSinRotation;
50 438 : *pdfY = dfEllipseX * dfSinRotation + dfEllipseY * dfCosRotation;
51 438 : }
52 :
53 : /************************************************************************/
54 : /* DGNStrokeArc() */
55 : /************************************************************************/
56 :
57 : /**
58 : * Generate a polyline approximation of an arc.
59 : *
60 : * Produce a series of equidistant (actually equi-angle) points along
61 : * an arc. Currently this only works for 2D arcs (and ellipses).
62 : *
63 : * @param hFile the DGN file to which the arc belongs (currently not used).
64 : * @param psArc the arc to be approximated.
65 : * @param nPoints the number of points to use to approximate the arc.
66 : * @param pasPoints the array of points into which to put the results.
67 : * There must be room for at least nPoints points.
68 : *
69 : * @return TRUE on success or FALSE on failure.
70 : */
71 :
72 6 : int DGNStrokeArc(CPL_UNUSED DGNHandle hFile, DGNElemArc *psArc, int nPoints,
73 : DGNPoint *pasPoints)
74 : {
75 6 : if (nPoints < 2)
76 0 : return FALSE;
77 :
78 6 : if (psArc->primary_axis == 0.0 || psArc->secondary_axis == 0.0)
79 : {
80 0 : CPLError(CE_Warning, CPLE_AppDefined,
81 : "Zero primary or secondary axis in DGNStrokeArc().");
82 0 : return FALSE;
83 : }
84 :
85 6 : const double dfAngleStep = psArc->sweepang / (nPoints - 1);
86 444 : for (int i = 0; i < nPoints; i++)
87 : {
88 438 : const double dfAngle = (psArc->startang + dfAngleStep * i) * DEG_TO_RAD;
89 :
90 438 : ComputePointOnArc2D(psArc->primary_axis, psArc->secondary_axis,
91 438 : psArc->rotation * DEG_TO_RAD, dfAngle,
92 438 : &(pasPoints[i].x), &(pasPoints[i].y));
93 438 : pasPoints[i].x += psArc->origin.x;
94 438 : pasPoints[i].y += psArc->origin.y;
95 438 : pasPoints[i].z = psArc->origin.z;
96 : }
97 :
98 6 : return TRUE;
99 : }
100 :
101 : /************************************************************************/
102 : /* DGNStrokeCurve() */
103 : /************************************************************************/
104 :
105 : /**
106 : * Generate a polyline approximation of an curve.
107 : *
108 : * Produce a series of equidistant points along a microstation curve element.
109 : * Currently this only works for 2D.
110 : *
111 : * @param hFile the DGN file to which the arc belongs (currently not used).
112 : * @param psCurve the curve to be approximated.
113 : * @param nPoints the number of points to use to approximate the curve.
114 : * @param pasPoints the array of points into which to put the results.
115 : * There must be room for at least nPoints points.
116 : *
117 : * @return TRUE on success or FALSE on failure.
118 : */
119 :
120 0 : int DGNStrokeCurve(CPL_UNUSED DGNHandle hFile, DGNElemMultiPoint *psCurve,
121 : int nPoints, DGNPoint *pasPoints)
122 : {
123 0 : const int nDGNPoints = psCurve->num_vertices;
124 :
125 0 : if (nDGNPoints < 6)
126 0 : return FALSE;
127 :
128 0 : if (nPoints < nDGNPoints - 4)
129 0 : return FALSE;
130 :
131 : /* -------------------------------------------------------------------- */
132 : /* Compute the Compute the slopes/distances of the segments. */
133 : /* -------------------------------------------------------------------- */
134 : double *padfMx =
135 0 : static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
136 : double *padfMy =
137 0 : static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
138 : double *padfD =
139 0 : static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
140 : double *padfTx =
141 0 : static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
142 : double *padfTy =
143 0 : static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
144 :
145 0 : double dfTotalD = 0.0;
146 :
147 0 : DGNPoint *pasDGNPoints = psCurve->vertices;
148 :
149 0 : for (int k = 0; k < nDGNPoints - 1; k++)
150 : {
151 : /* coverity[overrun-local] */
152 0 : padfD[k] = sqrt((pasDGNPoints[k + 1].x - pasDGNPoints[k].x) *
153 0 : (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) +
154 0 : (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) *
155 0 : (pasDGNPoints[k + 1].y - pasDGNPoints[k].y));
156 0 : if (padfD[k] == 0.0)
157 : {
158 0 : padfD[k] = 0.0001;
159 0 : padfMx[k] = 0.0;
160 0 : padfMy[k] = 0.0;
161 : }
162 : else
163 : {
164 0 : padfMx[k] = (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k];
165 0 : padfMy[k] = (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k];
166 : }
167 :
168 0 : if (k > 1 && k < nDGNPoints - 3)
169 0 : dfTotalD += padfD[k];
170 : }
171 :
172 : /* -------------------------------------------------------------------- */
173 : /* Compute the Tx, and Ty coefficients for each segment. */
174 : /* -------------------------------------------------------------------- */
175 0 : for (int k = 2; k < nDGNPoints - 2; k++)
176 : {
177 0 : if (fabs(padfMx[k + 1] - padfMx[k]) == 0.0 &&
178 0 : fabs(padfMx[k - 1] - padfMx[k - 2]) == 0.0)
179 : {
180 0 : padfTx[k] = (padfMx[k] + padfMx[k - 1]) / 2;
181 : }
182 : else
183 : {
184 0 : padfTx[k] = (padfMx[k - 1] * fabs(padfMx[k + 1] - padfMx[k]) +
185 0 : padfMx[k] * fabs(padfMx[k - 1] - padfMx[k - 2])) /
186 0 : (std::abs(padfMx[k + 1] - padfMx[k]) +
187 0 : std::abs(padfMx[k - 1] - padfMx[k - 2]));
188 : }
189 :
190 0 : if (fabs(padfMy[k + 1] - padfMy[k]) == 0.0 &&
191 0 : fabs(padfMy[k - 1] - padfMy[k - 2]) == 0.0)
192 : {
193 0 : padfTy[k] = (padfMy[k] + padfMy[k - 1]) / 2;
194 : }
195 : else
196 : {
197 0 : padfTy[k] = (padfMy[k - 1] * fabs(padfMy[k + 1] - padfMy[k]) +
198 0 : padfMy[k] * fabs(padfMy[k - 1] - padfMy[k - 2])) /
199 0 : (std::abs(padfMy[k + 1] - padfMy[k]) +
200 0 : std::abs(padfMy[k - 1] - padfMy[k - 2]));
201 : }
202 : }
203 :
204 : /* -------------------------------------------------------------------- */
205 : /* Determine a step size in D. We scale things so that we have */
206 : /* roughly equidistant steps in D, but assume we also want to */
207 : /* include every node along the way. */
208 : /* -------------------------------------------------------------------- */
209 0 : double dfStepSize = dfTotalD / (nPoints - (nDGNPoints - 4) - 1);
210 :
211 : /* ==================================================================== */
212 : /* Process each of the segments. */
213 : /* ==================================================================== */
214 0 : double dfD = dfStepSize;
215 0 : int iOutPoint = 0;
216 :
217 0 : for (int k = 2; k < nDGNPoints - 3; k++)
218 : {
219 : /* --------------------------------------------------------------------
220 : */
221 : /* Compute the "x" coefficients for this segment. */
222 : /* --------------------------------------------------------------------
223 : */
224 0 : const double dfCx = padfTx[k];
225 0 : const double dfBx =
226 0 : (3.0 * (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k] -
227 0 : 2.0 * padfTx[k] - padfTx[k + 1]) /
228 0 : padfD[k];
229 0 : const double dfAx =
230 0 : (padfTx[k] + padfTx[k + 1] -
231 0 : 2 * (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k]) /
232 0 : (padfD[k] * padfD[k]);
233 :
234 : /* --------------------------------------------------------------------
235 : */
236 : /* Compute the Y coefficients for this segment. */
237 : /* --------------------------------------------------------------------
238 : */
239 0 : const double dfCy = padfTy[k];
240 0 : const double dfBy =
241 0 : (3.0 * (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k] -
242 0 : 2.0 * padfTy[k] - padfTy[k + 1]) /
243 0 : padfD[k];
244 0 : const double dfAy =
245 0 : (padfTy[k] + padfTy[k + 1] -
246 0 : 2 * (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k]) /
247 0 : (padfD[k] * padfD[k]);
248 :
249 : /* --------------------------------------------------------------------
250 : */
251 : /* Add the start point for this segment. */
252 : /* --------------------------------------------------------------------
253 : */
254 0 : pasPoints[iOutPoint].x = pasDGNPoints[k].x;
255 0 : pasPoints[iOutPoint].y = pasDGNPoints[k].y;
256 0 : pasPoints[iOutPoint].z = 0.0;
257 0 : iOutPoint++;
258 :
259 : /* --------------------------------------------------------------------
260 : */
261 : /* Step along, adding intermediate points. */
262 : /* --------------------------------------------------------------------
263 : */
264 0 : while (dfD < padfD[k] && iOutPoint < nPoints - (nDGNPoints - k - 1))
265 : {
266 0 : pasPoints[iOutPoint].x = dfAx * dfD * dfD * dfD + dfBx * dfD * dfD +
267 0 : dfCx * dfD + pasDGNPoints[k].x;
268 0 : pasPoints[iOutPoint].y = dfAy * dfD * dfD * dfD + dfBy * dfD * dfD +
269 0 : dfCy * dfD + pasDGNPoints[k].y;
270 0 : pasPoints[iOutPoint].z = 0.0;
271 0 : iOutPoint++;
272 :
273 0 : dfD += dfStepSize;
274 : }
275 :
276 0 : dfD -= padfD[k];
277 : }
278 :
279 : /* -------------------------------------------------------------------- */
280 : /* Add the start point for this segment. */
281 : /* -------------------------------------------------------------------- */
282 0 : while (iOutPoint < nPoints)
283 : {
284 0 : pasPoints[iOutPoint].x = pasDGNPoints[nDGNPoints - 3].x;
285 0 : pasPoints[iOutPoint].y = pasDGNPoints[nDGNPoints - 3].y;
286 0 : pasPoints[iOutPoint].z = 0.0;
287 0 : iOutPoint++;
288 : }
289 :
290 : /* -------------------------------------------------------------------- */
291 : /* Cleanup. */
292 : /* -------------------------------------------------------------------- */
293 0 : CPLFree(padfMx);
294 0 : CPLFree(padfMy);
295 0 : CPLFree(padfD);
296 0 : CPLFree(padfTx);
297 0 : CPLFree(padfTy);
298 :
299 0 : return TRUE;
300 : }
301 :
302 : /************************************************************************/
303 : /* main() */
304 : /* */
305 : /* test mainline */
306 : /************************************************************************/
307 : #ifdef notdef
308 : int main(int argc, char **argv)
309 :
310 : {
311 : if (argc != 5)
312 : {
313 : printf( // ok
314 : "Usage: stroke primary_axis secondary_axis axis_rotation angle\n");
315 : exit(1);
316 : }
317 :
318 : const double dfPrimary = CPLAtof(argv[1]);
319 : const double dfSecondary = CPLAtof(argv[2]);
320 : const double dfAxisRotation = CPLAtof(argv[3]) / 180 * M_PI;
321 : const double dfAngle = CPLAtof(argv[4]) / 180 * M_PI;
322 :
323 : double dfX = 0.0;
324 : double dfY = 0.0;
325 : ComputePointOnArc2D(dfPrimary, dfSecondary, dfAxisRotation, dfAngle, &dfX,
326 : &dfY);
327 :
328 : printf("X=%.2f, Y=%.2f\n", dfX, dfY); // ok
329 :
330 : return 0;
331 : }
332 :
333 : #endif
|
