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