/layers-store/src/main/java/org/ala/layers/intersect/ComplexRegion.java

http://alageospatialportal.googlecode.com/ · Java · 637 lines · 412 code · 68 blank · 157 comment · 197 complexity · f42f680534029812dd6003067c7f80bf MD5 · raw file

  1. /**************************************************************************
  2. * Copyright (C) 2010 Atlas of Living Australia
  3. * All Rights Reserved.
  4. *
  5. * The contents of this file are subject to the Mozilla Public
  6. * License Version 1.1 (the "License"); you may not use this file
  7. * except in compliance with the License. You may obtain a copy of
  8. * the License at http://www.mozilla.org/MPL/
  9. *
  10. * Software distributed under the License is distributed on an "AS
  11. * IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
  12. * implied. See the License for the specific language governing
  13. * rights and limitations under the License.
  14. ***************************************************************************/
  15. package org.ala.layers.intersect;
  16. import java.util.ArrayList;
  17. /**
  18. * ComplexRegion is a collection of SimpleRegion, expect POLYGONs for now.
  19. *
  20. * treat as a shape file, overlapping regions cancel out presence.
  21. *
  22. * TODO: clockwise/anticlockwise identification
  23. *
  24. * @author Adam Collins
  25. */
  26. public class ComplexRegion extends SimpleRegion {
  27. static public SimpleRegion parseComplexRegion(String[] polygons) {
  28. ComplexRegion cr = new ComplexRegion();
  29. for (String s : polygons) {
  30. cr.addPolygon(SimpleRegion.parseSimpleRegion(s));
  31. }
  32. cr.useMask(-1, -1, -1);
  33. return cr;
  34. }
  35. /**
  36. * list of SimpleRegion members
  37. */
  38. ArrayList<SimpleRegion> simpleregions;
  39. /**
  40. * bounding box for all, see SimpleRegion boundingbox.
  41. */
  42. double[][] boundingbox_all;
  43. /**
  44. * value assigned
  45. */
  46. int value;
  47. /**
  48. * array for speeding up isWithin
  49. */
  50. byte[][] mask;
  51. Object[][] maskDepth;
  52. /**
  53. * mask height
  54. */
  55. int mask_height;
  56. /**
  57. * mask width
  58. */
  59. int mask_width;
  60. /**
  61. * mask multiplier for longitude inputs
  62. */
  63. double mask_long_multiplier;
  64. /**
  65. * mask mulitplier for latitude inputs
  66. */
  67. double mask_lat_multiplier;
  68. /**
  69. * maintain mapping for simpleregions belonging to the same polygon
  70. */
  71. ArrayList<Integer> polygons;
  72. /**
  73. * Constructor for empty ComplexRegion
  74. */
  75. public ComplexRegion() {
  76. super();
  77. simpleregions = new ArrayList();
  78. boundingbox_all = new double[2][2];
  79. value = -1;
  80. mask = null;
  81. polygons = new ArrayList();
  82. }
  83. /**
  84. * sets integer value stored
  85. * @param value_ as int
  86. */
  87. public void setValue(int value_) {
  88. value = value_;
  89. }
  90. /**
  91. * gets integer value stored
  92. * @return int
  93. */
  94. public int getValue() {
  95. return value;
  96. }
  97. /**
  98. * gets the bounding box for shapes in this ComplexRegion
  99. *
  100. * @return bounding box for ComplexRegion as double [][]
  101. */
  102. @Override
  103. public double[][] getBoundingBox() {
  104. return boundingbox_all;
  105. }
  106. /**
  107. * adds a new polygon
  108. *
  109. * note: if a mask is in use must call <code>useMask</code> again
  110. * @param points_ points = double[n][2]
  111. * where
  112. * n is number of points
  113. * [][0] is longitude
  114. * [][1] is latitude
  115. */
  116. public void addPolygon(SimpleRegion sr) {
  117. simpleregions.add(sr);
  118. /* update boundingbox_all */
  119. double[][] bb = sr.getBoundingBox();
  120. if (simpleregions.size() == 1 || boundingbox_all[0][0] > bb[0][0]) {
  121. boundingbox_all[0][0] = bb[0][0];
  122. }
  123. if (simpleregions.size() == 1 || boundingbox_all[1][0] < bb[1][0]) {
  124. boundingbox_all[1][0] = bb[1][0];
  125. }
  126. if (simpleregions.size() == 1 || boundingbox_all[0][1] > bb[0][1]) {
  127. boundingbox_all[0][1] = bb[0][1];
  128. }
  129. if (simpleregions.size() == 1 || boundingbox_all[1][1] < bb[1][1]) {
  130. boundingbox_all[1][1] = bb[1][1];
  131. }
  132. bounding_box = boundingbox_all;
  133. }
  134. /**
  135. * returns true when the point provided is within the ComplexRegion
  136. *
  137. * uses <code>mask</code> when available
  138. *
  139. * note: type UNDEFINED implies no boundary, always returns true.
  140. *
  141. * @param longitude
  142. * @param latitude
  143. * @return true iff point is within or on the edge of this ComplexRegion
  144. */
  145. @Override
  146. public boolean isWithin(double longitude, double latitude) {
  147. if (simpleregions.size() == 1) {
  148. return simpleregions.get(0).isWithin(longitude, latitude);
  149. }
  150. if (boundingbox_all[0][0] > longitude || boundingbox_all[1][0] < longitude
  151. || boundingbox_all[0][1] > latitude || boundingbox_all[1][1] < latitude) {
  152. return false;
  153. }
  154. short[] countsIn = new short[polygons.get(polygons.size()-1) + 1];
  155. //int count_in = 0; //count of regions overlapping the point
  156. if (mask != null) {
  157. /* use mask if exists */
  158. int long1 = (int) Math.floor((longitude - boundingbox_all[0][0]) * mask_long_multiplier);
  159. int lat1 = (int) Math.floor((latitude - boundingbox_all[0][1]) * mask_lat_multiplier);
  160. if (long1 == mask[0].length) {
  161. long1--;
  162. }
  163. if (lat1 == mask.length) {
  164. lat1--;
  165. }
  166. if (mask[lat1][long1] == SimpleRegion.GI_FULLY_PRESENT) {
  167. return true;
  168. } else if (mask[lat1][long1] == SimpleRegion.GI_UNDEFINED
  169. || mask[lat1][long1] == SimpleRegion.GI_ABSENCE) {
  170. return false;
  171. }
  172. //partial, try maskDepth and sum overlaps
  173. if (maskDepth != null && maskDepth[lat1][long1] != null) {
  174. int[] d = (int[]) maskDepth[lat1][long1];
  175. for (int i = 0; i < d.length; i++) {
  176. if (simpleregions.get(d[i]).isWithin(longitude, latitude)) {
  177. countsIn[polygons.get(d[i])]++;
  178. }
  179. }
  180. /* true iif within an odd number of regions for any polygon*/
  181. for (int i = 0; i < countsIn.length; i++) {
  182. if (countsIn[i] % 2 == 1) {
  183. return true;
  184. }
  185. }
  186. }
  187. }
  188. /* check for all SimpleRegions */
  189. for (int i = 0; i < simpleregions.size(); i++) {
  190. if (simpleregions.get(i).isWithin(longitude, latitude)) {
  191. countsIn[polygons.get(i)]++;
  192. }
  193. }
  194. /* true iif within an odd number of regions for any polygon*/
  195. for (int i = 0; i < countsIn.length; i++) {
  196. if (countsIn[i] % 2 == 1) {
  197. return true;
  198. }
  199. }
  200. return false;
  201. }
  202. /**
  203. * builds a grid (mask) to speed up isWithin.
  204. *
  205. * TODO: split out shapes with large numbers of points in GI_PARTIALLY_PRESENT grid cells.
  206. *
  207. * TODO: automatic(best) height/width specification
  208. *
  209. * @param width
  210. * @param height
  211. */
  212. public void useMask(int width, int height, int depthThreashold) {
  213. //calculate defaults for -1 inputs
  214. double[][] bb = getBoundingBox();
  215. int length = 0;
  216. for (SimpleRegion sr : simpleregions) {
  217. length += sr.getNumberOfPoints();
  218. }
  219. int w = (int) ((bb[1][0] - bb[0][0]) * 3);
  220. int h = (int) ((bb[1][1] - bb[0][1]) * 3);
  221. if (length > 5000) {
  222. w = 200;
  223. h = 200;
  224. }
  225. if (w > 200) {
  226. w = 200;
  227. }
  228. if (h > 200) {
  229. h = 200;
  230. }
  231. if (width == -1) {
  232. width = w;
  233. }
  234. if (height == -1) {
  235. height = h;
  236. }
  237. if (depthThreashold == -1) {
  238. depthThreashold = 100;
  239. }
  240. if(width < 3 || height < 3) {
  241. return;
  242. }
  243. int i, j;
  244. /* class variables assignment */
  245. mask_width = width;
  246. mask_height = height;
  247. mask_long_multiplier =
  248. mask_width / (double) (boundingbox_all[1][0] - boundingbox_all[0][0]);
  249. mask_lat_multiplier =
  250. mask_height / (double) (boundingbox_all[1][1] - boundingbox_all[0][1]);
  251. /* end result mask */
  252. mask = new byte[height][width];
  253. ArrayList<Integer>[][] md = null;
  254. if (simpleregions.size() > depthThreashold) {
  255. //use mask depth as well
  256. md = new ArrayList[height][width];
  257. }
  258. /* temp mask for current Polygon */
  259. byte[][] shapemask = new byte[height][width];
  260. /* temp mask for current SimpleRegion */
  261. byte[][] shapemaskregion = new byte[height][width];
  262. int k = 0;
  263. while(k < simpleregions.size()) {
  264. int p = k;
  265. for (; k < simpleregions.size()
  266. && (p == k || polygons.get(k - 1) == polygons.get(k)); k++) {
  267. SimpleRegion sr = simpleregions.get(k);
  268. sr.getOverlapGridCells(boundingbox_all[0][0], boundingbox_all[0][1], boundingbox_all[1][0], boundingbox_all[1][1], width, height, shapemaskregion);
  269. //shapemaskregion into shapemask
  270. for (i = 0; i < height; i++) {
  271. for (j = 0; j < width; j++) {
  272. if (shapemaskregion[i][j] == SimpleRegion.GI_PARTIALLY_PRESENT
  273. || shapemask[i][j] ==SimpleRegion.GI_PARTIALLY_PRESENT) {
  274. shapemask[i][j] = SimpleRegion.GI_PARTIALLY_PRESENT; //partially inside
  275. if (md != null) {
  276. if (md[i][j] == null) {
  277. md[i][j] = new ArrayList<Integer>();
  278. }
  279. md[i][j].add(k);
  280. }
  281. } else if (shapemaskregion[i][j] == SimpleRegion.GI_FULLY_PRESENT ) {
  282. if (shapemask[i][j] == SimpleRegion.GI_FULLY_PRESENT ) {
  283. shapemask[i][j] = SimpleRegion.GI_ABSENCE; //completely inside
  284. } else {
  285. shapemask[i][j] = SimpleRegion.GI_FULLY_PRESENT ;//completely outside (inside of a cutout region)
  286. }
  287. if (md != null) {
  288. if (md[i][j] == null) {
  289. md[i][j] = new ArrayList<Integer>();
  290. }
  291. md[i][j].add(k);
  292. }
  293. }
  294. /* reset shapemaskregion for next part */
  295. shapemaskregion[i][j] = 0;
  296. }
  297. }
  298. }
  299. //shapemask into mask
  300. for (i = 0; i < height; i++) {
  301. for (j = 0; j < width; j++) {
  302. if (shapemask[i][j] == SimpleRegion.GI_FULLY_PRESENT
  303. || mask[i][j] == SimpleRegion.GI_FULLY_PRESENT) {
  304. mask[i][j] = SimpleRegion.GI_FULLY_PRESENT; //partially inside
  305. } else if (shapemask[i][j] == SimpleRegion.GI_PARTIALLY_PRESENT) {
  306. mask[i][j] = SimpleRegion.GI_PARTIALLY_PRESENT;
  307. }
  308. /* reset shapemask for next part */
  309. shapemask[i][j] = 0;
  310. }
  311. }
  312. }
  313. //maskDepth to int[]
  314. if (md != null) {
  315. maskDepth = new Object[md.length][md[0].length];
  316. for (i = 0; i < height; i++) {
  317. for (j = 0; j < width; j++) {
  318. if (md[i][j] != null && mask[i][j] == SimpleRegion.GI_PARTIALLY_PRESENT ) {
  319. int[] d = new int[md[i][j].size()];
  320. for (k = 0; k < d.length; k++) {
  321. d[k] = md[i][j].get(k);
  322. }
  323. maskDepth[i][j] = d;
  324. }
  325. }
  326. }
  327. }
  328. }
  329. /**
  330. * determines overlap with a grid
  331. *
  332. * for type POLYGON
  333. * when <code>three_state_map</code> is not null populate it with one of:
  334. * GI_UNDEFINED
  335. * GI_PARTIALLY_PRESENT
  336. * GI_FULLY_PRESENT
  337. * GI_ABSENCE
  338. *
  339. * @param longitude1
  340. * @param latitude1
  341. * @param longitude2
  342. * @param latitude2
  343. * @param xres number of longitude segements as int
  344. * @param yres number of latitude segments as int
  345. * @return (x,y) as double [][2] for each grid cell at least partially falling
  346. * within the specified region of the specified resolution beginning at 0,0
  347. * for minimum longitude and latitude through to xres,yres for maximums *
  348. */
  349. @Override
  350. public int[][] getOverlapGridCells(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) {
  351. int i,j;
  352. // if no threestate map exists, create one
  353. if (three_state_map == null) {
  354. three_state_map = new byte[height][width];
  355. }
  356. byte[][] tmpMask = new byte[height][width];
  357. int k = 0;
  358. while(k < simpleregions.size()) {
  359. int p = k;
  360. for (; k < simpleregions.size()
  361. && (p == k || polygons.get(k - 1).equals(polygons.get(k))); k++) {
  362. SimpleRegion sr = simpleregions.get(k);
  363. sr.getOverlapGridCells_Acc(longitude1, latitude1, longitude2, latitude2, width, height, tmpMask);
  364. }
  365. fillAccMask(p, k-1, longitude1, latitude1, longitude2, latitude2, width, height, tmpMask, true);
  366. //tmpMask into mask
  367. for (i = 0; i < height; i++) {
  368. for (j = 0; j < width; j++) {
  369. if(tmpMask[i][j] == 2 || three_state_map[i][j] == 2) {
  370. three_state_map[i][j] = 2;
  371. } else if(tmpMask[i][j] == 1) {
  372. three_state_map[i][j] = 1;
  373. }
  374. /* reset shapemask for next part */
  375. tmpMask[i][j] = 0;
  376. }
  377. }
  378. }
  379. boolean cellsReturned = !noCellsReturned;
  380. if (cellsReturned) {
  381. int [][] data = new int[width * height][2];
  382. int p = 0;
  383. for(i=0;i<height;i++) {
  384. for(j=0;j<width;j++) {
  385. if (three_state_map[i][j] != GI_UNDEFINED) { //undefined == absence
  386. data[p][0] = j;
  387. data[p][1] = i;
  388. p++;
  389. }
  390. }
  391. }
  392. data = java.util.Arrays.copyOf(data, p);
  393. }
  394. return null;
  395. }
  396. @Override
  397. public boolean isWithin_EPSG900913(double longitude, double latitude) {
  398. short[] countsIn = new short[polygons.get(polygons.size()-1) + 1];
  399. /* check for all SimpleRegions */
  400. for (int i = 0; i < simpleregions.size(); i++) {
  401. if (simpleregions.get(i).isWithin_EPSG900913(longitude, latitude)) {
  402. countsIn[polygons.get(i)]++;
  403. }
  404. }
  405. /* true iif within an odd number of regions for any polygon*/
  406. for (int i = 0; i < countsIn.length; i++) {
  407. if (countsIn[i] % 2 == 1) {
  408. return true;
  409. }
  410. }
  411. return false;
  412. }
  413. @Override
  414. public int[][] getOverlapGridCells_EPSG900913(double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) {
  415. int i, j;
  416. int[][] output = null;
  417. byte[][] mask = three_state_map;
  418. // if no threestate map exists, create one
  419. if (mask == null) {
  420. mask = new byte[height][width];
  421. three_state_map = mask;
  422. }
  423. byte[][] tmpMask = new byte[height][width];
  424. int k = 0;
  425. while(k < simpleregions.size()) {
  426. int p = k;
  427. for (; k < simpleregions.size()
  428. && (p == k || polygons.get(k - 1).equals(polygons.get(k))); k++) {
  429. SimpleRegion sr = simpleregions.get(k);
  430. sr.getOverlapGridCells_Acc_EPSG900913(longitude1, latitude1, longitude2, latitude2, width, height, tmpMask);
  431. }
  432. fillAccMask_EPSG900913(k, p-1, longitude1, latitude1, longitude2, latitude2, width, height, tmpMask, true);
  433. //tmpMask into mask
  434. for (i = 0; i < height; i++) {
  435. for (j = 0; j < width; j++) {
  436. if(tmpMask[i][j] == 2 || mask[i][j] == 2) {
  437. mask[i][j] = 2;
  438. } else if(tmpMask[i][j] == 1) {
  439. mask[i][j] = 1;
  440. }
  441. /* reset shapemask for next part */
  442. tmpMask[i][j] = 0;
  443. }
  444. }
  445. }
  446. boolean cellsReturned = !noCellsReturned;
  447. if (cellsReturned) {
  448. int [][] data = new int[width * height][2];
  449. int p = 0;
  450. for(i=0;i<height;i++) {
  451. for(j=0;j<width;j++) {
  452. if (mask[i][j] != GI_UNDEFINED) { //undefined == absence
  453. data[p][0] = j;
  454. data[p][1] = i;
  455. p++;
  456. }
  457. }
  458. }
  459. data = java.util.Arrays.copyOf(data, p);
  460. return data;
  461. }
  462. return null;
  463. }
  464. void addSet(ArrayList<SimpleRegion> simpleRegions) {
  465. int nextSetNumber = (polygons.size() > 0) ? polygons.get(polygons.size() - 1) + 1 : 0;
  466. for (int i = 0; i < simpleRegions.size(); i++) {
  467. addPolygon(simpleRegions.get(i));
  468. polygons.add(nextSetNumber);
  469. }
  470. }
  471. public int[][] fillAccMask(int startPolygon, int endPolygon, double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) {
  472. double divx = (longitude2 - longitude1) / width;
  473. double divy = (latitude2 - latitude1) / height;
  474. int i, j;
  475. //do raster check
  476. int[][] data = null;
  477. boolean cellsReturned = !noCellsReturned;
  478. if (cellsReturned) {
  479. data = new int[width * height][2];
  480. }
  481. int p = 0;
  482. for (j = 0; j < three_state_map[0].length; j++) {
  483. for (i = 0; i < three_state_map.length; i++) {
  484. if (three_state_map[i][j] == GI_PARTIALLY_PRESENT) {
  485. //if it is partially present, do nothing
  486. } else if ((j == 0 || three_state_map[i][j - 1] == GI_PARTIALLY_PRESENT)) {
  487. if (i > 0
  488. && (three_state_map[i - 1][j] == GI_FULLY_PRESENT
  489. || three_state_map[i - 1][j] == GI_ABSENCE)) {
  490. //use same as LHS
  491. three_state_map[i][j] = three_state_map[i - 1][j];
  492. } else {
  493. int count = 0;
  494. for(int k=startPolygon;k<=endPolygon;k++) {
  495. count += isWithin(j * divx + divx / 2 + longitude1, i * divy + divy / 2 + latitude1)?1:0;
  496. }
  497. if(count % 2 == 1) {
  498. //if the previous was partially present, test
  499. three_state_map[i][j] = GI_FULLY_PRESENT;
  500. } //else absent
  501. } //else absent
  502. } else {
  503. //if the previous was fully present, repeat
  504. //if the previous was absent, repeat
  505. three_state_map[i][j] = three_state_map[i][j - 1];
  506. }
  507. //apply to cells;
  508. if (cellsReturned && three_state_map[i][j] != GI_UNDEFINED) { //undefined == absence
  509. data[p][0] = j;
  510. data[p][1] = i;
  511. p++;
  512. }
  513. }
  514. }
  515. if (data != null) {
  516. data = java.util.Arrays.copyOf(data, p);
  517. }
  518. return data;
  519. }
  520. public int[][] fillAccMask_EPSG900913(int startPolygon, int endPolygon, double longitude1, double latitude1, double longitude2, double latitude2, int width, int height, byte[][] three_state_map, boolean noCellsReturned) {
  521. double divx = (longitude2 - longitude1) / width;
  522. double divy = (latitude2 - latitude1) / height;
  523. int i, j;
  524. //do raster check
  525. int[][] data = null;
  526. boolean cellsReturned = !noCellsReturned;
  527. if (cellsReturned) {
  528. data = new int[width * height][2];
  529. }
  530. int p = 0;
  531. for (j = 0; j < three_state_map[0].length; j++) {
  532. for (i = 0; i < three_state_map.length; i++) {
  533. if (three_state_map[i][j] == GI_PARTIALLY_PRESENT) {
  534. //if it is partially present, do nothing
  535. } else if ((j == 0 || three_state_map[i][j - 1] == GI_PARTIALLY_PRESENT)) {
  536. if (i > 0
  537. && (three_state_map[i - 1][j] == GI_FULLY_PRESENT
  538. || three_state_map[i - 1][j] == GI_ABSENCE)) {
  539. //use same as LHS
  540. three_state_map[i][j] = three_state_map[i - 1][j];
  541. } else {
  542. int count = 0;
  543. for(int k=startPolygon;k<=endPolygon;k++) {
  544. count += isWithin_EPSG900913(j * divx + divx / 2 + longitude1, i * divy + divy / 2 + latitude1)?1:0;
  545. }
  546. if(count % 2 == 1) {
  547. //if the previous was partially present, test
  548. three_state_map[i][j] = GI_FULLY_PRESENT;
  549. } //else absent
  550. } //else absent
  551. } else {
  552. //if the previous was fully present, repeat
  553. //if the previous was absent, repeat
  554. three_state_map[i][j] = three_state_map[i][j - 1];
  555. }
  556. //apply to cells;
  557. if (cellsReturned && three_state_map[i][j] != GI_UNDEFINED) { //undefined == absence
  558. data[p][0] = j;
  559. data[p][1] = i;
  560. p++;
  561. }
  562. }
  563. }
  564. if (data != null) {
  565. data = java.util.Arrays.copyOf(data, p);
  566. }
  567. return data;
  568. }
  569. }