'''
Created on 13/06/2013
Perform tile by tile true ortho using
vector projection from each raster cell
onto 3D-model
@author: Tisham
'''
import argparse
import glob
import os
from index_tiles import parse_tile_def
from multiprocessing import Pool,cpu_count
from functools import partial
from cgkit.all import load,getScene,TriMeshGeom,vec3
import numpy as np
#from pylab import imshow,show
import osgeo.gdal as gdal
from rtree import index
import cv2 as cv
def generate_tri(obj_base,tri_mesh_geom):
p = index.Property()
p.dimension = 3
p.dat_extension = 'data'
p.idx_extension = 'idx'
idx3d = index.Index(obj_base,properties=p)
verts = tri_mesh_geom.verts
for i,face in enumerate(tri_mesh_geom.faces):
vert1 = verts[face[0]]
vert2 = verts[face[1]]
vert3 = verts[face[2]]
stack = np.vstack((vert1,vert2,vert3))
mintri = np.min(stack,axis=0)
maxtri = np.max(stack,axis=0)
idx3d.insert(i, np.hstack((mintri,maxtri)))
return idx3d
def copy_valid_tri(tri_geom, r_list, count):
"""
Copy only valid triangles using pre-computed kd-tree
"""
test_tri = TriMeshGeom()
test_tri.faces.resize(count)
test_tri.verts.resize(count * 3)
for k, face in enumerate(test_tri.faces):
thr_k = 3*k
test_tri.faces[k] = [thr_k, thr_k + 1,thr_k + 2]
test_tri.verts[thr_k] = tri_geom.verts[tri_geom.faces[r_list[k]][0]]
test_tri.verts[thr_k + 1] = tri_geom.verts[tri_geom.faces[r_list[k]][1]]
test_tri.verts[thr_k + 2] = tri_geom.verts[tri_geom.faces[r_list[k]][2]]
return test_tri
def clean_index(obj_base):
if os.path.exists(obj_base + ".data"):
os.remove(obj_base + ".data")
os.remove(obj_base + ".idx")
def extract_hit_vt(tri_vt,u,v, oa_face):
vt1 = np.array(vec3(tri_vt[oa_face * 3]))
vt2 = np.array(vec3(tri_vt[oa_face * 3 + 1]))
vt3 = np.array(vec3(tri_vt[oa_face * 3 + 2]))
hit_vt = (1.0-u-v)*vt1+u*vt2+v*vt3
return hit_vt
def process_dem_ortho(obj_file,tdef_folder,spacing,xoff,yoff):
"""
Generate maximum height map from given OBJ tile
"""
#Read tile def
obj_base = os.path.basename(obj_file).rstrip(".obj")
#if not redo_tiles.__contains__(obj_base): return
tile_def = os.path.join(tdef_folder,obj_base+"\\tile_overlapping.xml")
tile_bounds = parse_tile_def(tile_def)
#Parse OBJ file with CGKit
scene = getScene()
scene.clear()
load(obj_file)
try:
print scene.boundingBox()
tri_mesh = scene.worldObject("Mesh")
tri_geom = tri_mesh.geom
tri_vt = tri_mesh.geom.slot("st")
textures = []
tex_shapes = []
num_tex = tri_mesh.getNumMaterials()
for i in range(num_tex):
tex_file = tri_mesh.getMaterial(i).map_Kd.filename
tex_file = os.path.join(os.path.dirname(obj_file),tex_file)
texture = cv.imread(tex_file)
textures.append(cv.cvtColor(texture, cv.cv.CV_BGR2RGB))
tex_shapes.append(texture.shape)
#Remove previous index
clean_index(obj_base)
indexed_poly = generate_tri(obj_base,tri_mesh.geom)
xstart = tile_bounds[0]
xend = tile_bounds[3]
ystart = tile_bounds[1]
yend = tile_bounds[4]
xsize = int((xend - xstart)/spacing)
ysize = int((yend - ystart)/spacing)
half_s = spacing/2.0
#Assign raster dem using tile_def + spacing
dem_buffer = np.zeros((ysize,xsize))
#Assign raster ortho using tile_def + spacing
tex_buffer = np.zeros((ysize,xsize,3),dtype=np.ubyte)
#Go through every XY in raster and extract
#intersect z
for i in range(ysize):
for j in range(xsize):
x = xstart + spacing*j
y = yend - spacing*i
try:
r_list = list(indexed_poly.intersection((x-half_s,y-half_s,0,x+half_s,y+half_s,1000)))
count = len(r_list)
if(count>0):
test_tri = copy_valid_tri(tri_geom, r_list, count)
#hit_h, t_h, face_h, u_h, v_h = TriMeshGeom.intersectRay(test_tri, vec3(x,y,1000), vec3(0,0,-1))
hit_l, t_l, face_l, u_l, v_l = TriMeshGeom.intersectRay(test_tri, vec3(x,y,0), vec3(0,0,1))
if hit_l:
#z = 1000-t_h
z = t_l
dem_buffer[i,j] = z
#Use UV and extract intersect RGB
oa_face = r_list[face_l]
hit_vt = extract_hit_vt(tri_vt, u_l,v_l, oa_face)
color_sub = (0,0,0)
if num_tex == 1:
rows,cols,b = tex_shapes[0]
center = (hit_vt[0]*(cols-1),(1-hit_vt[1])*(rows-1))
color_sub = cv.getRectSubPix(textures[0], (1,1), center)
else:
matid = tri_geom.slot("matid")[oa_face]
rows,cols,b = tex_shapes[matid]
center = (hit_vt[0]*(cols-1),(1-hit_vt[1])*(rows-1))
color_sub = cv.getRectSubPix(textures[matid], (1,1), center)
#color = texture[(1-hit_vt[1])*(rows-1),hit_vt[0]*(cols-1),:]
tex_buffer[i,j,:] = color_sub
except Exception as e:
print "Failed",e
indexed_poly = None
#All intersections done
#Remove previous index
clean_index(obj_base)
#Save DEM file
"""
imshow(dem_buffer)
show()
"""
driver = gdal.GetDriverByName("GTiff")
out_ds = driver.Create(obj_file.replace(".obj","_dem.tif"),xsize,ysize,1,gdal.GDT_Float32,options=["COMPRESS=LZW","TFW=YES"])
out_ds.GetRasterBand(1).WriteArray(dem_buffer)
out_ds.SetGeoTransform([xstart+xoff,spacing,0,yend+yoff,0,-spacing])
out_ds.FlushCache()
out_ds = None
#Save Raster Image file
"""
imshow(tex_buffer)
show()
"""
driver = gdal.GetDriverByName("GTiff")
out_ds = driver.Create(obj_file.replace(".obj","_ortho.tif"),xsize,ysize,3,gdal.GDT_Byte,options=["COMPRESS=LZW","TFW=YES"] )
out_ds.GetRasterBand(1).WriteArray(tex_buffer[:,:,0])
out_ds.GetRasterBand(2).WriteArray(tex_buffer[:,:,1])
out_ds.GetRasterBand(3).WriteArray(tex_buffer[:,:,2])
out_ds.SetGeoTransform([xstart+xoff,spacing,0,yend+yoff,0,-spacing])
out_ds.FlushCache()
out_ds = None
except Exception as e:
print "Mesh not loaded",obj_file
print "Exception", e
if __name__ == "__main__":
parser = argparse.ArgumentParser(description="Parse OBJ, ray trace vertically and output tif tiles")
parser.add_argument("inobj",type=str,
help="Input OBJ folder")
parser.add_argument("intile",type=str,
help="Input Tile definition folder")
parser.add_argument("spacing",type=float,
help="output sample spacing")
parser.add_argument("xoffset",type=float,
help="XOffset to real co-ordinates")
parser.add_argument("yoffset",type=float,
help="YOffset to real co-ordinates")
args = parser.parse_args()
obj_folder = args.inobj
tdef_folder = args.intile
spacing = args.spacing
xoffset = args.xoffset
yoffset = args.yoffset
partial_ppd = partial(process_dem_ortho, tdef_folder=tdef_folder, spacing=spacing, xoff=xoffset,yoff=yoffset)
obj_files = glob.glob(os.path.join(obj_folder,"Tile*\\*[+-]???.obj"))
p = Pool(cpu_count()-1)
p.map(partial_ppd, obj_files)
#for obj_file in obj_files:
# partial_ppd(obj_file)