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