1 : /******************************************************************************
2 : * $Id: ilihelper.cpp 24609 2012-06-25 15:21:10Z pka $
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 24609 2012-06-25 15:21:10Z pka $");
35 :
36 :
37 15 : OGRPoint *getARCCenter(OGRPoint *ptStart, OGRPoint *ptArc, OGRPoint *ptEnd) {
38 : // FIXME precision
39 15 : double bx = ptStart->getX(); double by = ptStart->getY();
40 15 : double cx = ptArc->getX(); double cy = ptArc->getY();
41 15 : double dx = ptEnd->getX(); double dy = ptEnd->getY();
42 : double temp, bc, cd, det, x, y;
43 :
44 15 : temp = cx*cx+cy*cy;
45 15 : bc = (bx*bx + by*by - temp)/2.0;
46 15 : cd = (temp - dx*dx - dy*dy)/2.0;
47 15 : det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
48 :
49 15 : OGRPoint *center = new OGRPoint();
50 :
51 15 : if (fabs(det) < 1.0e-6) { // could not determin the determinante: too small
52 0 : return center;
53 : }
54 15 : det = 1/det;
55 15 : x = (bc*(cy-dy)-cd*(by-cy))*det;
56 15 : y = ((bx-cx)*cd-(cx-dx)*bc)*det;
57 :
58 15 : center->setX(x);
59 15 : center->setY(y);
60 :
61 15 : return center;
62 : }
63 :
64 15 : void interpolateArc(OGRLineString* line, OGRPoint *ptStart, OGRPoint *ptOnArc, OGRPoint *ptEnd, double arcIncr) {
65 15 : OGRPoint *center = getARCCenter(ptStart, ptOnArc, ptEnd);
66 :
67 15 : double cx = center->getX(); double cy = center->getY();
68 15 : double px = ptOnArc->getX(); double py = ptOnArc->getY();
69 15 : double r = sqrt((cx-px)*(cx-px)+(cy-py)*(cy-py));
70 :
71 : //assure minimal chord length (0.002m???)
72 15 : double myAlpha = 2.0*acos(1.0-0.002/r);
73 15 : if (myAlpha < arcIncr) {
74 0 : arcIncr = myAlpha;
75 : }
76 :
77 15 : double a1 = atan2(ptStart->getY() - cy, ptStart->getX() - cx);
78 15 : double a2 = atan2(py - cy, px - cx);
79 15 : double a3 = atan2(ptEnd->getY() - cy, ptEnd->getX() - cx);
80 :
81 : double sweep;
82 :
83 : // Clockwise
84 30 : if(a1 > a2 && a2 > a3) {
85 15 : 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 15 : int ptcount = ceil(fabs(sweep/arcIncr));
104 :
105 15 : if(sweep < 0) arcIncr *= -1.0;
106 :
107 15 : double angle = a1;
108 :
109 3061 : for(int i = 0; i < ptcount - 1; i++) {
110 3046 : angle += arcIncr;
111 :
112 3046 : if(arcIncr > 0.0 && angle > PI) angle -= 2*PI;
113 3046 : if(arcIncr < 0.0 && angle < -1*PI) angle -= 2*PI;
114 :
115 3046 : double x = cx + r*cos(angle);
116 3046 : double y = cy + r*sin(angle);
117 :
118 3046 : line->addPoint(x, y, 0);
119 :
120 :
121 3046 : if((angle < a2) && ((angle + arcIncr) > a2)) {
122 0 : line->addPoint(ptOnArc);
123 : }
124 :
125 3046 : if((angle > a2) && ((angle + arcIncr) < a2)) {
126 15 : line->addPoint(ptOnArc);
127 : }
128 :
129 : }
130 15 : line->addPoint(ptEnd);
131 15 : delete center;
132 15 : }
133 :
|
