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