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 CPL_CVSID("$Id: dgnstroke.cpp 1dcdb5c8f617002a6e034d4d98dfdf080f709325 2018-01-14 01:50:09Z Kurt Schwehr $")
33
34 constexpr double DEG_TO_RAD = M_PI / 180.0;
35
36 /************************************************************************/
37 /* ComputePointOnArc() */
38 /************************************************************************/
39
ComputePointOnArc2D(double dfPrimary,double dfSecondary,double dfAxisRotation,double dfAngle,double * pdfX,double * pdfY)40 static void ComputePointOnArc2D( double dfPrimary, double dfSecondary,
41 double dfAxisRotation, double dfAngle,
42 double *pdfX, double *pdfY )
43
44 {
45 // dfAxisRotation and dfAngle are supposed to be in Radians
46 const double dfCosRotation = cos(dfAxisRotation);
47 const double dfSinRotation = sin(dfAxisRotation);
48 const double dfEllipseX = dfPrimary * cos(dfAngle);
49 const double dfEllipseY = dfSecondary * sin(dfAngle);
50
51 *pdfX = dfEllipseX * dfCosRotation - dfEllipseY * dfSinRotation;
52 *pdfY = dfEllipseX * dfSinRotation + dfEllipseY * dfCosRotation;
53 }
54
55 /************************************************************************/
56 /* DGNStrokeArc() */
57 /************************************************************************/
58
59 /**
60 * Generate a polyline approximation of an arc.
61 *
62 * Produce a series of equidistant (actually equi-angle) points along
63 * an arc. Currently this only works for 2D arcs (and ellipses).
64 *
65 * @param hFile the DGN file to which the arc belongs (currently not used).
66 * @param psArc the arc to be approximated.
67 * @param nPoints the number of points to use to approximate the arc.
68 * @param pasPoints the array of points into which to put the results.
69 * There must be room for at least nPoints points.
70 *
71 * @return TRUE on success or FALSE on failure.
72 */
73
DGNStrokeArc(CPL_UNUSED DGNHandle hFile,DGNElemArc * psArc,int nPoints,DGNPoint * pasPoints)74 int DGNStrokeArc( CPL_UNUSED DGNHandle hFile,
75 DGNElemArc *psArc,
76 int nPoints, DGNPoint * pasPoints )
77 {
78 if( nPoints < 2 )
79 return FALSE;
80
81 if( psArc->primary_axis == 0.0 || psArc->secondary_axis == 0.0 )
82 {
83 CPLError( CE_Warning, CPLE_AppDefined,
84 "Zero primary or secondary axis in DGNStrokeArc()." );
85 return FALSE;
86 }
87
88 const double dfAngleStep = psArc->sweepang / (nPoints - 1);
89 for( int i = 0; i < nPoints; i++ )
90 {
91 const double dfAngle = (psArc->startang + dfAngleStep * i) * DEG_TO_RAD;
92
93 ComputePointOnArc2D( psArc->primary_axis,
94 psArc->secondary_axis,
95 psArc->rotation * DEG_TO_RAD,
96 dfAngle,
97 &(pasPoints[i].x),
98 &(pasPoints[i].y) );
99 pasPoints[i].x += psArc->origin.x;
100 pasPoints[i].y += psArc->origin.y;
101 pasPoints[i].z = psArc->origin.z;
102 }
103
104 return TRUE;
105 }
106
107 /************************************************************************/
108 /* DGNStrokeCurve() */
109 /************************************************************************/
110
111 /**
112 * Generate a polyline approximation of an curve.
113 *
114 * Produce a series of equidistant points along a microstation curve element.
115 * Currently this only works for 2D.
116 *
117 * @param hFile the DGN file to which the arc belongs (currently not used).
118 * @param psCurve the curve to be approximated.
119 * @param nPoints the number of points to use to approximate the curve.
120 * @param pasPoints the array of points into which to put the results.
121 * There must be room for at least nPoints points.
122 *
123 * @return TRUE on success or FALSE on failure.
124 */
125
DGNStrokeCurve(CPL_UNUSED DGNHandle hFile,DGNElemMultiPoint * psCurve,int nPoints,DGNPoint * pasPoints)126 int DGNStrokeCurve( CPL_UNUSED DGNHandle hFile,
127 DGNElemMultiPoint *psCurve,
128 int nPoints, DGNPoint * pasPoints )
129 {
130 const int nDGNPoints = psCurve->num_vertices;
131
132 if( nDGNPoints < 6 )
133 return FALSE;
134
135 if( nPoints < nDGNPoints - 4 )
136 return FALSE;
137
138 /* -------------------------------------------------------------------- */
139 /* Compute the Compute the slopes/distances of the segments. */
140 /* -------------------------------------------------------------------- */
141 double *padfMx = static_cast<double *>(
142 CPLMalloc(sizeof(double) * nDGNPoints));
143 double *padfMy = static_cast<double *>(
144 CPLMalloc(sizeof(double) * nDGNPoints));
145 double *padfD = static_cast<double *>(
146 CPLMalloc(sizeof(double) * nDGNPoints));
147 double *padfTx = static_cast<double *>(
148 CPLMalloc(sizeof(double) * nDGNPoints));
149 double *padfTy = static_cast<double *>(
150 CPLMalloc(sizeof(double) * nDGNPoints));
151
152 double dfTotalD = 0.0;
153
154 DGNPoint *pasDGNPoints = psCurve->vertices;
155
156 for( int k = 0; k < nDGNPoints-1; k++ )
157 {
158 /* coverity[overrun-local] */
159 padfD[k] = sqrt( (pasDGNPoints[k+1].x-pasDGNPoints[k].x)
160 * (pasDGNPoints[k+1].x-pasDGNPoints[k].x)
161 + (pasDGNPoints[k+1].y-pasDGNPoints[k].y)
162 * (pasDGNPoints[k+1].y-pasDGNPoints[k].y) );
163 if( padfD[k] == 0.0 )
164 {
165 padfD[k] = 0.0001;
166 padfMx[k] = 0.0;
167 padfMy[k] = 0.0;
168 }
169 else
170 {
171 padfMx[k] = (pasDGNPoints[k+1].x - pasDGNPoints[k].x) / padfD[k];
172 padfMy[k] = (pasDGNPoints[k+1].y - pasDGNPoints[k].y) / padfD[k];
173 }
174
175 if( k > 1 && k < nDGNPoints - 3 )
176 dfTotalD += padfD[k];
177 }
178
179 /* -------------------------------------------------------------------- */
180 /* Compute the Tx, and Ty coefficients for each segment. */
181 /* -------------------------------------------------------------------- */
182 for( int k = 2; k < nDGNPoints - 2; k++ )
183 {
184 if( fabs(padfMx[k+1] - padfMx[k]) == 0.0
185 && fabs(padfMx[k-1] - padfMx[k-2]) == 0.0 )
186 {
187 padfTx[k] = (padfMx[k] + padfMx[k-1]) / 2;
188 }
189 else
190 {
191 padfTx[k] = (padfMx[k-1] * fabs( padfMx[k+1] - padfMx[k])
192 + padfMx[k] * fabs( padfMx[k-1] - padfMx[k-2] ))
193 / (std::abs(padfMx[k+1] - padfMx[k]) +
194 std::abs(padfMx[k-1] - padfMx[k-2]));
195 }
196
197 if( fabs(padfMy[k+1] - padfMy[k]) == 0.0
198 && fabs(padfMy[k-1] - padfMy[k-2]) == 0.0 )
199 {
200 padfTy[k] = (padfMy[k] + padfMy[k-1]) / 2;
201 }
202 else
203 {
204 padfTy[k] = (padfMy[k-1] * fabs( padfMy[k+1] - padfMy[k])
205 + padfMy[k] * fabs( padfMy[k-1] - padfMy[k-2] ))
206 / (std::abs(padfMy[k+1] - padfMy[k]) +
207 std::abs(padfMy[k-1] - padfMy[k-2]));
208 }
209 }
210
211 /* -------------------------------------------------------------------- */
212 /* Determine a step size in D. We scale things so that we have */
213 /* roughly equidistant steps in D, but assume we also want to */
214 /* include every node along the way. */
215 /* -------------------------------------------------------------------- */
216 double dfStepSize = dfTotalD / (nPoints - (nDGNPoints - 4) - 1);
217
218 /* ==================================================================== */
219 /* Process each of the segments. */
220 /* ==================================================================== */
221 double dfD = dfStepSize;
222 int iOutPoint = 0;
223
224 for( int k = 2; k < nDGNPoints - 3; k++ )
225 {
226 /* -------------------------------------------------------------------- */
227 /* Compute the "x" coefficients for this segment. */
228 /* -------------------------------------------------------------------- */
229 const double dfCx = padfTx[k];
230 const double dfBx = (3.0 * (pasDGNPoints[k+1].x - pasDGNPoints[k].x) / padfD[k]
231 - 2.0 * padfTx[k] - padfTx[k+1]) / padfD[k];
232 const double dfAx = (padfTx[k] + padfTx[k+1]
233 - 2 * (pasDGNPoints[k+1].x - pasDGNPoints[k].x) / padfD[k])
234 / (padfD[k] * padfD[k]);
235
236 /* -------------------------------------------------------------------- */
237 /* Compute the Y coefficients for this segment. */
238 /* -------------------------------------------------------------------- */
239 const double dfCy = padfTy[k];
240 const double dfBy = (3.0 * (pasDGNPoints[k+1].y - pasDGNPoints[k].y) / padfD[k]
241 - 2.0 * padfTy[k] - padfTy[k+1]) / padfD[k];
242 const double dfAy = (padfTy[k] + padfTy[k+1]
243 - 2 * (pasDGNPoints[k+1].y - pasDGNPoints[k].y) / padfD[k])
244 / (padfD[k] * padfD[k]);
245
246 /* -------------------------------------------------------------------- */
247 /* Add the start point for this segment. */
248 /* -------------------------------------------------------------------- */
249 pasPoints[iOutPoint].x = pasDGNPoints[k].x;
250 pasPoints[iOutPoint].y = pasDGNPoints[k].y;
251 pasPoints[iOutPoint].z = 0.0;
252 iOutPoint++;
253
254 /* -------------------------------------------------------------------- */
255 /* Step along, adding intermediate points. */
256 /* -------------------------------------------------------------------- */
257 while( dfD < padfD[k] && iOutPoint < nPoints - (nDGNPoints-k-1) )
258 {
259 pasPoints[iOutPoint].x = dfAx * dfD * dfD * dfD
260 + dfBx * dfD * dfD
261 + dfCx * dfD
262 + pasDGNPoints[k].x;
263 pasPoints[iOutPoint].y = dfAy * dfD * dfD * dfD
264 + dfBy * dfD * dfD
265 + dfCy * dfD
266 + pasDGNPoints[k].y;
267 pasPoints[iOutPoint].z = 0.0;
268 iOutPoint++;
269
270 dfD += dfStepSize;
271 }
272
273 dfD -= padfD[k];
274 }
275
276 /* -------------------------------------------------------------------- */
277 /* Add the start point for this segment. */
278 /* -------------------------------------------------------------------- */
279 while( iOutPoint < nPoints )
280 {
281 pasPoints[iOutPoint].x = pasDGNPoints[nDGNPoints-3].x;
282 pasPoints[iOutPoint].y = pasDGNPoints[nDGNPoints-3].y;
283 pasPoints[iOutPoint].z = 0.0;
284 iOutPoint++;
285 }
286
287 /* -------------------------------------------------------------------- */
288 /* Cleanup. */
289 /* -------------------------------------------------------------------- */
290 CPLFree( padfMx );
291 CPLFree( padfMy );
292 CPLFree( padfD );
293 CPLFree( padfTx );
294 CPLFree( padfTy );
295
296 return TRUE;
297 }
298
299 /************************************************************************/
300 /* main() */
301 /* */
302 /* test mainline */
303 /************************************************************************/
304 #ifdef notdef
main(int argc,char ** argv)305 int main( int argc, char ** argv )
306
307 {
308 if( argc != 5 )
309 {
310 printf( // ok
311 "Usage: stroke primary_axis secondary_axis axis_rotation angle\n");
312 exit( 1 );
313 }
314
315 const double dfPrimary = CPLAtof(argv[1]);
316 const double dfSecondary = CPLAtof(argv[2]);
317 const double dfAxisRotation = CPLAtof(argv[3]) / 180 * M_PI;
318 const double dfAngle = CPLAtof(argv[4]) / 180 * M_PI;
319
320 double dfX = 0.0;
321 double dfY = 0.0;
322 ComputePointOnArc2D( dfPrimary, dfSecondary, dfAxisRotation, dfAngle,
323 &dfX, &dfY );
324
325 printf( "X=%.2f, Y=%.2f\n", dfX, dfY ); // ok
326
327 return 0;
328 }
329
330 #endif
331