Sanitized 3D Hough transform

pull/16/head
Cees Bassa 2019-04-22 16:11:46 +02:00
parent 029da7aacd
commit 783558e640
1 changed files with 30 additions and 18 deletions

View File

@ -5,8 +5,6 @@ import subprocess
from stvid.stio import fourframe
from stvid.stio import satid
from stvid.stio import observation
from tempfile import mkstemp
def generate_satellite_predictions(fname):
# Format command
@ -29,21 +27,24 @@ def find_hough3d_lines(fname, ntrkmin=20, dr=8):
# Compute selection mask
x, y, z, t, sig = ff.selection_mask(5.0, 40.0)
# Skip if not enough points
if len(t)<2:
return []
# Save points to temporary file
f, fpath = mkstemp()
with open(fpath, "w") as f:
with open("/tmp/hough.dat", "w") as f:
for i in range(len(t)):
f.write("%f,%f,%f\n" % (x[i], y[i], z[i]))
f.close()
# Run 3D Hough line-finding algorithm
command = "hough3dlines -dx %d -minvotes %d %s" % (dr, ntrkmin, fpath)
output = subprocess.check_output(command,
shell=True,
stderr=subprocess.STDOUT)
# Remove file
os.remove(fpath)
command = "hough3dlines -dx %d -minvotes %d %s" % (dr, ntrkmin, "/tmp/hough.dat")
try:
output = subprocess.check_output(command,
shell=True,
stderr=subprocess.STDOUT)
except:
return []
# Clean output (a bit cluncky)
cleaned_output = output.replace("npoints=", "")
@ -51,19 +52,30 @@ def find_hough3d_lines(fname, ntrkmin=20, dr=8):
cleaned_output = cleaned_output.replace("), b=(", " ")
cleaned_output = cleaned_output.replace(")", "")
cleaned_output = cleaned_output.replace(",", " ")
print(cleaned_output)
# Generate identifications
lines = []
s = cleaned_output.split()
for i in range(len(s)//7):
# Extract points
x0, y0, z0 = float(s[1+7*i]), float(s[2+7*i]), float(s[3+7*i])
dx, dy, dz = float(s[4+7*i]), float(s[5+7*i]), float(s[6+7*i])
xmin = x0-z0*dx/dz
xmax = x0+(ff.nz-z0)*dx/dz
ymin = y0-z0*dy/dz
ymax = y0+(ff.nz-z0)*dy/dz
line = "%.23s %8.3f %8.3f %8.3f %8.3f %8.5f %s unidentified sunlit" %\
(ff.nfd, xmin, ymin, xmax, ymax, ff.texp, 99999)
# Reconstruct start and end points
xmin = x0-z0*dx/(dz+1e-9)
xmax = x0+(ff.nz-z0)*dx/(dz+1e-9)
ymin = y0-z0*dy/(dz+1e-9)
ymax = y0+(ff.nz-z0)*dy/(dz+1e-9)
# Output line
line = "%.23s %8.3f %8.3f %8.3f %8.3f %8.5f %s unidentified sunlit\n" %\
(ff.nfd, xmin, ymin, xmax, ymax, ff.texp, 99999-i)
lines.append(line)
# Store identifications
fp = open(fname + ".id", "a")
for line in lines:
fp.write(line)
fp.close()
return [satid(line) for line in lines]