This project is read-only.

Work with LiDAR

Nov 11, 2010 at 10:08 PM

Dear contributors,

This is Ping, I would like to add some functions related to LiDAR data for DotSpatial, the function including reading LAS, writing LAS, display LiDAR points, watershed analysis based on LiDAR points.etc...

Is there a team that I can join in for learning and discussing this issue?

I will be appreciated if you can give me some information for LiDAR in DotSpatial.

Ping

Jan 9, 2011 at 8:38 PM

What do you have working currently for this functionality?  I've been wanting to add a LAS reader class soon as well.

Jan 10, 2011 at 4:15 PM
This is Ping Yang, I had done some research on this function, I
attached a source file in this email, I hope it can be helpful.
Jan 10, 2011 at 5:57 PM
Edited Jan 10, 2011 at 6:07 PM

Did you?  I couldn't see anything attached (or a way to do so) or a link to code or source files.

Maybe you could just briefly describe what you've done so far or a link to source

Jan 10, 2011 at 6:06 PM
Now I paste the source code as follows(it is in C#):
********************************************************************
using System;
using System.Collections.Generic;
using System.Runtime.InteropServices;
using System.Text;
using System.IO;
using System.Drawing;

namespace LasLibrary
{
[StructLayout(LayoutKind.Sequential, Pack = 1,CharSet= CharSet.Ansi)]
public struct LAS_FILE_HEADER
{
public Int32 fileSignature; //“LASF”
public UInt16 fileSourceID;
public UInt16 globalEncoding;
public UInt32 projectIDGUIDdata1;
public UInt16 projectIDGUIDdata2;
public UInt16 projectIDGUIDdata3;
[MarshalAs(UnmanagedType.ByValArray, SizeConst = 8)]
public Byte[] projectIDGUIDdata4;
public Byte versionMajor;
public Byte versionMinor;
[MarshalAs(UnmanagedType.ByValTStr, SizeConst = 32)]
public String systemIdentifier;
[MarshalAs(UnmanagedType.ByValTStr, SizeConst = 32)]
public string generatingSoftware;
public UInt16 fileCreationDayofYear;
public UInt16 fileCreationYear;
public UInt16 headerSize;
public UInt32 offsetToPointData;
public UInt32 numberOfVariableLengthRecords;
public Byte pointDataFormatID; // (0-99 for spec)
public UInt16 pointDataRecordLength;
public UInt32 numberOfPointRecords;
[MarshalAs(UnmanagedType.ByValArray, SizeConst = 5)]
public UInt32[] numberOfPointsByReturn;
public Double xScaleFactor;
public Double yScaleFactor;
public Double zScaleFactor;
public Double xOffset;
public Double yOffset;
public Double zOffset;
public Double maxX;
public Double minX;
public Double maxY;
public Double minY;
public Double maxZ;
public Double minZ;
}

[StructLayout(LayoutKind.Sequential, Pack = 1, CharSet = CharSet.Ansi)]
public struct LAS_POINT_FORMAT_0
{
public Int32 X;
public Int32 Y;
public Int32 Z;
public UInt16 intensity;
public Byte pntFlags; // BitArray of Return Number 0:2, Number of Returns 3:5, Scan Direction Flag 6:, Edge of Flight Line 7:
public Byte classification;
public SByte scanAngleRank;
public Byte userData;
public UInt16 pointSrcID;
}

[StructLayout(LayoutKind.Sequential, Pack = 1, CharSet = CharSet.Ansi)]
public struct LAS_POINT_FORMAT_1
{
public Int32 X;
public Int32 Y;
public Int32 Z;
public UInt16 intensity;
public Byte pntFlags; // BitArray of Return Number 0:2, Number of Returns 3:5, Scan Direction Flag 6:, Edge of Flight Line 7:
public Byte classification;
public SByte scanAngleRank;
public Byte userData;
public UInt16 pointSrcID;
public Double GPSTime;
}

[StructLayout(LayoutKind.Sequential, Pack = 1, CharSet = CharSet.Ansi)]
public struct LAS_POINT_FORMAT_2
{
public Int32 X;
public Int32 Y;
public Int32 Z;
public UInt16 intensity;
public Byte pntFlags; // BitArray of Return Number 0:2, Number of Returns 3:5, Scan Direction Flag 6:, Edge of Flight Line 7:
public Byte classification;
public SByte scanAngleRank;
public Byte userData;
public UInt16 pointSrcID;
public UInt16 red;
public UInt16 green;
public UInt16 blue;
}

[StructLayout(LayoutKind.Sequential, Pack = 1, CharSet = CharSet.Ansi)]
public struct LAS_POINT_FORMAT_3
{
public Int32 X;
public Int32 Y;
public Int32 Z;
public UInt16 intensity;
public Byte pntFlags; // BitArray of Return Number 0:2, Number of Returns 3:5, Scan Direction Flag 6:, Edge of Flight Line 7:
public Byte classification;
public SByte scanAngleRank;
public Byte userData;
public UInt16 pointSrcID;
public Double GPSTime;
public UInt16 red;
public UInt16 green;
public UInt16 blue;
}

public struct POINT_RECORD
{
public bool skipPnt;
public Double X;
public Double Y;
public Double Z;
public UInt16 intensity;
public Color pntColor;
}

public class LasReader
{
UInt32 _curRec;
FileStream fs;
LAS_FILE_HEADER fh;
int fhStructSize;
int lpf0Size;
int lpf1Size;
int lpf2Size;
int lpf3Size;
bool pntRecLenMatchType; //Does the pnt record length match the current type of record length
Int32 recAdjustment; //How many bytes to move when an adjustment is required

public LasReader(String fileName)
{
//Is file name valid?
//Does file exist
if (!File.Exists(fileName))
throw (new FileNotFoundException());

lpf0Size = Marshal.SizeOf(typeof(LAS_POINT_FORMAT_0));
lpf1Size = Marshal.SizeOf(typeof(LAS_POINT_FORMAT_1));
lpf2Size = Marshal.SizeOf(typeof(LAS_POINT_FORMAT_2));
lpf3Size = Marshal.SizeOf(typeof(LAS_POINT_FORMAT_3));
fhStructSize = Marshal.SizeOf(typeof(LAS_FILE_HEADER));
fs = new FileStream(fileName, FileMode.Open, FileAccess.Read);
ReadFileHeader();

// Does the point record length math the header size for the type
pntRecLenMatchType = false;
switch (fh.pointDataFormatID)
{
case 0:
if (fh.pointDataRecordLength == lpf0Size)
pntRecLenMatchType = true;
else
recAdjustment = fh.pointDataRecordLength - lpf0Size;
break;
case 1:
if (fh.pointDataRecordLength == lpf1Size)
pntRecLenMatchType = true;
else
recAdjustment = fh.pointDataRecordLength - lpf1Size;
break;
case 2:
if (fh.pointDataRecordLength == lpf2Size)
pntRecLenMatchType = true;
else
recAdjustment = fh.pointDataRecordLength - lpf2Size;
break;
case 3:
if (fh.pointDataRecordLength == lpf3Size)
pntRecLenMatchType = true;
else
recAdjustment = fh.pointDataRecordLength - lpf3Size;
break;
}

// Set the file to the first record
fs.Seek(fh.offsetToPointData, SeekOrigin.Begin);
_curRec = 0;

}

public bool ParsableFormat
{
get
{
if (fh.pointDataFormatID == 0 ||
fh.pointDataFormatID == 1 ||
fh.pointDataFormatID == 2 ||
fh.pointDataFormatID == 3)
return true;
else
return false;
}
}

public bool HasColor
{
get
{
if (fh.pointDataFormatID == 2 ||
fh.pointDataFormatID == 3)
return true;
else
return false;
}
}

public UInt32 RecordCount
{
get{ return fh.numberOfPointRecords;}
}
public POINT_RECORD ReadPnt(UInt32 i)
//public POINT_RECORD this[UInt32 i]
{
// get
// {
//Starting at the
if (_curRec != i)
{
_curRec = i + 1;
fs.Seek(fh.pointDataRecordLength * i + fh.offsetToPointData, SeekOrigin.Begin);
}
else
_curRec++;
return ReadCurrentPoint();
// }
}
public POINT_RECORD ReadPnt(UInt32 i, UInt16 classification)
//public POINT_RECORD this[UInt32 i, UInt16 classification]
{
// get
// {
//Starting at the
if (_curRec != i)
{
_curRec = i + 1;
fs.Seek(fh.pointDataRecordLength * i + fh.offsetToPointData, SeekOrigin.Begin);
}
else
_curRec++;
return ReadCurrentPoint(classification);
// }
}

~LasReader()
{
fs.Close();
}

private void ReadFileHeader()
{
byte[] buff = new byte[fhStructSize];
Int32 sr = fs.Read(buff, 0, fhStructSize);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
fh = (LAS_FILE_HEADER)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_FILE_HEADER));
handle.Free();
}

private POINT_RECORD ReadCurrentPoint()
{
POINT_RECORD pnt = new POINT_RECORD();
// What type of point record are we reading
switch (fh.pointDataFormatID)
{
case 0:
{
byte[] buff = new byte[lpf0Size];
Int32 sr = fs.Read(buff, 0, lpf0Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_0 pntRec = (LAS_POINT_FORMAT_0)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_0));
handle.Free();
if (TestPntSkip(pntRec.classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
}
break;
case 1:
{
byte[] buff = new byte[lpf1Size];
Int32 sr = fs.Read(buff, 0, lpf1Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_1 pntRec = (LAS_POINT_FORMAT_1)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_1));
handle.Free();
if (TestPntSkip(pntRec.classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
}
break;
case 2:
{
byte[] buff = new byte[lpf2Size];
Int32 sr = fs.Read(buff, 0, lpf2Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_2 pntRec = (LAS_POINT_FORMAT_2)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_2));
handle.Free();
if (TestPntSkip(pntRec.classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
pnt.pntColor = Color.FromArgb(pntRec.red, pntRec.green, pntRec.blue);
}
break;
case 3:
{
byte[] buff = new byte[lpf3Size];
Int32 sr = fs.Read(buff, 0, lpf3Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_3 pntRec = (LAS_POINT_FORMAT_3)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_3));
handle.Free();
if (TestPntSkip(pntRec.classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
pnt.pntColor = Color.FromArgb(pntRec.red, pntRec.green, pntRec.blue);
}
break;
}
if (!pntRecLenMatchType)
fs.Seek(recAdjustment, SeekOrigin.Current);

return pnt;
}

private POINT_RECORD ReadCurrentPoint(UInt16 classification)
{
POINT_RECORD pnt = new POINT_RECORD();
// What type of point record are we reading
switch (fh.pointDataFormatID)
{
case 0:
{
byte[] buff = new byte[lpf0Size];
Int32 sr = fs.Read(buff, 0, lpf0Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_0 pntRec = (LAS_POINT_FORMAT_0)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_0));
handle.Free();
if (TestPntSkip(pntRec.classification,classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
}
break;
case 1:
{
byte[] buff = new byte[lpf1Size];
Int32 sr = fs.Read(buff, 0, lpf1Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_1 pntRec = (LAS_POINT_FORMAT_1)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_1));
handle.Free();
if (TestPntSkip(pntRec.classification, classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
}
break;
case 2:
{
byte[] buff = new byte[lpf2Size];
Int32 sr = fs.Read(buff, 0, lpf2Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_2 pntRec = (LAS_POINT_FORMAT_2)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_2));
handle.Free();
if (TestPntSkip(pntRec.classification, classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
pnt.pntColor = Color.FromArgb(pntRec.red, pntRec.green, pntRec.blue);
}
break;
case 3:
{
byte[] buff = new byte[lpf3Size];
Int32 sr = fs.Read(buff, 0, lpf3Size);
GCHandle handle = GCHandle.Alloc(buff, GCHandleType.Pinned);
LAS_POINT_FORMAT_3 pntRec = (LAS_POINT_FORMAT_3)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(LAS_POINT_FORMAT_3));
handle.Free();
if (TestPntSkip(pntRec.classification, classification))
{
pnt.skipPnt = true;
break;
}
pnt.skipPnt = false;
pnt.X = pntRec.X * fh.xScaleFactor + fh.xOffset;
pnt.Y = pntRec.Y * fh.yScaleFactor + fh.yOffset;
pnt.Z = pntRec.Z * fh.zScaleFactor + fh.zOffset;
pnt.intensity = pntRec.intensity;
pnt.pntColor = Color.FromArgb(pntRec.red, pntRec.green, pntRec.blue);
}
break;
}
if (!pntRecLenMatchType)
fs.Seek(recAdjustment, SeekOrigin.Current);

return pnt;
}
private bool TestPntSkip(byte testField)
{
if ((testField & 0x80) == 0x80)
return true;
else
return false;

}
private bool TestPntSkip(byte testField, UInt16 classification)
{
if ((testField & 0x80) == 0x80)
return true;
else
{
bool bRetVal = true; // By default we will skip
switch ((testField & 0x1F))
{
case 0: // Unclassified
case 1: // Unclassified
if ((classification & 0x01) == 0x001 || classification == 0)
bRetVal = false;
break;
case 2: // Ground
if ((classification & 0x002) == 0x002)
bRetVal = false;
break;
case 3: // Low Vegetation
if ((classification & 0x004) == 0x004)
bRetVal = false;
break;
case 4: // Medium Vegetation
if ((classification & 0x008) == 0x008)
bRetVal = false;
break;
case 5: // High Vegetation
if ((classification & 0x010) == 0x010)
bRetVal = false;
break;
case 6: // Buildings
if ((classification & 0x020) == 0x020)
bRetVal = false;
break;
case 7: // Low Point(noise)
if ((classification & 0x040) == 0x040)
bRetVal = false;
break;
case 8: // Model Key-point
if ((classification & 0x080) == 0x080)
bRetVal = false;
break;
case 9: // Water
if ((classification & 0x100) == 0x100)
bRetVal = false;
break;
default:
break;
}
return bRetVal;
}
}
}
}

**********************************************************
Here used Marshal for management of the memory.

Ping Yang

Jan 25, 2011 at 10:40 PM

It looks sort of similar to what I have...can I make a couple of suggestions?

1) Use boring, primitive BinaryReader functions to fill structs (ReadInt32, ReadDouble, etc), vs the Marshal.PtrToStructure approach, which is cooler though :-)  You should notice a significant performance increase.  But maybe you've experienced this differently?

2) Maybe use a factory pattern, even a 'simple' factory would get rid of all those switch statements, duplicate less code and allow for a easier time of upgrading the library once LAS gets versioned again.

Jan 26, 2011 at 6:15 PM

Good suggestion!

I would more like to see the factory version, would you please upgrade it?

Jan 26, 2011 at 8:04 PM

Before I go writing code, can I ask how do you usually use the lidar points - are they indexed, stored in a data structure or a list, or do you iterate through the contents on-disk? Do you have a sample 'main' method or something as an example of how you'd use the class?

Also does this live in Dotspatial or on Codeplex somewhere?  Looks good.

Jan 26, 2011 at 8:18 PM

Actually I am using a list to store the data now, however, I am thinking an index method to manage the points in 3D representation. you can code a main method or you can use this class simply by creating a bottom event in a form with the following code: 

 String LASfile = "C:\\alaska24.las";          

 LasLibrary.LasReader Reader = new LasReader(LASfile); 

 

Jan 27, 2011 at 5:35 PM

OK I looked a bit closer at the code and it looks like the ReadPoint function doesn't return the specific type of LAS record being used, but a lightweight 'point_record' struct, is this correct?