1 /******************************************************************************
2 * $Id: ntfstroke.cpp 10645 2007-01-18 02:22:39Z warmerdam $
3 *
4 * Project: NTF Translator
5 * Purpose: NTF Arc to polyline stroking code. This code is really generic,
6 * and might be moved into an OGR module at some point in the
7 * future.
8 * Author: Frank Warmerdam, warmerdam@pobox.com
9 *
10 ******************************************************************************
11 * Copyright (c) 2001, Frank Warmerdam
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining a
14 * copy of this software and associated documentation files (the "Software"),
15 * to deal in the Software without restriction, including without limitation
16 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 * and/or sell copies of the Software, and to permit persons to whom the
18 * Software is furnished to do so, subject to the following conditions:
19 *
20 * The above copyright notice and this permission notice shall be included
21 * in all copies or substantial portions of the Software.
22 *
23 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 * DEALINGS IN THE SOFTWARE.
30 ****************************************************************************/
31
32 #include <stdarg.h>
33 #include "ntf.h"
34 #include "cpl_conv.h"
35 #include "cpl_string.h"
36
37 CPL_CVSID("$Id: ntfstroke.cpp 10645 2007-01-18 02:22:39Z warmerdam $");
38
39 #ifndef PI
40 #define PI 3.14159265358979323846
41 #endif
42
43 /************************************************************************/
44 /* NTFArcCenterFromEdgePoints() */
45 /* */
46 /* Compute the center of an arc/circle from three edge points. */
47 /************************************************************************/
48
NTFArcCenterFromEdgePoints(double x_c0,double y_c0,double x_c1,double y_c1,double x_c2,double y_c2,double * x_center,double * y_center)49 int NTFArcCenterFromEdgePoints( double x_c0, double y_c0,
50 double x_c1, double y_c1,
51 double x_c2, double y_c2,
52 double *x_center, double *y_center )
53
54 {
55
56 /* -------------------------------------------------------------------- */
57 /* Handle a degenerate case that occurs in OSNI products by */
58 /* making some assumptions. If the first and third points are */
59 /* the same assume they are intended to define a full circle, */
60 /* and that the second point is on the opposite side of the */
61 /* circle. */
62 /* -------------------------------------------------------------------- */
63 if( x_c0 == x_c2 && y_c0 == y_c2 )
64 {
65 *x_center = (x_c0 + x_c1) * 0.5;
66 *y_center = (y_c0 + y_c1) * 0.5;
67
68 return TRUE;
69 }
70
71 /* -------------------------------------------------------------------- */
72 /* Compute the inverse of the slopes connecting the first and */
73 /* second points. Also compute the center point of the two */
74 /* lines ... the point our crossing line will go through. */
75 /* -------------------------------------------------------------------- */
76 double m1, x1, y1;
77
78 if( (y_c1 - y_c0) != 0.0 )
79 m1 = (x_c0 - x_c1) / (y_c1 - y_c0);
80 else
81 m1 = 1e+10;
82
83 x1 = (x_c0 + x_c1) * 0.5;
84 y1 = (y_c0 + y_c1) * 0.5;
85
86 /* -------------------------------------------------------------------- */
87 /* Compute the same for the second point compared to the third */
88 /* point. */
89 /* -------------------------------------------------------------------- */
90 double m2, x2, y2;
91
92 if( (y_c2 - y_c1) != 0.0 )
93 m2 = (x_c1 - x_c2) / (y_c2 - y_c1);
94 else
95 m2 = 1e+10;
96
97 x2 = (x_c1 + x_c2) * 0.5;
98 y2 = (y_c1 + y_c2) * 0.5;
99
100 /* -------------------------------------------------------------------- */
101 /* Turn these into the Ax+By+C = 0 form of the lines. */
102 /* -------------------------------------------------------------------- */
103 double a1, a2, b1, b2, c1, c2;
104
105 a1 = m1;
106 a2 = m2;
107
108 b1 = -1.0;
109 b2 = -1.0;
110
111 c1 = (y1 - m1*x1);
112 c2 = (y2 - m2*x2);
113
114 /* -------------------------------------------------------------------- */
115 /* Compute the intersection of the two lines through the center */
116 /* of the circle, using Kramers rule. */
117 /* -------------------------------------------------------------------- */
118 double det_inv;
119
120 if( a1*b2 - a2*b1 == 0.0 )
121 return FALSE;
122
123 det_inv = 1 / (a1*b2 - a2*b1);
124
125 *x_center = (b1*c2 - b2*c1) * det_inv;
126 *y_center = (a2*c1 - a1*c2) * det_inv;
127
128 return TRUE;
129 }
130
131 /************************************************************************/
132 /* NTFStrokeArcToOGRGeometry_Points() */
133 /************************************************************************/
134
135 OGRGeometry *
NTFStrokeArcToOGRGeometry_Points(double dfStartX,double dfStartY,double dfAlongX,double dfAlongY,double dfEndX,double dfEndY,int nVertexCount)136 NTFStrokeArcToOGRGeometry_Points( double dfStartX, double dfStartY,
137 double dfAlongX, double dfAlongY,
138 double dfEndX, double dfEndY,
139 int nVertexCount )
140
141 {
142 double dfStartAngle, dfEndAngle, dfAlongAngle;
143 double dfCenterX, dfCenterY, dfRadius;
144
145 if( !NTFArcCenterFromEdgePoints( dfStartX, dfStartY, dfAlongX, dfAlongY,
146 dfEndX, dfEndY, &dfCenterX, &dfCenterY ) )
147 return NULL;
148
149 if( dfStartX == dfEndX && dfStartY == dfEndY )
150 {
151 dfStartAngle = 0.0;
152 dfEndAngle = 360.0;
153 }
154 else
155 {
156 double dfDeltaX, dfDeltaY;
157
158 dfDeltaX = dfStartX - dfCenterX;
159 dfDeltaY = dfStartY - dfCenterY;
160 dfStartAngle = atan2(dfDeltaY,dfDeltaX) * 180.0 / PI;
161
162 dfDeltaX = dfAlongX - dfCenterX;
163 dfDeltaY = dfAlongY - dfCenterY;
164 dfAlongAngle = atan2(dfDeltaY,dfDeltaX) * 180.0 / PI;
165
166 dfDeltaX = dfEndX - dfCenterX;
167 dfDeltaY = dfEndY - dfCenterY;
168 dfEndAngle = atan2(dfDeltaY,dfDeltaX) * 180.0 / PI;
169
170 #ifdef notdef
171 if( dfStartAngle > dfAlongAngle && dfAlongAngle > dfEndAngle )
172 {
173 double dfTempAngle;
174
175 dfTempAngle = dfStartAngle;
176 dfStartAngle = dfEndAngle;
177 dfEndAngle = dfTempAngle;
178 }
179 #endif
180
181 while( dfAlongAngle < dfStartAngle )
182 dfAlongAngle += 360.0;
183
184 while( dfEndAngle < dfAlongAngle )
185 dfEndAngle += 360.0;
186
187 if( dfEndAngle - dfStartAngle > 360.0 )
188 {
189 double dfTempAngle;
190
191 dfTempAngle = dfStartAngle;
192 dfStartAngle = dfEndAngle;
193 dfEndAngle = dfTempAngle;
194
195 while( dfEndAngle < dfStartAngle )
196 dfStartAngle -= 360.0;
197 }
198 }
199
200 dfRadius = sqrt( (dfCenterX - dfStartX) * (dfCenterX - dfStartX)
201 + (dfCenterY - dfStartY) * (dfCenterY - dfStartY) );
202
203 return NTFStrokeArcToOGRGeometry_Angles( dfCenterX, dfCenterY,
204 dfRadius,
205 dfStartAngle, dfEndAngle,
206 nVertexCount );
207 }
208
209 /************************************************************************/
210 /* NTFStrokeArcToOGRGeometry_Angles() */
211 /************************************************************************/
212
213 OGRGeometry *
NTFStrokeArcToOGRGeometry_Angles(double dfCenterX,double dfCenterY,double dfRadius,double dfStartAngle,double dfEndAngle,int nVertexCount)214 NTFStrokeArcToOGRGeometry_Angles( double dfCenterX, double dfCenterY,
215 double dfRadius,
216 double dfStartAngle, double dfEndAngle,
217 int nVertexCount )
218
219 {
220 OGRLineString *poLine = new OGRLineString;
221 double dfArcX, dfArcY, dfSlice;
222 int iPoint;
223
224 nVertexCount = MAX(2,nVertexCount);
225 dfSlice = (dfEndAngle-dfStartAngle)/(nVertexCount-1);
226
227 poLine->setNumPoints( nVertexCount );
228
229 for( iPoint=0; iPoint < nVertexCount; iPoint++ )
230 {
231 double dfAngle;
232
233 dfAngle = (dfStartAngle + iPoint * dfSlice) * PI / 180.0;
234
235 dfArcX = dfCenterX + cos(dfAngle) * dfRadius;
236 dfArcY = dfCenterY + sin(dfAngle) * dfRadius;
237
238 poLine->setPoint( iPoint, dfArcX, dfArcY );
239 }
240
241 return poLine;
242 }
243
244
245