1 | package UtilsPlugin;
|
---|
2 |
|
---|
3 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
4 |
|
---|
5 | import java.awt.event.ActionEvent;
|
---|
6 | import java.util.ArrayList;
|
---|
7 | import java.util.Collection;
|
---|
8 | import java.util.Collections;
|
---|
9 | import java.util.HashSet;
|
---|
10 | import java.util.LinkedList;
|
---|
11 | import java.util.List;
|
---|
12 |
|
---|
13 | import org.openstreetmap.josm.Main;
|
---|
14 | import org.openstreetmap.josm.command.ChangeCommand;
|
---|
15 | import org.openstreetmap.josm.command.Command;
|
---|
16 | import org.openstreetmap.josm.command.DeleteCommand;
|
---|
17 | import org.openstreetmap.josm.command.SequenceCommand;
|
---|
18 | import org.openstreetmap.josm.data.coor.LatLon;
|
---|
19 | import org.openstreetmap.josm.data.osm.Node;
|
---|
20 | import org.openstreetmap.josm.data.osm.OsmPrimitive;
|
---|
21 | import org.openstreetmap.josm.data.osm.Way;
|
---|
22 | import org.openstreetmap.josm.data.osm.visitor.CollectBackReferencesVisitor;
|
---|
23 |
|
---|
24 | import org.openstreetmap.josm.data.osm.DataSet;
|
---|
25 | import org.openstreetmap.josm.actions.JosmAction;
|
---|
26 |
|
---|
27 | public class SimplifyWayAction extends JosmAction {
|
---|
28 | public SimplifyWayAction() {
|
---|
29 | super(tr("Simplify Way"), "simplify",
|
---|
30 | tr("Delete unnecessary nodes from a way."), 0, 0, true);
|
---|
31 | }
|
---|
32 |
|
---|
33 | public void actionPerformed(ActionEvent e) {
|
---|
34 | Collection<OsmPrimitive> selection = Main.ds.getSelected();
|
---|
35 |
|
---|
36 | if (selection.size() == 1 && selection.iterator().next() instanceof Way) {
|
---|
37 | simplifyWay((Way) selection.iterator().next());
|
---|
38 | }
|
---|
39 | }
|
---|
40 |
|
---|
41 | public void simplifyWay(Way w) {
|
---|
42 | double threshold = Double.parseDouble(
|
---|
43 | Main.pref.get("simplify-way.max-error", "50"));
|
---|
44 |
|
---|
45 | Way wnew = new Way(w);
|
---|
46 |
|
---|
47 | int toI = wnew.nodes.size() - 1;
|
---|
48 | for (int i = wnew.nodes.size() - 1; i >= 0; i--) {
|
---|
49 | CollectBackReferencesVisitor backRefsV =
|
---|
50 | new CollectBackReferencesVisitor(Main.ds, false);
|
---|
51 | backRefsV.visit(wnew.nodes.get(i));
|
---|
52 | boolean used = false;
|
---|
53 | if (backRefsV.data.size() == 1) {
|
---|
54 | used = Collections.frequency(
|
---|
55 | w.nodes, wnew.nodes.get(i)) > 1;
|
---|
56 | } else {
|
---|
57 | backRefsV.data.remove(w);
|
---|
58 | used = !backRefsV.data.isEmpty();
|
---|
59 | }
|
---|
60 |
|
---|
61 | if (used) {
|
---|
62 | if (toI - i >= 2) {
|
---|
63 | ArrayList<Node> ns = new ArrayList<Node>();
|
---|
64 | simplifyWayRange(wnew, i, toI, ns, threshold);
|
---|
65 | for (int j = toI-1; j > i; j--) wnew.nodes.remove(j);
|
---|
66 | wnew.nodes.addAll(i+1, ns);
|
---|
67 | }
|
---|
68 | toI = i;
|
---|
69 | }
|
---|
70 | }
|
---|
71 |
|
---|
72 | HashSet<Node> delNodes = new HashSet<Node>();
|
---|
73 | delNodes.addAll(w.nodes);
|
---|
74 | delNodes.removeAll(wnew.nodes);
|
---|
75 |
|
---|
76 | if (wnew.nodes.size() != w.nodes.size()) {
|
---|
77 | Collection<Command> cmds = new LinkedList<Command>();
|
---|
78 | cmds.add(new ChangeCommand(w, wnew));
|
---|
79 | cmds.add(new DeleteCommand(delNodes));
|
---|
80 | Main.main.undoRedo.add(
|
---|
81 | new SequenceCommand(tr("Simplify Way (remove {0} nodes)",
|
---|
82 | delNodes.size()),
|
---|
83 | cmds));
|
---|
84 | Main.map.repaint();
|
---|
85 | }
|
---|
86 | }
|
---|
87 |
|
---|
88 | /*
|
---|
89 | * Takes an interval [from,to] and adds nodes from the set (from,to) to
|
---|
90 | * ns.
|
---|
91 | */
|
---|
92 | public void simplifyWayRange(Way wnew, int from, int to, ArrayList<Node> ns, double thr) {
|
---|
93 | Node fromN = wnew.nodes.get(from), toN = wnew.nodes.get(to);
|
---|
94 |
|
---|
95 | int imax = -1;
|
---|
96 | double xtemax = 0;
|
---|
97 | for (int i = from+1; i < to; i++) {
|
---|
98 | Node n = wnew.nodes.get(i);
|
---|
99 | double xte = radtometers(linedist(
|
---|
100 | fromN.coor.lat(), fromN.coor.lon(),
|
---|
101 | n.coor.lat(), n.coor.lon(),
|
---|
102 | toN.coor.lat(), toN.coor.lon()));
|
---|
103 | if (xte > xtemax) {
|
---|
104 | xtemax = xte;
|
---|
105 | imax = i;
|
---|
106 | }
|
---|
107 | }
|
---|
108 |
|
---|
109 | if (imax != -1 && xtemax >= thr) {
|
---|
110 | simplifyWayRange(wnew, from, imax, ns, thr);
|
---|
111 | ns.add(wnew.nodes.get(imax));
|
---|
112 | simplifyWayRange(wnew, imax, to, ns, thr);
|
---|
113 | }
|
---|
114 | }
|
---|
115 |
|
---|
116 | /* ----------------------------------------------------------------------
|
---|
117 | * Everything below this comment was converted from C to Java by Frederik
|
---|
118 | * Ramm. The original sources are the files grtcirc.c and smplrout.c from
|
---|
119 | * the gpsbabel source code (www.gpsbabel.org), which is under GPL. The
|
---|
120 | * relevant code portions have been written by Robert Lipe.
|
---|
121 | *
|
---|
122 | * Method names have been left unchanged where possible.
|
---|
123 | */
|
---|
124 |
|
---|
125 | public static double EARTH_RAD = 6378137.0;
|
---|
126 | public static double radmiles = EARTH_RAD*100.0/2.54/12.0/5280.0;
|
---|
127 |
|
---|
128 | public static double[] crossproduct(double[] v1, double[] v2) {
|
---|
129 | double[] rv = new double[3];
|
---|
130 | rv[0] = v1[1]*v2[2]-v2[1]*v1[2];
|
---|
131 | rv[1] = v1[2]*v2[0]-v2[2]*v1[0];
|
---|
132 | rv[2] = v1[0]*v2[1]-v1[1]*v2[0];
|
---|
133 | return rv;
|
---|
134 | }
|
---|
135 |
|
---|
136 | public static double dotproduct(double[] v1, double[] v2) {
|
---|
137 | return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
|
---|
138 | }
|
---|
139 |
|
---|
140 | public static double radtomiles(double rads) {
|
---|
141 | return (rads*radmiles);
|
---|
142 | }
|
---|
143 |
|
---|
144 | public static double radtometers(double rads) {
|
---|
145 | return (rads * EARTH_RAD);
|
---|
146 | }
|
---|
147 |
|
---|
148 | public static double veclen(double[] vec) {
|
---|
149 | return Math.sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
|
---|
150 | }
|
---|
151 |
|
---|
152 | public static double gcdist(double lat1, double lon1, double lat2, double lon2)
|
---|
153 | {
|
---|
154 | double res;
|
---|
155 | double sdlat, sdlon;
|
---|
156 |
|
---|
157 | sdlat = Math.sin((lat1 - lat2) / 2.0);
|
---|
158 | sdlon = Math.sin((lon1 - lon2) / 2.0);
|
---|
159 |
|
---|
160 | res = Math.sqrt(sdlat * sdlat + Math.cos(lat1) * Math.cos(lat2) * sdlon * sdlon);
|
---|
161 |
|
---|
162 | if (res > 1.0) {
|
---|
163 | res = 1.0;
|
---|
164 | } else if (res < -1.0) {
|
---|
165 | res = -1.0;
|
---|
166 | }
|
---|
167 |
|
---|
168 | res = Math.asin(res);
|
---|
169 | return 2.0 * res;
|
---|
170 | }
|
---|
171 |
|
---|
172 | static double linedist(double lat1, double lon1, double lat2, double lon2, double lat3, double lon3) {
|
---|
173 |
|
---|
174 | double dot;
|
---|
175 |
|
---|
176 | /* degrees to radians */
|
---|
177 | lat1 = Math.toRadians(lat1); lon1 = Math.toRadians(lon1);
|
---|
178 | lat2 = Math.toRadians(lat2); lon2 = Math.toRadians(lon2);
|
---|
179 | lat3 = Math.toRadians(lat3); lon3 = Math.toRadians(lon3);
|
---|
180 |
|
---|
181 | /* polar to ECEF rectangular */
|
---|
182 | double[] v1 = new double[3];
|
---|
183 | double[] v2 = new double[3];
|
---|
184 | double[] v3 = new double[3];
|
---|
185 | v1[0] = Math.cos(lon1)*Math.cos(lat1); v1[1] = Math.sin(lat1); v1[2] = Math.sin(lon1)*Math.cos(lat1);
|
---|
186 | v2[0] = Math.cos(lon2)*Math.cos(lat2); v2[1] = Math.sin(lat2); v2[2] = Math.sin(lon2)*Math.cos(lat2);
|
---|
187 | v3[0] = Math.cos(lon3)*Math.cos(lat3); v3[1] = Math.sin(lat3); v3[2] = Math.sin(lon3)*Math.cos(lat3);
|
---|
188 |
|
---|
189 | /* 'va' is the axis; the line that passes through the center of the earth
|
---|
190 | * and is perpendicular to the great circle through point 1 and point 2
|
---|
191 | * It is computed by taking the cross product of the '1' and '2' vectors.*/
|
---|
192 | double[] va = crossproduct(v1, v2);
|
---|
193 | double la = veclen(va);
|
---|
194 |
|
---|
195 | if (la != 0) {
|
---|
196 | va[0] /= la;
|
---|
197 | va[1] /= la;
|
---|
198 | va[2] /= la;
|
---|
199 |
|
---|
200 | /* dot is the component of the length of '3' that is along the axis.
|
---|
201 | * What's left is a non-normalized vector that lies in the plane of
|
---|
202 | * 1 and 2. */
|
---|
203 |
|
---|
204 | dot = dotproduct(v3, va);
|
---|
205 |
|
---|
206 | double[] vp = new double[3];
|
---|
207 | vp[0]=v3[0]-dot*va[0];
|
---|
208 | vp[1]=v3[1]-dot*va[1];
|
---|
209 | vp[2]=v3[2]-dot*va[2];
|
---|
210 |
|
---|
211 | double lp = veclen(vp);
|
---|
212 |
|
---|
213 | if (lp != 0) {
|
---|
214 |
|
---|
215 | /* After this, 'p' is normalized */
|
---|
216 | vp[0] /= lp;
|
---|
217 | vp[1] /= lp;
|
---|
218 | vp[2] /= lp;
|
---|
219 |
|
---|
220 | double[] cp1 = crossproduct(v1, vp);
|
---|
221 | double dp1 = dotproduct(cp1, va);
|
---|
222 |
|
---|
223 | double[] cp2 = crossproduct(v2, vp);
|
---|
224 | double dp2 = dotproduct(cp2, va);
|
---|
225 |
|
---|
226 | if ( dp1 >= 0 && dp2 >= 0 ) {
|
---|
227 | /* rather than call gcdist and all its sines and cosines and
|
---|
228 | * worse, we can get the angle directly. It's the arctangent
|
---|
229 | * of the length of the component of vector 3 along the axis
|
---|
230 | * divided by the length of the component of vector 3 in the
|
---|
231 | * plane. We already have both of those numbers.
|
---|
232 | *
|
---|
233 | * atan2 would be overkill because lp and Math.abs are both
|
---|
234 | * known to be positive. */
|
---|
235 | return Math.atan(Math.abs(dot)/lp);
|
---|
236 | }
|
---|
237 |
|
---|
238 | /* otherwise, get the distance from the closest endpoint */
|
---|
239 | double c1 = dotproduct(v1, vp);
|
---|
240 | double c2 = dotproduct(v2, vp);
|
---|
241 | dp1 = Math.abs(dp1);
|
---|
242 | dp2 = Math.abs(dp2);
|
---|
243 |
|
---|
244 | /* This is a hack. d$n$ is proportional to the sine of the angle
|
---|
245 | * between point $n$ and point p. That preserves orderedness up
|
---|
246 | * to an angle of 90 degrees. c$n$ is proportional to the cosine
|
---|
247 | * of the same angle; if the angle is over 90 degrees, c$n$ is
|
---|
248 | * negative. In that case, we flop the sine across the y=1 axis
|
---|
249 | * so that the resulting value increases as the angle increases.
|
---|
250 | *
|
---|
251 | * This only works because all of the points are on a unit sphere. */
|
---|
252 |
|
---|
253 | if (c1 < 0) {
|
---|
254 | dp1 = 2 - dp1;
|
---|
255 | }
|
---|
256 | if (c2 < 0) {
|
---|
257 | dp2 = 2 - dp2;
|
---|
258 | }
|
---|
259 |
|
---|
260 | if (Math.abs(dp1) < Math.abs(dp2)) {
|
---|
261 | return gcdist(lat1,lon1,lat3,lon3);
|
---|
262 | } else {
|
---|
263 | return gcdist(lat2,lon2,lat3,lon3);
|
---|
264 | }
|
---|
265 | } else {
|
---|
266 | /* lp is 0 when 3 is 90 degrees from the great circle */
|
---|
267 | return Math.PI/2;
|
---|
268 | }
|
---|
269 | } else {
|
---|
270 | /* la is 0 when 1 and 2 are either the same point or 180 degrees apart */
|
---|
271 | dot = dotproduct(v1, v2);
|
---|
272 | if (dot >= 0) {
|
---|
273 | return gcdist(lat1,lon1,lat3,lon3);
|
---|
274 | } else {
|
---|
275 | return 0;
|
---|
276 | }
|
---|
277 | }
|
---|
278 | }
|
---|
279 | }
|
---|