1 : /******************************************************************************
2 : * $Id: ilihelper.cpp 17904 2009-10-26 19:22:46Z chaitanya $
3 : *
4 : * Project: Interlis 1/2 Translator
5 : * Purpose: Helper functions for Interlis reader
6 : * Author: Pirmin Kalberer, Sourcepole AG
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Pirmin Kalberer, Sourcepole AG
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 :
31 : #include "ilihelper.h"
32 :
33 :
34 : CPL_CVSID("$Id: ilihelper.cpp 17904 2009-10-26 19:22:46Z chaitanya $");
35 :
36 :
37 0 : OGRPoint *getARCCenter(OGRPoint *ptStart, OGRPoint *ptArc, OGRPoint *ptEnd) {
38 : // FIXME precision
39 0 : double bx = ptStart->getX(); double by = ptStart->getY();
40 0 : double cx = ptArc->getX(); double cy = ptArc->getY();
41 0 : double dx = ptEnd->getX(); double dy = ptEnd->getY();
42 : double temp, bc, cd, det, x, y;
43 :
44 0 : temp = cx*cx+cy*cy;
45 0 : bc = (bx*bx + by*by - temp)/2.0;
46 0 : cd = (temp - dx*dx - dy*dy)/2.0;
47 0 : det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
48 :
49 0 : OGRPoint *center = new OGRPoint();
50 :
51 0 : if (fabs(det) < 1.0e-6) { // could not determin the determinante: too small
52 0 : return center;
53 : }
54 0 : det = 1/det;
55 0 : x = (bc*(cy-dy)-cd*(by-cy))*det;
56 0 : y = ((bx-cx)*cd-(cx-dx)*bc)*det;
57 :
58 0 : center->setX(x);
59 0 : center->setY(y);
60 :
61 0 : return center;
62 : }
63 :
64 0 : void interpolateArc(OGRLineString* line, OGRPoint *ptStart, OGRPoint *ptOnArc, OGRPoint *ptEnd, double arcIncr) {
65 0 : OGRPoint *center = getARCCenter(ptStart, ptOnArc, ptEnd);
66 :
67 0 : double cx = center->getX(); double cy = center->getY();
68 0 : double px = ptOnArc->getX(); double py = ptOnArc->getY();
69 0 : double r = sqrt((cx-px)*(cx-px)+(cy-py)*(cy-py));
70 :
71 0 : double myAlpha = 2.0*acos(1.0-0.002/r);
72 :
73 0 : if (myAlpha < arcIncr) {
74 0 : arcIncr = myAlpha;
75 : }
76 :
77 0 : double a1 = atan2(ptStart->getY() - cy, ptStart->getX() - cx);
78 0 : double a2 = atan2(py - cy, px - cx);
79 0 : double a3 = atan2(ptEnd->getY() - cy, ptEnd->getX() - cx);
80 :
81 : double sweep;
82 :
83 : // Clockwise
84 0 : if(a1 > a2 && a2 > a3) {
85 0 : sweep = a3 - a1;
86 : }
87 : // Counter-clockwise
88 0 : else if(a1 < a2 && a2 < a3) {
89 0 : sweep = a3 - a1;
90 : }
91 : // Clockwise, wrap
92 0 : else if((a1 < a2 && a1 > a3) || (a2 < a3 && a1 > a3)) {
93 0 : sweep = a3 - a1 + 2*PI;
94 : }
95 : // Counter-clockwise, wrap
96 0 : else if((a1 > a2 && a1 < a3) || (a2 > a3 && a1 < a3)) {
97 0 : sweep = a3 - a1 - 2*PI;
98 : }
99 : else {
100 0 : sweep = 0.0;
101 : }
102 :
103 0 : int ptcount = ceil(fabs(sweep/arcIncr));
104 :
105 0 : if(sweep < 0) arcIncr *= -1.0;
106 :
107 0 : double angle = a1;
108 :
109 0 : for(int i = 0; i < ptcount - 1; i++) {
110 0 : angle += arcIncr;
111 :
112 0 : if(arcIncr > 0.0 && angle > PI) angle -= 2*PI;
113 0 : if(arcIncr < 0.0 && angle < -1*PI) angle -= 2*PI;
114 :
115 0 : double x = cx + r*cos(angle);
116 0 : double y = cy + r*sin(angle);
117 :
118 0 : line->addPoint(x, y, 0);
119 :
120 :
121 0 : if((angle < a2) && ((angle + arcIncr) > a2)) {
122 0 : line->addPoint(ptOnArc);
123 : }
124 :
125 0 : if((angle > a2) && ((angle + arcIncr) < a2)) {
126 0 : line->addPoint(ptOnArc);
127 : }
128 :
129 : }
130 0 : line->addPoint(ptEnd);
131 0 : delete center;
132 0 : }
133 :
|
