Data hosted with ♥ by Pastebin.com - Download Raw - See Original
  1. '''
  2. Created on 13/06/2013
  3.  
  4. Perform tile by tile true ortho using
  5. vector projection from each raster cell
  6. onto 3D-model
  7.  
  8. @author: Tisham
  9. '''
  10.  
  11. import argparse
  12. import glob
  13. import os
  14. from index_tiles import parse_tile_def
  15. from multiprocessing import Pool,cpu_count
  16. from functools import partial
  17. from cgkit.all import load,getScene,TriMeshGeom,vec3
  18. import numpy as np
  19. #from pylab import imshow,show
  20. import osgeo.gdal as gdal
  21. from rtree import index
  22. import cv2 as cv
  23.  
  24. def generate_tri(obj_base,tri_mesh_geom):
  25.     p = index.Property()
  26.     p.dimension = 3
  27.     p.dat_extension = 'data'
  28.     p.idx_extension = 'idx'
  29.    
  30.     idx3d = index.Index(obj_base,properties=p)
  31.     verts = tri_mesh_geom.verts
  32.     for i,face in enumerate(tri_mesh_geom.faces):
  33.         vert1 = verts[face[0]]
  34.         vert2 = verts[face[1]]
  35.         vert3 = verts[face[2]]
  36.         stack = np.vstack((vert1,vert2,vert3))
  37.         mintri = np.min(stack,axis=0)
  38.         maxtri = np.max(stack,axis=0)
  39.         idx3d.insert(i, np.hstack((mintri,maxtri)))
  40.     return idx3d
  41.        
  42.  
  43.  
  44. def copy_valid_tri(tri_geom, r_list, count):
  45.     """
  46.    Copy only valid triangles using pre-computed kd-tree
  47.    """
  48.     test_tri = TriMeshGeom()
  49.     test_tri.faces.resize(count)
  50.     test_tri.verts.resize(count * 3)
  51.     for k, face in enumerate(test_tri.faces):
  52.         thr_k = 3*k
  53.         test_tri.faces[k] = [thr_k, thr_k + 1,thr_k + 2]
  54.         test_tri.verts[thr_k] = tri_geom.verts[tri_geom.faces[r_list[k]][0]]
  55.         test_tri.verts[thr_k + 1] = tri_geom.verts[tri_geom.faces[r_list[k]][1]]
  56.         test_tri.verts[thr_k + 2] = tri_geom.verts[tri_geom.faces[r_list[k]][2]]
  57.    
  58.     return test_tri
  59.  
  60.  
  61. def clean_index(obj_base):
  62.     if os.path.exists(obj_base + ".data"):
  63.         os.remove(obj_base + ".data")
  64.         os.remove(obj_base + ".idx")
  65.  
  66.  
  67. def extract_hit_vt(tri_vt,u,v, oa_face):
  68.    
  69.     vt1 = np.array(vec3(tri_vt[oa_face * 3]))
  70.     vt2 = np.array(vec3(tri_vt[oa_face * 3 + 1]))
  71.     vt3 = np.array(vec3(tri_vt[oa_face * 3 + 2]))
  72.    
  73.     hit_vt = (1.0-u-v)*vt1+u*vt2+v*vt3
  74.    
  75.     return hit_vt
  76.  
  77. def process_dem_ortho(obj_file,tdef_folder,spacing,xoff,yoff):
  78.     """
  79.    Generate maximum height map from given OBJ tile
  80.    """
  81.     #Read tile def
  82.     obj_base = os.path.basename(obj_file).rstrip(".obj")
  83.    
  84.     #if not redo_tiles.__contains__(obj_base): return
  85.    
  86.     tile_def = os.path.join(tdef_folder,obj_base+"\\tile_overlapping.xml")
  87.     tile_bounds = parse_tile_def(tile_def)
  88.    
  89.     #Parse OBJ file with CGKit
  90.     scene =  getScene()
  91.     scene.clear()
  92.     load(obj_file)
  93.    
  94.     try:
  95.         print scene.boundingBox()
  96.         tri_mesh = scene.worldObject("Mesh")
  97.        
  98.         tri_geom = tri_mesh.geom
  99.         tri_vt = tri_mesh.geom.slot("st")
  100.        
  101.         textures = []
  102.         tex_shapes = []
  103.         num_tex = tri_mesh.getNumMaterials()
  104.        
  105.         for i in range(num_tex):
  106.             tex_file = tri_mesh.getMaterial(i).map_Kd.filename
  107.             tex_file = os.path.join(os.path.dirname(obj_file),tex_file)
  108.             texture = cv.imread(tex_file)
  109.             textures.append(cv.cvtColor(texture, cv.cv.CV_BGR2RGB))
  110.             tex_shapes.append(texture.shape)
  111.        
  112.         #Remove previous index
  113.         clean_index(obj_base)
  114.        
  115.         indexed_poly = generate_tri(obj_base,tri_mesh.geom)
  116.        
  117.         xstart = tile_bounds[0]
  118.         xend = tile_bounds[3]
  119.         ystart = tile_bounds[1]
  120.         yend = tile_bounds[4]
  121.        
  122.         xsize = int((xend - xstart)/spacing)
  123.         ysize = int((yend - ystart)/spacing)
  124.         half_s = spacing/2.0
  125.        
  126.         #Assign raster dem using tile_def + spacing
  127.         dem_buffer = np.zeros((ysize,xsize))
  128.        
  129.         #Assign raster ortho using tile_def + spacing
  130.         tex_buffer = np.zeros((ysize,xsize,3),dtype=np.ubyte)
  131.        
  132.         #Go through every XY in raster and extract
  133.         #intersect z
  134.         for i in range(ysize):
  135.             for j in range(xsize):
  136.                 x = xstart + spacing*j
  137.                 y = yend - spacing*i
  138.                
  139.                 try:
  140.                     r_list = list(indexed_poly.intersection((x-half_s,y-half_s,0,x+half_s,y+half_s,1000)))
  141.                     count = len(r_list)
  142.                     if(count>0):
  143.                         test_tri = copy_valid_tri(tri_geom, r_list, count)
  144.                        
  145.                         #hit_h, t_h, face_h, u_h, v_h = TriMeshGeom.intersectRay(test_tri, vec3(x,y,1000), vec3(0,0,-1))
  146.                         hit_l, t_l, face_l, u_l, v_l = TriMeshGeom.intersectRay(test_tri, vec3(x,y,0), vec3(0,0,1))
  147.                         if hit_l:
  148.                             #z = 1000-t_h
  149.                             z = t_l
  150.                             dem_buffer[i,j] = z
  151.                             #Use UV and extract intersect RGB
  152.                             oa_face = r_list[face_l]
  153.                             hit_vt = extract_hit_vt(tri_vt, u_l,v_l, oa_face)
  154.                            
  155.                             color_sub = (0,0,0)
  156.                             if num_tex == 1:
  157.                                 rows,cols,b = tex_shapes[0]
  158.                                 center = (hit_vt[0]*(cols-1),(1-hit_vt[1])*(rows-1))
  159.                                 color_sub = cv.getRectSubPix(textures[0], (1,1), center)
  160.                             else:
  161.                                 matid = tri_geom.slot("matid")[oa_face]
  162.                                 rows,cols,b = tex_shapes[matid]
  163.                                 center = (hit_vt[0]*(cols-1),(1-hit_vt[1])*(rows-1))
  164.                                 color_sub = cv.getRectSubPix(textures[matid], (1,1), center)
  165.                             #color = texture[(1-hit_vt[1])*(rows-1),hit_vt[0]*(cols-1),:]
  166.                             tex_buffer[i,j,:] = color_sub
  167.                 except Exception as e:
  168.                     print "Failed",e
  169.        
  170.         indexed_poly = None
  171.        
  172.         #All intersections done
  173.         #Remove previous index
  174.         clean_index(obj_base)
  175.        
  176.         #Save DEM file
  177.         """
  178.        imshow(dem_buffer)
  179.        show()
  180.        """
  181.        
  182.         driver = gdal.GetDriverByName("GTiff")
  183.         out_ds = driver.Create(obj_file.replace(".obj","_dem.tif"),xsize,ysize,1,gdal.GDT_Float32,options=["COMPRESS=LZW","TFW=YES"])
  184.         out_ds.GetRasterBand(1).WriteArray(dem_buffer)
  185.         out_ds.SetGeoTransform([xstart+xoff,spacing,0,yend+yoff,0,-spacing])
  186.         out_ds.FlushCache()
  187.         out_ds = None
  188.        
  189.         #Save Raster Image file
  190.         """
  191.        imshow(tex_buffer)
  192.        show()
  193.        """
  194.         driver = gdal.GetDriverByName("GTiff")
  195.         out_ds = driver.Create(obj_file.replace(".obj","_ortho.tif"),xsize,ysize,3,gdal.GDT_Byte,options=["COMPRESS=LZW","TFW=YES"] )
  196.         out_ds.GetRasterBand(1).WriteArray(tex_buffer[:,:,0])
  197.         out_ds.GetRasterBand(2).WriteArray(tex_buffer[:,:,1])
  198.         out_ds.GetRasterBand(3).WriteArray(tex_buffer[:,:,2])
  199.         out_ds.SetGeoTransform([xstart+xoff,spacing,0,yend+yoff,0,-spacing])
  200.         out_ds.FlushCache()
  201.         out_ds = None
  202.        
  203.     except Exception as e:
  204.         print "Mesh not loaded",obj_file
  205.         print "Exception", e
  206.  
  207.  
  208. if __name__ == "__main__":
  209.     parser = argparse.ArgumentParser(description="Parse OBJ, ray trace vertically and output tif tiles")
  210.     parser.add_argument("inobj",type=str,
  211.                         help="Input OBJ folder")
  212.     parser.add_argument("intile",type=str,
  213.                         help="Input Tile definition folder")
  214.     parser.add_argument("spacing",type=float,
  215.                         help="output sample spacing")
  216.     parser.add_argument("xoffset",type=float,
  217.                         help="XOffset to real co-ordinates")
  218.     parser.add_argument("yoffset",type=float,
  219.                         help="YOffset to real co-ordinates")
  220.    
  221.     args = parser.parse_args()
  222.    
  223.     obj_folder = args.inobj
  224.     tdef_folder = args.intile
  225.     spacing = args.spacing
  226.     xoffset = args.xoffset
  227.     yoffset = args.yoffset
  228.    
  229.     partial_ppd = partial(process_dem_ortho, tdef_folder=tdef_folder, spacing=spacing, xoff=xoffset,yoff=yoffset)
  230.    
  231.     obj_files = glob.glob(os.path.join(obj_folder,"Tile*\\*[+-]???.obj"))
  232.    
  233.     p = Pool(cpu_count()-1)
  234.     p.map(partial_ppd, obj_files)
  235.    
  236.     #for obj_file in obj_files:
  237.     #    partial_ppd(obj_file)