1 | # Script to read a lidar file and output a (numeric) histogram for use in filtering |
---|
2 | # Author: Ben Taylor |
---|
3 | # Date: 15th April 2008 |
---|
4 | # Copyright 2008 Plymouth Marine Laboratory - you may use this software as you wish, |
---|
5 | # but no warranty of any kind is given |
---|
6 | |
---|
7 | import re |
---|
8 | import math |
---|
9 | import sys |
---|
10 | |
---|
11 | HEIGHT_COL = 3 |
---|
12 | NORTH_COL = 2 |
---|
13 | EAST_COL = 1 |
---|
14 | |
---|
15 | # Function to print an (approximate) ascii histogram |
---|
16 | # Uses log of histogram values |
---|
17 | def print_histogram(disp_histogram, bounds, totalpoints): |
---|
18 | keylist = disp_histogram.keys() |
---|
19 | keylist.sort() |
---|
20 | print "Region boundaries:" |
---|
21 | print "n: " + str(bounds["n"]) |
---|
22 | print "s: " + str(bounds["s"]) |
---|
23 | print "e: " + str(bounds["e"]) |
---|
24 | print "w: " + str(bounds["w"]) |
---|
25 | print "-------------------" |
---|
26 | print "Log scale x5 histogram" |
---|
27 | print "Total points: " + str(totalpoints) |
---|
28 | for key in keylist: |
---|
29 | value = disp_histogram[key] |
---|
30 | logval = math.log10(value) |
---|
31 | dispval = int(math.floor(5*logval)) # Scale value by 5 |
---|
32 | print str(key) + ": " + ("#"*dispval) + " (" + str(value) + ")" |
---|
33 | # end for |
---|
34 | # end function |
---|
35 | |
---|
36 | # Function to work out recommended upper/lower bounds for the calculated histogram |
---|
37 | def print_cutoffs(disp_histogram, totalpoints): |
---|
38 | |
---|
39 | # -1 denotes not yet set (negative value are legal but will be multiples of 10 if valid) |
---|
40 | lbound = -1 |
---|
41 | ubound = -1 |
---|
42 | |
---|
43 | # Sort list of keys |
---|
44 | keylist = disp_histogram.keys() |
---|
45 | keylist.sort() |
---|
46 | |
---|
47 | # Run through list. First bracket containing more than 1/10000th of all points is minimum, less 10m |
---|
48 | for key in keylist: |
---|
49 | value = disp_histogram[key] |
---|
50 | |
---|
51 | if (lbound == -1): |
---|
52 | if (value > 0.0001 * totalpoints): |
---|
53 | lbound = key - 10 |
---|
54 | break |
---|
55 | # end if |
---|
56 | # end if |
---|
57 | # end for |
---|
58 | |
---|
59 | # Run through list again backwards. First bracket containing more than 1/10000th of all points is maximum, |
---|
60 | # plus 10m (20 = width of bracket + 10m) |
---|
61 | keylist.reverse() |
---|
62 | for key in keylist: |
---|
63 | value = disp_histogram[key] |
---|
64 | |
---|
65 | if (ubound == -1): |
---|
66 | if (value > 0.0001 * totalpoints): |
---|
67 | ubound = key + 20 |
---|
68 | break |
---|
69 | # end if |
---|
70 | # end if |
---|
71 | # end for |
---|
72 | |
---|
73 | print "Recommended bounds:" |
---|
74 | print "Upper: " + str(ubound) |
---|
75 | print "Lower: " + str(lbound) |
---|
76 | print "-------------------" |
---|
77 | # end function |
---|
78 | |
---|
79 | if (len(sys.argv) != 2): |
---|
80 | print "Usage: python lidar_histogram.py input_filename" |
---|
81 | sys.exit(1) |
---|
82 | # end if |
---|
83 | |
---|
84 | filepath = sys.argv[1] |
---|
85 | |
---|
86 | try: |
---|
87 | datafile = open(filepath) |
---|
88 | except: |
---|
89 | print "Could not open input file: " + filepath |
---|
90 | sys.exit(1) |
---|
91 | # end try |
---|
92 | |
---|
93 | histogram = dict() # Empty dictionary to hold output |
---|
94 | boundaries = dict() # Empty dictionary to hold map boundaries |
---|
95 | |
---|
96 | eof = False |
---|
97 | |
---|
98 | pointcounter = 0 |
---|
99 | |
---|
100 | # For each file line, read it, split on whitespace, read the entries for height, northing and easting, |
---|
101 | # increment the relevant 10m height bracket in the histogram and update the map boundaries if appropriate |
---|
102 | while (not(eof)): |
---|
103 | dataline = datafile.readline() |
---|
104 | if (dataline == ""): |
---|
105 | eof = True |
---|
106 | break |
---|
107 | # end if |
---|
108 | |
---|
109 | pointcounter = pointcounter + 1 |
---|
110 | |
---|
111 | linesplit = re.split("\s+", dataline) |
---|
112 | lineheight = float(linesplit[HEIGHT_COL]) |
---|
113 | linenorth = float(linesplit[NORTH_COL]) |
---|
114 | lineeast = float(linesplit[EAST_COL]) |
---|
115 | |
---|
116 | heightindex = int(math.floor(lineheight) - (math.floor(lineheight) % 10)) # Get height bracket (round down to nearest 10m) |
---|
117 | |
---|
118 | if (heightindex in histogram): |
---|
119 | histogram[heightindex] = histogram[heightindex] + 1 |
---|
120 | else: |
---|
121 | histogram[heightindex] = 1 |
---|
122 | # end if |
---|
123 | |
---|
124 | if not ("n" in boundaries): |
---|
125 | boundaries["n"] = int(math.ceil(linenorth)) |
---|
126 | # end if |
---|
127 | |
---|
128 | if not ("s" in boundaries): |
---|
129 | boundaries["s"] = int(math.floor(linenorth)) |
---|
130 | # end if |
---|
131 | |
---|
132 | if not ("e" in boundaries): |
---|
133 | boundaries["e"] = int(math.ceil(lineeast)) |
---|
134 | # end if |
---|
135 | |
---|
136 | if not ("w" in boundaries): |
---|
137 | boundaries["w"] = int(math.floor(lineeast)) |
---|
138 | # end if |
---|
139 | |
---|
140 | if (linenorth > boundaries["n"]): |
---|
141 | boundaries["n"] = int(math.ceil(linenorth)) |
---|
142 | elif (linenorth < boundaries["s"]): |
---|
143 | boundaries["s"] = int(math.floor(linenorth)) |
---|
144 | # end if |
---|
145 | |
---|
146 | if (lineeast > boundaries["e"]): |
---|
147 | boundaries["e"] = int(math.ceil(lineeast)) |
---|
148 | elif (lineeast < boundaries["w"]): |
---|
149 | boundaries["w"] = int(math.floor(lineeast)) |
---|
150 | # end if |
---|
151 | |
---|
152 | # end while |
---|
153 | |
---|
154 | datafile.close() |
---|
155 | |
---|
156 | print_cutoffs(histogram, pointcounter) |
---|
157 | print_histogram(histogram, boundaries, pointcounter) |
---|