Using Sun's denoising algorithm on topographic data

Introduction

Preparation

Windows users

Installing GDAL

• Install FWtools from http://fwtools.maptools.org/

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 users

Installing GDAL

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 users

Installing GDAL

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

• 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

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

• Manchester University School of Earth Atmospheric and Environmental Science
• Cardiff University Mesh Filtering Page
• Engineering and Physical Sciences Research Council