-
Notifications
You must be signed in to change notification settings - Fork 47
/
get_unique_craters.py
204 lines (170 loc) · 6.85 KB
/
get_unique_craters.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
#!/usr/bin/env python
"""Unique Crater Distribution Functions
Functions for extracting craters from model target predictions and filtering
out duplicates.
"""
from __future__ import absolute_import, division, print_function
import numpy as np
import h5py
import sys
import utils.template_match_target as tmt
import utils.processing as proc
import utils.transform as trf
#from keras.models import load_model
#########################
def get_model_preds(CP):
"""Reads in or generates model predictions.
Parameters
----------
CP : dict
Containins directory locations for loading data and storing
predictions.
Returns
-------
craters : h5py
Model predictions.
"""
n_imgs, dtype = CP['n_imgs'], CP['datatype']
data = h5py.File(CP['dir_data'], 'r')
Data = {
dtype: [data['input_images'][:n_imgs].astype('float32'),
data['target_masks'][:n_imgs].astype('float32')]
}
data.close()
proc.preprocess(Data)
model = load_model(CP['dir_model'])
preds = model.predict(Data[dtype][0])
# save
h5f = h5py.File(CP['dir_preds'], 'w')
h5f.create_dataset(dtype, data=preds)
print("Successfully generated and saved model predictions.")
return preds
#########################
def add_unique_craters(craters, craters_unique, thresh_longlat2, thresh_rad):
"""Generates unique crater distribution by filtering out duplicates.
Parameters
----------
craters : array
Crater tuples from a single image in the form (long, lat, radius).
craters_unique : array
Master array of unique crater tuples in the form (long, lat, radius)
thresh_longlat2 : float.
Hyperparameter that controls the minimum squared longitude/latitude
difference between craters to be considered unique entries.
thresh_rad : float
Hyperparaeter that controls the minimum squared radius difference
between craters to be considered unique entries.
Returns
-------
craters_unique : array
Modified master array of unique crater tuples with new crater entries.
"""
k2d = 180. / (np.pi * 1737.4) # km to deg
Long, Lat, Rad = craters_unique.T
for j in range(len(craters)):
lo, la, r = craters[j].T
la_m = (la + Lat) / 2.
minr = np.minimum(r, Rad) # be liberal when filtering dupes
# duplicate filtering criteria
dL = (((Long - lo) / (minr * k2d / np.cos(np.pi * la_m / 180.)))**2
+ ((Lat - la) / (minr * k2d))**2)
dR = np.abs(Rad - r) / minr
index = (dR < thresh_rad) & (dL < thresh_longlat2)
if len(np.where(index == True)[0]) == 0:
craters_unique = np.vstack((craters_unique, craters[j]))
return craters_unique
#########################
def estimate_longlatdiamkm(dim, llbd, distcoeff, coords):
"""First-order estimation of long/lat, and radius (km) from
(Orthographic) x/y position and radius (pix).
For images transformed from ~6000 pixel crops of the 30,000 pixel
LROC-Kaguya DEM, this results in < ~0.4 degree latitude, <~0.2
longitude offsets (~2% and ~1% of the image, respectively) and ~2% error in
radius. Larger images thus may require an exact inverse transform,
depending on the accuracy demanded by the user.
Parameters
----------
dim : tuple or list
(width, height) of input images.
llbd : tuple or list
Long/lat limits (long_min, long_max, lat_min, lat_max) of image.
distcoeff : float
Ratio between the central heights of the transformed image and original
image.
coords : numpy.ndarray
Array of crater x coordinates, y coordinates, and pixel radii.
Returns
-------
craters_longlatdiamkm : numpy.ndarray
Array of crater longitude, latitude and radii in km.
"""
# Expand coords.
long_pix, lat_pix, radii_pix = coords.T
# Determine radius (km).
km_per_pix = 1. / trf.km2pix(dim[1], llbd[3] - llbd[2], dc=distcoeff)
radii_km = radii_pix * km_per_pix
# Determine long/lat.
deg_per_pix = km_per_pix * 180. / (np.pi * 1737.4)
long_central = 0.5 * (llbd[0] + llbd[1])
lat_central = 0.5 * (llbd[2] + llbd[3])
# Iterative method for determining latitude.
lat_deg_firstest = lat_central - deg_per_pix * (lat_pix - dim[1] / 2.)
latdiff = abs(lat_central - lat_deg_firstest)
# Protect against latdiff = 0 situation.
latdiff[latdiff < 1e-7] = 1e-7
lat_deg = lat_central - (deg_per_pix * (lat_pix - dim[1] / 2.) *
(np.pi * latdiff / 180.) /
np.sin(np.pi * latdiff / 180.))
# Determine longitude using determined latitude.
long_deg = long_central + (deg_per_pix * (long_pix - dim[0] / 2.) /
np.cos(np.pi * lat_deg / 180.))
# Return combined long/lat/radius array.
return np.column_stack((long_deg, lat_deg, radii_km))
def extract_unique_craters(CP, craters_unique):
"""Top level function that extracts craters from model predictions,
converts craters from pixel to real (degree, km) coordinates, and filters
out duplicate detections across images.
Parameters
----------
CP : dict
Crater Parameters needed to run the code.
craters_unique : array
Empty master array of unique crater tuples in the form
(long, lat, radius).
Returns
-------
craters_unique : array
Filled master array of unique crater tuples.
"""
# Load/generate model preds
try:
preds = h5py.File(CP['dir_preds'], 'r')[CP['datatype']]
print("Loaded model predictions successfully")
except:
print("Couldnt load model predictions, generating")
preds = get_model_preds(CP)
# need for long/lat bounds
P = h5py.File(CP['dir_data'], 'r')
llbd, pbd, distcoeff = ('longlat_bounds', 'pix_bounds',
'pix_distortion_coefficient')
#r_moon = 1737.4
dim = (float(CP['dim']), float(CP['dim']))
N_matches_tot = 0
for i in range(CP['n_imgs']):
id = proc.get_id(i)
coords = tmt.template_match_t(preds[i])
# convert, add to master dist
if len(coords) > 0:
new_craters_unique = estimate_longlatdiamkm(
dim, P[llbd][id], P[distcoeff][id][0], coords)
N_matches_tot += len(coords)
# Only add unique (non-duplicate) craters
if len(craters_unique) > 0:
craters_unique = add_unique_craters(new_craters_unique,
craters_unique,
CP['llt2'], CP['rt'])
else:
craters_unique = np.concatenate((craters_unique,
new_craters_unique))
np.save(CP['dir_result'], craters_unique)
return craters_unique