SHOW:
|
|
- or go back to the newest paste.
| 1 | import time | |
| 2 | - | from functools import partial |
| 2 | + | |
| 3 | - | from itertools import chain, imap |
| 3 | + | |
| 4 | ||
| 5 | ||
| 6 | def iter_points(layer): | |
| 7 | for feature in layer: | |
| 8 | exterior_ring = feature.GetGeometryRef().GetGeometryRef(0) | |
| 9 | for i in xrange(exterior_ring.GetPointCount()): | |
| 10 | yield exterior_ring.GetPoint(i) | |
| 11 | ||
| 12 | ||
| 13 | def clean_points(points): | |
| 14 | # | |
| 15 | # Delete the last point of each triangle (we only need three point for our | |
| 16 | # calculations) | |
| 17 | # | |
| 18 | # TODO: The code below has the effect of the original code but it doesn't | |
| 19 | # match the description above and it repeats points which will get | |
| 20 | # filtered out in a later processing step anyway. I guess the | |
| 21 | # `drop_every_nth()` function below is actually wanted here instead. | |
| 22 | # | |
| 23 | points = iter(points) | |
| 24 | try: | |
| 25 | while True: | |
| 26 | yield next(points) | |
| 27 | next(points) | |
| 28 | point = next(points) | |
| 29 | yield point | |
| 30 | yield point | |
| 31 | next(points) | |
| 32 | except StopIteration: | |
| 33 | pass # Ignored intentionally. | |
| 34 | ||
| 35 | ||
| 36 | def drop_every_nth(nth, iterable): | |
| 37 | return (item for i, item in enumerate(iterable, 1) if i % nth != 0) | |
| 38 | - | def drop_every_nth(iterable, nth): |
| 38 | + | |
| 39 | - | return (item for i, item in enumerate(iterable) if i % nth != 0) |
| 39 | + | |
| 40 | def distance(point_a, point_b): | |
| 41 | # | |
| 42 | - | def calculate_mid_value(a, b): |
| 42 | + | # TODO: To be the usual euklidean distance formula it lacks a `math.sqrt()` |
| 43 | - | return min(a, b) + abs(b - a) / 2.0 |
| 43 | + | # call! |
| 44 | # | |
| 45 | return sum((a - b)**2 for a, b in zip(point_a, point_b)) | |
| 46 | - | def partition_points(points, bounding_box): |
| 46 | + | |
| 47 | - | bbox_xmin, bbox_xmax, bbox_ymin, bbox_ymax = bounding_box |
| 47 | + | |
| 48 | - | vertical_middle = calculate_mid_value(bbox_xmin, bbox_xmax) |
| 48 | + | |
| 49 | - | horizontal_middle = calculate_mid_value(bbox_ymin, bbox_ymax) |
| 49 | + | |
| 50 | - | result = [[list(), list()], [list(), list()]] |
| 50 | + | # TODO: The reversing is only neccessary if the points have to be looked |
| 51 | # at in the order they are given to this function. If going through | |
| 52 | - | x, y, _ = point |
| 52 | + | # them in reverse would be okay we could spare this step. |
| 53 | # | |
| 54 | - | # TODO: Points *on* the lines won't get filtered out here. If that's |
| 54 | + | points = list(reversed(points)) |
| 55 | - | # a bug we need an extra test here. |
| 55 | + | while points: |
| 56 | point = points.pop() | |
| 57 | - | result[y < horizontal_middle][x < vertical_middle].append(point) |
| 57 | + | point_count = len(points) |
| 58 | points = [p for p in points if distance(point, p) >= minimal_distance] | |
| 59 | - | for row in result: |
| 59 | + | if len(points) == point_count: |
| 60 | - | for cell in row: |
| 60 | + | |
| 61 | - | yield cell |
| 61 | + | |
| 62 | ||
| 63 | def create_geometry(points): | |
| 64 | geometry = ogr.Geometry(ogr.wkbMultiPoint) | |
| 65 | wkb_point = ogr.Geometry(ogr.wkbPoint) | |
| 66 | for point in points: | |
| 67 | wkb_point.AddPoint(*point) | |
| 68 | geometry.AddGeometry(wkb_point) | |
| 69 | - | seen = set() |
| 69 | + | |
| 70 | ||
| 71 | - | if ( |
| 71 | + | |
| 72 | - | point not in seen |
| 72 | + | |
| 73 | - | and all(distance(point, p) >= minimal_distance for p in points) |
| 73 | + | |
| 74 | - | ): |
| 74 | + | |
| 75 | - | seen.add(point) |
| 75 | + | |
| 76 | identifier=__name__, | |
| 77 | title='Pollmueller Tin-Reducer', | |
| 78 | version = '0.1', | |
| 79 | storeSupported='true', | |
| 80 | statusSupported='true', | |
| 81 | abstract=( | |
| 82 | 'reduces the points of a given triangulated irregular network' | |
| 83 | ' by a variable factor' | |
| 84 | ), | |
| 85 | grassLocation=False | |
| 86 | ) | |
| 87 | self.data_output = self.addComplexOutput( | |
| 88 | identifier='output', | |
| 89 | title='Output vector data', | |
| 90 | formats=[{'mimeType': 'text/xml'}]
| |
| 91 | ) | |
| 92 | self.time_output = self.addLiteralOutput( | |
| 93 | identifier='output2', | |
| 94 | title='computation time (to compare different approaches)' | |
| 95 | ) | |
| 96 | self.data_filename_input = self.addComplexInput( | |
| 97 | identifier='data', | |
| 98 | title='Input vector data', | |
| 99 | formats=[{'mimeType': 'text/xml'}]
| |
| 100 | ) | |
| 101 | self.minimal_distance_input = self.addLiteralInput( | |
| 102 | identifier='input2', title='minimal distance between points' | |
| 103 | ) | |
| 104 | ||
| 105 | def execute(self): | |
| 106 | # | |
| 107 | # TODO: Different operating systems have different semantics and | |
| 108 | # resolutions for `time()` and `clock()`. Make sure `clock()` is | |
| 109 | # the correct one here. | |
| 110 | # | |
| 111 | start_timestamp = time.clock() | |
| 112 | ||
| 113 | shape_file = ogr.Open(self.data_filename_input.getValue()) | |
| 114 | if shape_file is None: | |
| 115 | raise Exception('failed to open shape file.')
| |
| 116 | # | |
| 117 | # Get the first layer. | |
| 118 | # | |
| 119 | layer = shape_file.GetLayer(0) | |
| 120 | layer.ResetReading() | |
| 121 | points = filter_close_points( | |
| 122 | self.minimal_distance_input.getValue(), | |
| 123 | clean_points(iter_points(layer)) | |
| 124 | ) | |
| 125 | gml = create_geometry(points).ExportToGML() | |
| 126 | ||
| 127 | elapsed_time = time.clock() - start_timestamp | |
| 128 | self.data_output.setValue(gml) | |
| 129 | self.time_output.setValue(elapsed_time) |