PageRenderTime 38ms CodeModel.GetById 12ms app.highlight 21ms RepoModel.GetById 0ms app.codeStats 0ms

/alaspatial/src/main/java/org/ala/spatial/util/ComplexRegion.java

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