Skip to content

Render shadows #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
danlipsa wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Render shadows #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a22c6d3
Fix render to PNG
danlipsa Nov 8, 2018
bc0e791
Add script to create shadow mask.
danlipsa Nov 27, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 16 additions & 1 deletion danesfield/gdal_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
import re
import ogr
import osr

from vtk.util import numpy_support

def gdal_bounding_box(raster, outProj=None):
"""
Expand Down Expand Up @@ -148,3 +148,18 @@ def read_offset(fileName, offset):
if match:
for i in range(3):
offset[i] = float(match.group(1+i))

def vtk_to_numpy_order(aFlatVtk, dimensions):
'''
Convert a 2D array from VTK order to numpy order
'''
# VTK to numpy
aFlat = numpy_support.vtk_to_numpy(aFlatVtk)
# VTK X,Y corresponds to numpy cols,rows. VTK stores as
# in Fortran order.
aTranspose = numpy.reshape(aFlat, dimensions, "F")
# changes from cols, rows to rows,cols.
a = numpy.transpose(aTranspose)
# numpy rows increase as you go down, Y for VTK images increases as you go up
a = numpy.flip(a, 0)
return a
37 changes: 22 additions & 15 deletions danesfield/ortho.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,20 +31,22 @@ def circ_structure(n):

def orthorectify(args_source_image, args_dsm, args_destination_image,
args_occlusion_thresh=1.0, args_denoise_radius=2,
args_raytheon_rpc=None, args_dtm=None):
args_raytheon_rpc=None, args_dtm=None, args_convert_to_latlon=True):
"""
Orthorectify an image given the DSM

Args:
source_image: Source image file name
dsm: Digital surface model (DSM) image file name
destination_image: Orthorectified image file name
occlusion-thresh: Threshold on height difference for detecting
occlusion_thresh: Threshold on height difference for detecting
and masking occluded regions (in meters)
denoise-radius: Apply morphological operations with this radius
denoise_radius: Apply morphological operations with this radius
to the DSM reduce speckled noise
raytheon-rpc: Raytheon RPC file name. If not provided
raytheon_rpc: Raytheon RPC file name or model object. If not provided
the RPC is read from the source_image
dtm: Optional DTM parameter used to replace nodata areas in the
orthorectified image

Returns:
COMPLETE_DSM_INTERSECTION = 0
Expand All @@ -60,9 +62,12 @@ def orthorectify(args_source_image, args_dsm, args_destination_image,
sourceBand = sourceImage.GetRasterBand(1)

if (args_raytheon_rpc):
# read the RPC from raytheon file
print("Reading RPC from Raytheon file: {}".format(args_raytheon_rpc))
model = raytheon_rpc.read_raytheon_rpc_file(args_raytheon_rpc)
if (isinstance(args_raytheon_rpc, str)):
# read the RPC from raytheon file
print("Reading RPC from Raytheon file: {}".format(args_raytheon_rpc))
model = raytheon_rpc.read_raytheon_rpc_file(args_raytheon_rpc)
else:
model = args_raytheon_rpc
else:
# read the RPC from RPC Metadata in the image file
print("Reading RPC Metadata from {}".format(args_source_image))
Expand Down Expand Up @@ -157,12 +162,13 @@ def orthorectify(args_source_image, args_dsm, args_destination_image,
print("Driver {} does not supports Create().".format(driver))
return ERROR

# convert coordinates to Long/Lat
srs = osr.SpatialReference(wkt=projection)
proj_srs = srs.ExportToProj4()
inProj = pyproj.Proj(proj_srs)
outProj = pyproj.Proj('+proj=longlat +datum=WGS84')
arrayX, arrayY = pyproj.transform(inProj, outProj, arrayX, arrayY)
if (args_convert_to_latlon):
# convert coordinates to Long/Lat
srs = osr.SpatialReference(wkt=projection)
proj_srs = srs.ExportToProj4()
inProj = pyproj.Proj(proj_srs)
outProj = pyproj.Proj('+proj=longlat +datum=WGS84')
arrayX, arrayY = pyproj.transform(inProj, outProj, arrayX, arrayY)

# Sort the points by height so that higher points project last
if (args_occlusion_thresh > 0):
Expand Down Expand Up @@ -240,10 +246,11 @@ def orthorectify(args_source_image, args_dsm, args_destination_image,
print("Processing band {} ...".format(bandIndex))
sourceBand = sourceImage.GetRasterBand(bandIndex)
nodata_value = sourceBand.GetNoDataValue()
# for now use zero as a no-data value if one is not specified
# for now use -9999 as a no-data value if one is not specified
# it would probably be better to add a mask (alpha) band instead
if nodata_value is None:
nodata_value = 0
nodata_value = -9999
print("nodata: {}".format(nodata_value))
if numpy.any(cropSize < 1):
# read one value for data type
sourceRaster = sourceBand.ReadAsArray(
Expand Down
50 changes: 36 additions & 14 deletions tools/buildings_to_dsm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,10 +197,10 @@ def main(args):
buildingsScalarRange = p2cBuildings.GetOutput().GetCellData().GetScalars().GetRange()

if (args.debug):
polyWriter = vtk.vtkXMLPolyDataWriter()
polyWriter.SetFileName("p2c.vtp")
polyWriter.SetInputConnection(p2cBuildings.GetOutputPort())
polyWriter.Write()
writer = vtk.vtkXMLPolyDataWriter()
writer.SetFileName("p2c.vtp")
writer.SetInputConnection(p2cBuildings.GetOutputPort())
writer.Write()

buildingsMapper = vtk.vtkPolyDataMapper()
buildingsMapper.SetInputDataObject(p2cBuildings.GetOutput())
Expand All @@ -217,7 +217,26 @@ def main(args):
dtmReader = vtk.vtkGDALRasterReader()
dtmReader.SetFileName(args.input_dtm)
dtmReader.Update()
dtmVtk = dtmReader.GetOutput()
dtmVtkCell = dtmReader.GetOutput()

# convert from cell to point data
origin = dtmVtkCell.GetOrigin()
spacing = dtmVtkCell.GetSpacing()
dims = list(dtmVtkCell.GetDimensions())
dims[0] = dims[0] - 1
dims[1] = dims[1] - 1
dtmVtk = vtk.vtkUniformGrid()
dtmVtk.SetDimensions(dims)
dtmVtk.SetOrigin(origin)
dtmVtk.SetSpacing(spacing)
data = dtmVtkCell.GetCellData().GetScalars()
dtmVtk.GetPointData().SetScalars(data)

if (args.debug):
writer = vtk.vtkXMLImageDataWriter()
writer.SetFileName("dtm.vti")
writer.SetInputDataObject(dtmVtk)
writer.Write()

# Convert the terrain into a polydata.
surface = vtk.vtkImageDataGeometryFilter()
Expand All @@ -236,6 +255,12 @@ def main(args):
warp.Update()
dsmScalarRange = warp.GetOutput().GetPointData().GetScalars().GetRange()

if (args.debug):
writer = vtk.vtkXMLPolyDataWriter()
writer.SetFileName("dtm_warped.vtp")
writer.SetInputConnection(warp.GetOutputPort())
writer.Write()

dtmMapper = vtk.vtkPolyDataMapper()
dtmMapper.SetInputConnection(warp.GetOutputPort())
dtmActor = vtk.vtkActor()
Expand Down Expand Up @@ -264,6 +289,10 @@ def main(args):
ren.RemoveActor(dtmActor)

renWin.Render()
# iren = vtk.vtkRenderWindowInteractor()
# iren.SetRenderWindow(renWin)
# iren.Start();

windowToImageFilter = vtk.vtkWindowToImageFilter()
windowToImageFilter.SetInput(renWin)
windowToImageFilter.SetInputBufferTypeToRGBA()
Expand Down Expand Up @@ -310,15 +339,8 @@ def main(args):
valuePass.ReleaseGraphicsResources(renWin)

print("Writing the DSM ...")
elevationFlat = numpy_support.vtk_to_numpy(elevationFlatVtk)
# VTK X,Y corresponds to numpy cols,rows. VTK stores arrays
# in Fortran order.
elevationTranspose = numpy.reshape(
elevationFlat, [dtm.RasterXSize, dtm.RasterYSize], "F")
# changes from cols, rows to rows,cols.
elevation = numpy.transpose(elevationTranspose)
# numpy rows increase as you go down, Y for VTK images increases as you go up
elevation = numpy.flip(elevation, 0)
elevation = gdal_utils.vtk_to_numpy_order(elevationFlatVtk,
[dtm.RasterXSize, dtm.RasterYSize])
if args.buildings_only:
dsmElevation = elevation
else:
Expand Down
Loading