forked from hwiyoung/Orthophoto_Maps
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Orthophoto.py
93 lines (73 loc) · 3.81 KB
/
Orthophoto.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
import os
import numpy as np
import cv2
import time
from module.ExifData import *
from module.EoData import readEO, geographic2plane, Rot3D
from module.Boundary import boundary
from module.BackprojectionResample import projectedCoord, backProjection, resample, createGeoTiff
if __name__ == '__main__':
ground_height = 65 # unit: m
sensor_width = 6.3 # unit: mm
for root, dirs, files in os.walk('./Data'):
for file in files:
image_start_time = time.time()
start_time = time.time()
filename = os.path.splitext(file)[0]
extension = os.path.splitext(file)[1]
file_path = root + '/' + file
if extension == '.JPG' or extension == '.jpg':
print('Read the image - ' + file)
image = cv2.imread(file_path, -1)
# 1. Extract EXIF data from a image
focal_length, orientation = getExif(file_path) # unit: m, _
# 2. Restore the image based on orientation information
restored_image = restoreOrientation(image, orientation)
# 3. Convert pixel values into temperature
image_rows = restored_image.shape[0]
image_cols = restored_image.shape[1]
pixel_size = sensor_width / image_cols # unit: mm/px
pixel_size = pixel_size / 1000 # unit: m/px
end_time = time.time()
print("--- %s seconds ---" % (time.time() - start_time))
read_time = end_time - start_time
else:
print('Read EOP - ' + file)
print('Latitude | Longitude | Height | Omega | Phi | Kappa')
eo = readEO(file_path)
eo = geographic2plane(eo)
R = Rot3D(eo)
# 4. Extract a projected boundary of the image
bbox = boundary(restored_image, eo, R, ground_height, pixel_size, focal_length)
print("--- %s seconds ---" % (time.time() - start_time))
# 5. Compute GSD & Boundary size
# GSD
gsd = (pixel_size * (eo[2] - ground_height)) / focal_length # unit: m/px
# Boundary size
boundary_cols = int((bbox[1, 0] - bbox[0, 0]) / gsd)
boundary_rows = int((bbox[3, 0] - bbox[2, 0]) / gsd)
# 6. Compute coordinates of the projected boundary
print('projectedCoord')
start_time = time.time()
proj_coords = projectedCoord(bbox, boundary_rows, boundary_cols, gsd, eo, ground_height)
print("--- %s seconds ---" % (time.time() - start_time))
# Image size
image_size = np.reshape(restored_image.shape[0:2], (2, 1))
# 6. Back-projection into camera coordinate system
print('backProjection')
start_time = time.time()
backProj_coords = backProjection(proj_coords, R, focal_length, pixel_size, image_size)
print("--- %s seconds ---" % (time.time() - start_time))
# 7. Resample the pixels
print('resample')
start_time = time.time()
b, g, r, a = resample(backProj_coords, boundary_rows, boundary_cols, image)
print("--- %s seconds ---" % (time.time() - start_time))
# 8. Create GeoTiff
print('Save the image in GeoTiff')
start_time = time.time()
dst = './' + filename
createGeoTiff(b, g, r, a, bbox, gsd, boundary_rows, boundary_cols, dst)
print("--- %s seconds ---" % (time.time() - start_time))
print('*** Processing time per each image')
print("--- %s seconds ---" % (time.time() - image_start_time + read_time))