1"""
2Generate DTM in Blender
3
4Using the Displace and Decimate modifiers to generate a terrain 3D model from
5(geo)images (png, tif, geotif, etc)
6
7usage: blender -b -P blend_dtm.py -- dsm.tif image.jpg [terrain_resolution [strength]]
8
9The DSM must be Float32, Blender does not support Float64, If you get the error:
10    imb_loadtiff: Failed to read tiff image.
11    IMB_ibImageFromMemory: unknown fileformat
12Convert first using: gdal_translate -ot Float32 dsm.tif dsm32.tif
13
14You can safely ignore the "unknown field with tag" TIFFReadDirectory Warning.
15"""
16import os
17import sys
18import json
19import subprocess
20import bpy # Blender Python API
21
22TERRAIN_RESOLUTION=0.50 # in meter
23
24def new_texture():
25    bpy.ops.texture.new()
26    return bpy.data.textures[-1]
27
28def new_material():
29    bpy.ops.material.new()
30    return bpy.data.materials[-1]
31
32def open_image(filepath):
33    bpy.ops.image.open(filepath=filepath)
34    return bpy.data.images[-1]
35
36def new_image_texture(material, image, name='img'):
37    """ Create a new texture from an image (which we open)
38
39    we need a material here, to set the texture active
40    otherwise we can't change its type (to image in this case)
41    """
42    material.active_texture = new_texture()
43    material.active_texture.name = name
44    material.active_texture.type = 'IMAGE'
45    material.active_texture.image = image
46    return material.active_texture
47
48def new_grid(name='Grid', x_subdivisions=10, y_subdivisions=10, radius=1, scale=(1,1,1)):
49    bpy.ops.object.select_all(action='DESELECT')
50    # if blender version < 2.68 : radius -> size
51    if bpy.app.version < (2,68):
52        bpy.ops.mesh.primitive_grid_add(x_subdivisions=x_subdivisions, \
53        y_subdivisions=y_subdivisions, size=radius)
54    else:
55        bpy.ops.mesh.primitive_grid_add(x_subdivisions=x_subdivisions, \
56            y_subdivisions=y_subdivisions, radius=radius)
57    bpy.context.object.scale = scale
58    bpy.context.object.name = name
59    return bpy.context.object
60
61def select_only(obj):
62    bpy.ops.object.select_all(action='DESELECT')
63    obj.select = True
64    bpy.context.scene.objects.active = obj
65
66def subdivide(obj, number_cuts=1):
67    select_only(obj)
68    bpy.ops.object.mode_set(mode='EDIT')
69    bpy.ops.mesh.subdivide(number_cuts=number_cuts)
70    bpy.ops.object.mode_set(mode='OBJECT')
71    return obj
72
73def new_cube(name='Cube', scale=(1,1,1)):
74    bpy.ops.object.select_all(action='DESELECT')
75    bpy.ops.mesh.primitive_cube_add()
76    bpy.context.object.scale = scale
77    bpy.context.object.name = name
78    return bpy.context.object
79
80def set_viewport(viewport_shade='WIREFRAME', clip_end=1000):
81    """ Set the default view mode
82
83    :param viewport_shade: enum in ['BOUNDBOX', 'WIREFRAME', 'SOLID', 'TEXTURED'], default 'WIREFRAME'
84    """
85    for area in bpy.context.window.screen.areas:
86        if area.type == 'VIEW_3D':
87            for space in area.spaces:
88                if space.type == 'VIEW_3D':
89                    space.viewport_shade = viewport_shade
90                    space.clip_end = clip_end
91                    space.grid_scale = 10
92                    space.grid_lines = 50
93
94def save(filepath=None, check_existing=False, compress=True):
95    """ Save .blend file
96
97    :param filepath: File Path
98    :type  filepath: string, (optional, default: current file)
99    :param check_existing: Check and warn on overwriting existing files
100    :type  check_existing: boolean, (optional, default: False)
101    :param compress: Compress, Write compressed .blend file
102    :type  compress: boolean, (optional, default: True)
103    """
104    if not filepath:
105        filepath = bpy.data.filepath
106    bpy.ops.wm.save_mainfile(filepath=filepath, check_existing=check_existing,
107            compress=compress)
108
109def setup():
110    """ Setup the scene """
111    # Delete the default objects
112    bpy.ops.object.select_all(action='SELECT')
113    bpy.ops.object.delete()
114    # Add light (sun)
115    bpy.ops.object.lamp_add(type='SUN', location=(0, 0, 300))
116    # Set Game mode
117    bpy.context.scene.render.engine = 'BLENDER_GAME'
118    # make sure OpenGL shading language shaders (GLSL) is the
119    # material mode to use for rendering
120    bpy.context.scene.game_settings.material_mode = 'GLSL'
121    # Set the unit system to use for button display (in edit mode) to metric
122    bpy.context.scene.unit_settings.system = 'METRIC'
123    # Select the type of Framing to Extend,
124    # Show the entire viewport in the display window,
125    # viewing more horizontally or vertically.
126    bpy.context.scene.game_settings.frame_type = 'EXTEND'
127    # see the texture in realtime (in the 3D viewport)
128    set_viewport('TEXTURED', 10000)
129    # Set the color at the horizon to dark azure
130    bpy.context.scene.world.horizon_color = (0.05, 0.22, 0.4)
131    # Display framerate and profile information of the simulation
132    bpy.context.scene.game_settings.show_framerate_profile = True
133
134def usage():
135    sys.stderr.write(__doc__)
136
137def ext_exec(cmd, python=None):
138    if not python:
139        python = 'python%i.%i'%(sys.version_info.major, sys.version_info.minor)
140    return subprocess.getoutput('%s -c"%s"' % (python, cmd) )
141
142def fix_python_path(python=None):
143    pythonpath = ext_exec("import os,sys;print(os.pathsep.join(sys.path))")
144    sys.path.extend(pythonpath.split(os.pathsep))
145
146def get_gdalinfo(filepath):
147    """ Get the GeoTransform from the :param:filepath DEM using gdal
148
149    In a north up image, padfTransform[1] is the pixel width,
150    and padfTransform[5] is the pixel height. The upper left corner of the
151    upper left pixel is at position (padfTransform[0],padfTransform[3]).
152
153    The default transform is (0,1,0,0,0,1)
154    """
155    # Run externally since `gdal.Open` crashes Blender :-/
156    meta = ext_exec("import json,gdal,gdalconst;"
157        "g=gdal.Open('%s',gdalconst.GA_ReadOnly);"
158        "print(json.dumps({'transform':g.GetGeoTransform(),"
159            "'minmax':g.GetRasterBand(1).ComputeRasterMinMax(),"
160            "'meta':g.GetMetadata()}))"%filepath, 'python2')
161    return json.loads(meta) # {'transform':(0,1,0,0,0,1)}
162
163def main(argv=[]):
164    args = []
165    if '--' in argv:
166        args = argv[argv.index('--')+1:]
167
168    if len(args) < 2:
169        usage()
170        return 1
171
172    strength = 1
173    if len(args) > 2:
174        terrain_resolution = float(args[2])
175        if len(args) > 3:
176            strength = float(args[3])
177    else:
178        terrain_resolution = TERRAIN_RESOLUTION
179
180    image_dem = open_image(args[0])
181    image_img = open_image(args[1])
182    if image_dem.size[:] != image_img.size[:]:
183        print('[WARN] The size of the 2 image differ')
184
185    fix_python_path()
186    gdalinfo = get_gdalinfo(args[0])
187    geot = gdalinfo['transform']
188    meta = gdalinfo['meta']
189    xsize = image_dem.size[0] * abs(geot[1]) # in meters
190    ysize = image_dem.size[1] * abs(geot[5]) # in meters
191
192    translation = [0.0, 0.0, 0.0]
193    if 'CUSTOM_X_ORIGIN' in meta and 'CUSTOM_Y_ORIGIN' in meta:
194        custom_x_origin = float(meta['CUSTOM_X_ORIGIN'])
195        custom_y_origin = float(meta['CUSTOM_Y_ORIGIN'])
196        print("[gdal] got custom (%f, %f)"%(custom_x_origin, custom_y_origin))
197        center_x_utm = geot[0] + image_dem.size[0] * geot[1] / 2
198        center_y_utm = geot[3] + image_dem.size[1] * geot[5] / 2
199        translation = [ center_x_utm - custom_x_origin,
200                        center_y_utm - custom_y_origin, 0.0 ]
201    if gdalinfo['minmax']:
202        translation[2] = - round(gdalinfo['minmax'][0] + 1)
203    setup()
204
205    #########################################################################
206    # Add our ground object
207    ground = new_grid('Ground', xsize/terrain_resolution, \
208        ysize/terrain_resolution, 1, (xsize/2.0, ysize/2.0, 1))
209
210    #########################################################################
211    # Add material
212    ground.active_material = new_material()
213    ground.active_material.specular_intensity = 0
214    ground.active_material.name = 'Ground'
215    material = ground.active_material
216
217    #########################################################################
218    # Add texture
219    texture_img = new_image_texture(material, image_img, 'img')
220
221    material.active_texture_index += 1
222
223    if bpy.app.version > (2,63):
224        # since 2.64, default colorspace is sRGB which breaks decimate modifier
225        image_dem.colorspace_settings.name = 'Linear'
226    texture_dem = new_image_texture(material, image_dem, 'dem')
227    # do not show on the material (use only later for the terrain modifier)
228    material.use_textures[material.active_texture_index] = False
229
230    #bpy.ops.object.shade_smooth()
231    # displace from image_dem (gdal)
232    add_displace_modifier(ground, texture_dem, strength=strength, apply=True)
233    # unlink dem after apply (reduce size)
234    material.texture_slots.clear(1)
235    image_dem.user_clear()
236    bpy.data.images.remove(image_dem)
237    texture_dem.user_clear()
238    bpy.data.textures.remove(texture_dem)
239    #add_decimate_modifier(ground)
240    #add_water_cube((xsize+1, ysize+1, max(image_dem.pixels)/2))
241
242    # move to: custom origin - center (utm from gdalinfo)
243    ground.location = translation
244
245    bpy.ops.file.pack_all() # bpy.ops.image.pack(as_png=True)
246    save('/tmp/dtm.blend')
247    # bpy.ops.view3d.view_selected(use_all_regions=True) # XXX no poll()
248    print("\n----\nsaved in /tmp/dtm.blend\n----\n")
249    return 0
250
251def add_water_cube(scale=(1,1,1)):
252    """ Add cube 'water' to hide the bottom
253
254    (save memory, don't render many vertices)
255    """
256    cube = new_cube('Water', scale)
257    cube.active_material = new_material()
258    cube.active_material.name = 'Water'
259    cube.active_material.use_shadeless = True
260    cube.active_material.diffuse_color = (0.3, 0.4, 0.8)
261    return cube
262
263def add_displace_modifier(obj, texture, strength=1, direction='Z', apply=False):
264    """ Add displace modifier
265
266    http://wiki.blender.org/index.php/Doc:2.6/Manual/Modifiers/Deform/Displace
267    """
268    select_only(obj)
269    bpy.ops.object.modifier_add(type='DISPLACE')
270    modifier = bpy.context.object.modifiers[-1]
271    modifier.texture = texture
272    modifier.strength = strength
273    modifier.direction = direction
274    if apply:
275        bpy.ops.object.modifier_apply(modifier=modifier.name)
276    return modifier
277
278def add_decimate_modifier(obj, ratio=0.5, apply=False, show_viewport=False):
279    """ Add displace modifier
280
281    http://wiki.blender.org/index.php/Doc:2.6/Manual/Modifiers/Generate/Decimate
282    """
283    select_only(obj)
284    bpy.ops.object.modifier_add(type='DECIMATE')
285    modifier = bpy.context.object.modifiers[-1]
286    modifier.show_viewport = show_viewport # no preview : faster
287    modifier.ratio = ratio
288    if apply:
289        bpy.ops.object.modifier_apply(modifier=modifier.name)
290    return modifier
291
292def add_hillshade(demfile, material):
293    """ Add hillshade """
294    hillshade_png = '%s.hillshade.png'%demfile
295    if not os.path.isfile(hillshade_png):
296        subprocess.getstatusoutput("gdaldem hillshade %s %s -of png"%(demfile, hillshade_png))
297    image_hsh = open_image(hillshade_png)
298    material.active_texture_index += 1
299    texture_hsh = new_image_texture(material, image_hsh, 'hsh')
300    texture_hsh.use_normal_map = True
301    texture_slot_hsh = material.texture_slots[material.active_texture_index]
302    texture_slot_hsh.use_map_color_diffuse = False
303    texture_slot_hsh.use_map_normal = True
304    texture_slot_hsh.normal_factor = -0.5
305
306if __name__ == '__main__':
307    main(sys.argv)
308