1 /*
2  * This file is part of ELKI:
3  * Environment for Developing KDD-Applications Supported by Index-Structures
4  *
5  * Copyright (C) 2018
6  * ELKI Development Team
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU Affero General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU Affero General Public License for more details.
17  *
18  * You should have received a copy of the GNU Affero General Public License
19  * along with this program. If not, see <http://www.gnu.org/licenses/>.
20  */
21 package de.lmu.ifi.dbs.elki.application.experiments;
22 
23 import java.awt.image.BufferedImage;
24 import java.io.File;
25 import java.io.IOException;
26 
27 import javax.imageio.ImageIO;
28 
29 import de.lmu.ifi.dbs.elki.application.AbstractApplication;
30 import de.lmu.ifi.dbs.elki.data.DoubleVector;
31 import de.lmu.ifi.dbs.elki.data.ModifiableHyperBoundingBox;
32 import de.lmu.ifi.dbs.elki.logging.Logging;
33 import de.lmu.ifi.dbs.elki.logging.progress.FiniteProgress;
34 import de.lmu.ifi.dbs.elki.math.geodesy.EarthModel;
35 import de.lmu.ifi.dbs.elki.math.geodesy.SphereUtil;
36 import de.lmu.ifi.dbs.elki.math.geodesy.SphericalVincentyEarthModel;
37 import de.lmu.ifi.dbs.elki.utilities.Alias;
38 import de.lmu.ifi.dbs.elki.utilities.documentation.Reference;
39 import de.lmu.ifi.dbs.elki.utilities.optionhandling.OptionID;
40 import de.lmu.ifi.dbs.elki.utilities.optionhandling.parameterization.Parameterization;
41 import de.lmu.ifi.dbs.elki.utilities.optionhandling.parameters.EnumParameter;
42 import de.lmu.ifi.dbs.elki.utilities.optionhandling.parameters.IntParameter;
43 import de.lmu.ifi.dbs.elki.utilities.optionhandling.parameters.ObjectParameter;
44 import net.jafama.FastMath;
45 
46 /**
47  * Visualization function for Cross-track, Along-track, and minimum distance
48  * function.
49  * <p>
50  * TODO: make origin point / rectangle configurable.
51  * <p>
52  * Reference:
53  * <p>
54  * Erich Schubert, Arthur Zimek, Hans-Peter Kriegel<br>
55  * Geodetic Distance Queries on R-Trees for Indexing Geographic Data<br>
56  * Int. Symp. Advances in Spatial and Temporal Databases (SSTD'2013)
57  *
58  * @author Niels Dörre
59  * @author Erich Schubert
60  * @since 0.5.5
61  */
62 @Alias({ "de.lmu.ifi.dbs.elki.application.geo.VisualizeGeodesicDistances" })
63 @Reference(authors = "Erich Schubert, Arthur Zimek, Hans-Peter Kriegel", //
64     title = "Geodetic Distance Queries on R-Trees for Indexing Geographic Data", //
65     booktitle = "Int. Symp. Advances in Spatial and Temporal Databases (SSTD'2013)", //
66     url = "https://doi.org/10.1007/978-3-642-40235-7_9", //
67     bibkey = "DBLP:conf/ssd/SchubertZK13")
68 public class VisualizeGeodesicDistances extends AbstractApplication {
69   /**
70    * Get a logger for this class.
71    */
72   private final static Logging LOG = Logging.getLogger(VisualizeGeodesicDistances.class);
73 
74   /**
75    * Visualization mode.
76    *
77    * @author Erich Schubert
78    */
79   public enum Mode {
80     /** Cross track distance */
81     XTD,
82     /** Along track distance */
83     ATD,
84     /** Mindist */
85     MINDIST
86   }
87 
88   /**
89    * Holds the file to print results to.
90    */
91   private File out;
92 
93   /**
94    * Image size.
95    */
96   protected int width = 2000, height = 1000;
97 
98   /**
99    * Number of steps for shades.
100    */
101   protected int steps = 10;
102 
103   /**
104    * Visualization mode.
105    */
106   protected Mode mode = Mode.XTD;
107 
108   /**
109    * Earth model.
110    */
111   protected EarthModel model;
112 
113   /**
114    * Constructor.
115    *
116    * @param out Output filename
117    * @param steps Number of steps in the color map
118    * @param mode Visualization mode
119    * @param model Earth model
120    */
VisualizeGeodesicDistances(File out, int resolution, int steps, Mode mode, EarthModel model)121   public VisualizeGeodesicDistances(File out, int resolution, int steps, Mode mode, EarthModel model) {
122     super();
123     this.width = resolution;
124     this.height = resolution >> 1;
125     this.out = out;
126     this.steps = steps;
127     this.mode = mode;
128     this.model = model;
129   }
130 
131   @Override
run()132   public void run() {
133     // Format: Latitude, Longitude
134     // München:
135     DoubleVector stap = DoubleVector.wrap(new double[] { 48.133333, 11.566667 });
136     // New York:
137     DoubleVector endp = DoubleVector.wrap(new double[] { 40.712778, -74.005833 });
138     // Bavaria:
139     ModifiableHyperBoundingBox bb = new ModifiableHyperBoundingBox(new double[] { 47.27011150, 8.97634970 }, new double[] { 50.56471420, 13.83963710 });
140     // Bavaria slice on lat
141     // bb = new ModifiableHyperBoundingBox(new double[] { 47.27011150, -80 }, //
142     // new double[] { 50.56471420, 80 });
143     // Bavaria slice on lon
144     // bb = new ModifiableHyperBoundingBox(new double[] { -10, 8.97634970 }, //
145     // new double[] { 50, 13.83963710 });
146 
147     BufferedImage img = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB);
148 
149     final double max = model.getEquatorialRadius() * Math.PI;
150 
151     // Red: left off-course, green: right off-course
152     int red = 0xffff0000;
153     int green = 0xff00ff00;
154 
155     FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("columns", width, LOG) : null;
156     for(int x = 0; x < width; x++) {
157       final double lon = x * 360. / width - 180.;
158       for(int y = 0; y < height; y++) {
159         final double lat = y * -180. / height + 90.;
160         switch(mode){
161         case ATD: {
162           final double atd = model.getEquatorialRadius() * SphereUtil.alongTrackDistanceDeg(stap.doubleValue(0), stap.doubleValue(1), endp.doubleValue(0), endp.doubleValue(1), lat, lon);
163           if(atd < 0) {
164             img.setRGB(x, y, colorMultiply(red, -atd / max, false));
165           }
166           else {
167             img.setRGB(x, y, colorMultiply(green, atd / max, false));
168           }
169           break;
170         }
171         case XTD: {
172           final double ctd = model.getEquatorialRadius() * SphereUtil.crossTrackDistanceDeg(stap.doubleValue(0), stap.doubleValue(1), endp.doubleValue(0), endp.doubleValue(1), lat, lon);
173           if(ctd < 0) {
174             img.setRGB(x, y, colorMultiply(red, -ctd / max, false));
175           }
176           else {
177             img.setRGB(x, y, colorMultiply(green, ctd / max, false));
178           }
179           break;
180         }
181         case MINDIST: {
182           final double dist = model.minDistDeg(lat, lon, bb.getMin(0), bb.getMin(1), bb.getMax(0), bb.getMax(1));
183           if(dist < 0) {
184             img.setRGB(x, y, colorMultiply(red, -dist / max, true));
185           }
186           else {
187             img.setRGB(x, y, colorMultiply(green, dist / max, true));
188           }
189           break;
190         }
191         }
192       }
193       LOG.incrementProcessed(prog);
194     }
195     LOG.ensureCompleted(prog);
196 
197     try {
198       ImageIO.write(img, "png", out);
199     }
200     catch(IOException e) {
201       LOG.exception(e);
202     }
203   }
204 
colorMultiply(int col, double reldist, boolean ceil)205   private int colorMultiply(int col, double reldist, boolean ceil) {
206     if(steps > 0) {
207       if(!ceil) {
208         reldist = FastMath.round(reldist * steps) / steps;
209       }
210       else {
211         reldist = FastMath.ceil(reldist * steps) / steps;
212       }
213     }
214     else if(steps < 0 && reldist > 0.) {
215       double s = reldist * -steps;
216       double off = Math.abs(s - FastMath.round(s));
217       double factor = -steps * 1. / 1000; // height;
218       if(off < factor) { // Blend with black:
219         factor = (off / factor);
220         int a = (col >> 24) & 0xFF;
221         a = (int) (a * FastMath.sqrt(reldist)) & 0xFF;
222         a = (int) ((1 - factor) * 0xFF + factor * a);
223         int r = (int) (factor * ((col >> 16) & 0xFF));
224         int g = (int) (factor * ((col >> 8) & 0xFF));
225         int b = (int) (factor * (col & 0xFF));
226         return a << 24 | r << 16 | g << 8 | b;
227       }
228     }
229     int a = (col >> 24) & 0xFF, r = (col >> 16) & 0xFF, g = (col >> 8) & 0xFF,
230         b = (col) & 0xFF;
231     a = (int) (a * FastMath.sqrt(reldist)) & 0xFF;
232     return a << 24 | r << 16 | g << 8 | b;
233   }
234 
235   /**
236    * Main method for application.
237    *
238    * @param args Parameters
239    */
main(String[] args)240   public static void main(String[] args) {
241     runCLIApplication(VisualizeGeodesicDistances.class, args);
242   }
243 
244   /**
245    * Parameterization class.
246    *
247    * @author Erich Schubert
248    */
249   public static class Parameterizer extends AbstractApplication.Parameterizer {
250     /**
251      * Number of steps in the distance map.
252      */
253     public static final OptionID STEPS_ID = new OptionID("geodistvis.steps", "Number of steps for the distance map. Use negative numbers to get contour lines.");
254 
255     /**
256      * Image resolution.
257      */
258     public static final OptionID RESOLUTION_ID = new OptionID("geodistvis.resolution", "Horizontal resolution for the image map (vertical resolution is horizonal / 2).");
259 
260     /**
261      * Visualization mode.
262      */
263     public static final OptionID MODE_ID = new OptionID("geodistvis.mode", "Visualization mode.");
264 
265     /**
266      * Holds the file to print results to.
267      */
268     protected File out = null;
269 
270     /**
271      * Number of steps in the color map.
272      */
273     protected int steps = 0;
274 
275     /**
276      * Horizontal resolution.
277      */
278     protected int resolution = 2000;
279 
280     /**
281      * Visualization mode.
282      */
283     protected Mode mode = Mode.XTD;
284 
285     /**
286      * Earth model to use.
287      */
288     protected EarthModel model;
289 
290     @Override
makeOptions(Parameterization config)291     protected void makeOptions(Parameterization config) {
292       super.makeOptions(config);
293       out = super.getParameterOutputFile(config, "Output image file name.");
294       IntParameter stepsP = new IntParameter(STEPS_ID) //
295           .setOptional(true);
296       if(config.grab(stepsP)) {
297         steps = stepsP.intValue();
298       }
299       IntParameter resolutionP = new IntParameter(RESOLUTION_ID, 2000);
300       if(config.grab(resolutionP)) {
301         resolution = resolutionP.intValue();
302       }
303       EnumParameter<Mode> modeP = new EnumParameter<>(MODE_ID, Mode.class, Mode.XTD);
304       if(config.grab(modeP)) {
305         mode = modeP.getValue();
306       }
307       ObjectParameter<EarthModel> modelP = new ObjectParameter<>(EarthModel.MODEL_ID, EarthModel.class, SphericalVincentyEarthModel.class);
308       if(config.grab(modelP)) {
309         model = modelP.instantiateClass(config);
310       }
311     }
312 
313     @Override
makeInstance()314     protected VisualizeGeodesicDistances makeInstance() {
315       return new VisualizeGeodesicDistances(out, resolution, steps, mode, model);
316     }
317   }
318 }
319