Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- '''
- 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)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement