Manchester University

EPSRC

Using Sun's denoising algorithm on topographic data


Introduction  |  Preparation  |  Denoising rasters  |  Denoising point clouds  |  Reprojection  |  References  |  Links


Introduction

Sun et al, (2007) developed a feature-preserving mesh denoising algorithm that smooths the surfaces of computer models of three dimensional objects such as those used in computer-aided design and graphics.

Sun's denoising algorithm is freely available and easy to use.  The code now has native support for topographic data in projected coordinate systems, including both raster data (in ESRI ASCII Grid format) and point clouds (as XYZ coordinates).  This allows it to be used with digital elevation models, removing noise while preserving sharp features in areas such as mountain ranges.

This page provides instructions for using the algorithm on topographic data.  Most topographic data can be converted into a suitable format using the GDAL utilities.  These are a group of freely-available open-source programs that convert GIS data between different projection systems and formats (e.g. Arc, Erdas IMAGINE, GMT, Geotiff, KML).  The installation and use of the GDAL utilities on various platforms is also described.

Sun's denoising algorithm can also be run directly from within GRASS GIS using the r.denoise addon script.
Example: SRTM data from California
Iterative denoising of SRTM data

Download this example dataset: my_dem.tif

Preparation

Microsoft Windows

Windows users

Installing GDAL
The easiest way to install the GDAL utilities in Windows is via FWTools, which are a collection of freely-available, open source GIS progams including the GDAL utilities as well as projection tools and a simple GIS data viewer.
Installing Sun's denoising algorithm
  • Download and unzip WinMDenoise.zip from Cardiff University Mesh Filtering group website.  The Windows executable is MDenoise.exe.
  • Create a working directory on your computer and copy your data and MDenoise.exe into it.  

Linux

Linux users

Installing GDAL
The GDAL utilities are available in the repositories of many common distributions (e.g. Ubuntu) and can be installed via your package manager.  Further information, or instructions for compilation from source are available on the GDAL website.
Installing Sun's denoising algorithm
  • Download and unzip mdsource.zip from the Cardiff University Mesh Filtering group website.
  • Enter the mdenoise directory and compile the source code:

    g++ -o mdenoise mdenoise.cpp triangle.c

  • Link/copy mdenoise to your working directory, or to a directory in your PATH

Apple Mac

Apple Mac users

Installing GDAL
The GDAL utilities have been compiled and made available for download as binaries from the kyngchaos website.  Further information, or instructions for compilation from source are available on the GDAL website.
Installing Sun's denoising algorithm
  • Download and unzip mdsource.zip from the Cardiff University Mesh Filtering group website.
  • Enter the mdenoise directory and compile the source code:

    g++ -o mdenoise mdenoise.cpp triangle.c

  • Link/copy mdenoise to your working directory, or to a directory in your PATH

Denoising rasters

1) Checking the installation

Both GDAL and Sun's denoising algorithm are run from the command line.  Windows users should access the command line by running the FWTools shell.  This also gives access to the GDAL utilities.  Linux / Mac users can open a normal terminal.
  • Navigate to your working directory:

    cd "C:\Documents and Settings\Username\My Documents\working\"

  • Testrun the mdenoise program to see a list of options:

    MDenoise.exe

    (Note: Linux and Mac users substitute mdenoise for MDenoise.exe)

  • Test gdalinfo to see the range of formats available:

    gdalinfo --formats

2) Conversion to ESRI ASCII Grid format

  • Data should be in a projected coordinate system (e.g. UTM).  If this is not the case, read the reprojection section.
  • For processing with Sun's denoising algorithm, the data must be converted to ESRI ASCII grid format.  This is done with the gdal_translate utility:

    gdal_translate -of AAIGrid my_dem.tif my_dem.asc

  • This will create the files my_dem.asc with the data and my_dem.prj with the projection information.

3) Denoising

  • Sun's denoising algorithm is run from the command line. The amount of denoising is controlled by two parameters, which are set when the command is entered.  The -n flag controls the number of normal-updating iterations, and -t controls the threshold.  
  • More normal-updating iterations lead to greater smoothing, and higher threshold values preserve sharper features (but allow noise to pass if set too high).
  • Suitable parameter values depend on the data and the noise but are usually in the range n = 1-15 and t = 0.90-0.99.    
  • For SRTM data, such as the example dataset, n = 5, t = 0.99 gives good results.
  • Run mdenoise, specifying the settings that you require:

    MDenoise.exe -i my_dem.asc -n 5 -t 0.99 -o my_dem_DN.asc

    (Note: Linux and Mac users substitute mdenoise for MDenoise.exe)
  • Two new files, my_dem_DN.asc with the denoised data and my_dem_DN.prj with the projection information will be created.

4) Returning to the original format

  • The GDAL utilities are used to return the data to their original format:

    gdal_translate -of GTiff my_dem_DN.asc my_dem_DN.tif

  • Different format codes from the output of gdalinfo --formats can be substituted for 'GTiff'.
  • The denoised file can then be reopened with the GIS software.

Denoising point clouds

  • Point clouds should be supplied as ASCII text files, with the extension .xyz.  
  • The file must contain x, y, z data, separated by white space, in the first three columns.  Header lines (beginning with a # symbol) are ignored, as are additional columns (e.g. intensity, R, G, B, etc.).
  • Run mdenoise, specifying the settings that you require:

    MDenoise.exe -i my_points.xyz -n 5 -t 0.99 -o my_points_DN.xyz

    (Note: Linux and Mac users substitute mdenoise for MDenoise.exe)
  • During denoising, datapoints are moved in three dimensions.  If you want to limit the movements to just vertical changes, add the '-z' flag to your command.

Reprojection


Sun's denoising algorithm requires that the data are in a projected coordinate system.  If necessary, it is possible to reproject data using the GDAL utilities.  The example dataset is available in geographic coordinates (lat long): my_dem_ll.tif

1) Checking that the coordinate system is projected (metric)

  • Find out projection information: 

    gdalinfo my_dem_ll.tif
  • The resulting output should contain a section like this: 

    Coordinate System is:
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]]

  • If the result has UNITs of "degree", then you need to change to a projected coordinate system with units of metres. This can be done using the gdalwarp utility. 
  • Coordinate systems can be specified to the GDAL software using a short code (the EPSG code).  In this example, the current projection of the data is represented by the final EPSG code, 4326.

2) Reprojecting your data

  • You can find a suitable system for your data by searching on http://www.epsg-registry.org.  In the Type box, select "Coordinate Reference System", then "Projected CRS".  Enter the coordinates or name of your area and hit "Search".
  • Choose a coordinate system from the list produced.  Systems of WGS 84/UTM zone XXX type are suitable, and are compatible with most GPS systems.
  • Reproject your data using the gdalwarp utility, substituting in your own EPSG code:

    gdalwarp -t_srs "EPSG:32610" -rc my_dem_ll.tif my_dem_utm.tif

  • The reprojected data can now be denoised as before.

3) Returning to the original projection

  • The GDAL utilities are used to return the denoised data to their original format, using the original EPSG code.

    gdalwarp -t_srs "EPSG:4326" -rc my_dem_utm_DN.tif my_dem_ll_DN.tif

  • Point clouds can reprojected with the cs2cs program, which is part of the proj4 libraries and normally installed alongside GDAL (see instructions above).  Like GDAL, cs2cs can use EPSG codes.
  • Note that any form of reprojection of raster grids involves interpolation, rounding, and may introduce aliasing.  
  • In situations where aliasing is a problem, it can be avoided by exporting the raster as a point cloud, reprojecting, denoising the .xyz file with the '-z' flag, then replacing the old z values with the new ones.

References

  • Sun X, Rosin PL, Martin RR, Langbein FC (2007) Fast and Effective Feature-Preserving Mesh Denoising. IEEE Transactions on Visualisation and Computer Graphics, 13(5):925-938 doi:10.1109/TVCG.2007.1065

  • Stevenson JA, Sun X, Mitchell NC (2009) Despeckling SRTM and other topographic data with a denoising algorithm. Geomorphology doi:10.1016/j.geomorph.2009.07.006

Links


For comments / questions / further information, contact johnalexanderstevenson@yahoo.co.uk