LAS点云的读写以及投影信息的读写（LASlib）

本文使用LASlib库完成点云及投影信息的读写。具体安装方式可自行网上搜索，推荐使用vcpkg安装，非常方便。

对于库的使用下面直接上代码：

头文件：

#include <LASlib/lasreader.hpp>
#include <LASlib/laswriter.hpp>

1、las文件的读

先是LASreadOpener对象构造，设置文件名，然后进行点云的打开。

    LASreadOpener lasreadopener;
    lasreadopener.set_file_name(input_file.c_str());
    LASreader* lasreader = lasreadopener.open();

读取las头文件信息：

    //打印头文件信息
    F64 x_offset = lasreader->header.x_offset;
	F64 y_offset = lasreader->header.y_offset;
	F64 z_offset = lasreader->header.z_offset;
	F64 x_scale_factor = lasreader->header.x_scale_factor;
	F64 y_scale_factor = lasreader->header.y_scale_factor;
	F64 z_scale_factor = lasreader->header.z_scale_factor;
	U8 point_data_format = lasreader->header.point_data_format;
	U16 point_data_record_length = lasreader->header.point_data_record_length;
	
	U16 file_source_ID = lasreader->header.file_source_ID;
	U16 global_encoding = lasreader->header.global_encoding;
	//读取.las文件的版本号
	U8 version_major = lasreader->header.version_major;
	U8 version_minor = lasreader->header.version_minor;

查看las文件的投影信息，可以参考该博客。

    string input_file = argv[1];

	LASreadOpener lasreadopener;
	lasreadopener.set_file_name(input_file.c_str());
	LASreader* lasreader = lasreadopener.open();

	printf("system_identifier : %s\n", lasreader->header.system_identifier);
	printf("generating_software : %s\n", lasreader->header.generating_software);

	char* ds = (char*)lasreader->header.vlrs->description;
	printf("vlrs description: %s\n", ds);

	char* ui = (char*)lasreader->header.vlrs->user_id;
	printf("user_id: %s\n", ui);
	printf("point_data_format:%d\n", lasreader->header.point_data_format);

	printf("minX:%f maxX:%f\n", lasreader->header.min_x, lasreader->header.max_x);
	printf("minY:%f maxY:%f\n", lasreader->header.min_y, lasreader->header.max_y);
	printf("minZ:%f maxZ:%f\n", lasreader->header.min_z, lasreader->header.max_z);

	printf("number_of_point_records:%d\n", lasreader->header.number_of_point_records);

	//顺序读取Geokeys
	int numof_keys2 = lasreader->header.vlr_geo_keys->number_of_keys;
	cout << "numof_keys2:" << numof_keys2 << endl;
	for (int i = 0; i < numof_keys2; i++)
	{
		//顺序读取GeoKey内容
		U16 id = (lasreader->header.vlr_geo_key_entries + i)->key_id;
		U16 Location = (lasreader->header.vlr_geo_key_entries + i)->tiff_tag_location;
		U16 count = (lasreader->header.vlr_geo_key_entries + i)->count;
		U16 value_offset = (lasreader->header.vlr_geo_key_entries + i)->value_offset;

		//vlr_geo_double_params
		F64 value;
		if (Location == 34736)
		{
			value = lasreader->header.vlr_geo_double_params[value_offset];
			printf("keysID:%d ,Location:%d,Count: %d ,Value_offset:%d ,Value:%f \n", id, Location, count, value_offset, value);
		}
		//vlr_geo_ascii_params
		else if (Location == 34737)
		{
			CHAR value_char[32];
			for (int j = 0; j < count; j++)
			{
				value_char[j] = lasreader->header.vlr_geo_ascii_params[value_offset + j];
			}
			printf("keysID:%d ,Location:%d,Count: %d ,Value_offset:%d ,ValueChar:%s \n", id, Location, count, value_offset, value_char);
		}
		else
			printf("keysID:%d ,Location:%d,Count: %d ,Value_offset:%d\n", id, Location, count, value_offset);
	}

打印结果如下：

system_identifier : libLAS
generating_software : libLAS 1.8.1
vlrs description: GeoTIFF GeoKeyDirectoryTag
user_id: LASF_Projection
point_data_format:3
minX:772891.614624 maxX:774110.863464
minY:2572529.591858 maxY:2573509.307495
minZ:-468.167542 maxZ:190.500195
number_of_point_records:36862258
numof_keys2:7
keysID:1024 ,Location:0,Count: 1 ,Value_offset:1
keysID:1025 ,Location:0,Count: 1 ,Value_offset:1
keysID:1026 ,Location:34737,Count: 22 ,Value_offset:0 ,ValueChar:WGS 84 / UTM zone 49N|
keysID:2049 ,Location:34737,Count: 7 ,Value_offset:22 ,ValueChar:WGS 84|/ UTM zone 49N|
keysID:2054 ,Location:0,Count: 1 ,Value_offset:9102
keysID:3072 ,Location:0,Count: 1 ,Value_offset:32649
keysID:3076 ,Location:0,Count: 1 ,Value_offset:9001

2、las文件的写

las文件的写使用的是LASwriteOpener，和LASreadOpener 类似，先构造LASwriteOpener对象，然后再设置输出路径。

    LASwriteOpener laswriteopener;
	laswriteopener.set_file_name(out_laz.c_str());

下面是设置las的头文件信息，代码中的具体参数值要根据实际情况数据来设定。

    LASheader las_header_writer;
	las_header_writer.global_encoding = global_encoding;
	las_header_writer.x_scale_factor = resolution;
	las_header_writer.y_scale_factor = resolution;
	las_header_writer.z_scale_factor = resolution;
	las_header_writer.x_offset = x_offset;
	las_header_writer.y_offset = y_offset;
	las_header_writer.z_offset = z_offset;
	las_header_writer.point_data_format = point_data_format;
	las_header_writer.point_data_record_length = point_data_record_length;
	las_header_writer.version_minor = version_minor;
	las_header_writer.version_major = version_major;
	las_header_writer.file_source_ID = file_source_ID;
	
	char system_identifier[32] = "PDAL";
	char generating_software[32] = "PDAL 2.1.0";
	char user_id[16] = "LASF_Projection";
	char description[32] = " GeoTiff GeoKeyDirectoryTag";
	
	las_header_writer.vlrs = new LASvlr();
	strcpy(las_header_writer.system_identifier, system_identifier);
	strcpy(las_header_writer.generating_software, generating_software);
	strcpy(las_header_writer.vlrs->description, description);
	strcpy(las_header_writer.vlrs->user_id, user_id);

上面是设置las文件基本的头文件信息，如果需要增加其他的信息，比如添加las文件的投影信息，则需要在可变长记录区(VLRs)中添加。下面是添加las文件的投影信息的例子。代码中的key_id以及value_offset代表的意义可以在这里查看。

    LASvlr_key_entry geo_keys[5];
	geo_keys[0].key_id = 1024; // GTModelTypeGeoKey
	geo_keys[0].tiff_tag_location = 0;
	geo_keys[0].count = 1;
	geo_keys[0].value_offset = 1; // ModelTypeProjected

	geo_keys[1].key_id = 3072; // ProjectedCSTypeGeoKey
	geo_keys[1].tiff_tag_location = 0;
	geo_keys[1].count = 1;
	geo_keys[1].value_offset = epsg_code; // PCS_WGS84_UTM
	
	// horizontal units
	geo_keys[2].key_id = 3076; // ProjLinearUnitsGeoKey
	geo_keys[2].tiff_tag_location = 0;
	geo_keys[2].count = 1;
	geo_keys[2].value_offset = 9001; // meters

	// vertical units
	geo_keys[3].key_id = 4099; // VerticalUnitsGeoKey
	geo_keys[3].tiff_tag_location = 0;
	geo_keys[3].count = 1;
	geo_keys[3].value_offset = 9001; // meters

	// vertical datum
	geo_keys[4].key_id = 4096; // VerticalCSTypeGeoKey
	geo_keys[4].tiff_tag_location = 0;
	geo_keys[4].count = 1;
	geo_keys[4].value_offset = 5030; // WGS84

	las_header_writer.set_geo_keys(5, geo_keys);

接下来是构造LASwriter对象，写入点数据到该对象中。

    // open laswriter
	LASwriter* laswriter = laswriteopener.open(&las_header_writer);
	if (laswriter == 0)
	{
		printf("ERROR: could not open laswriter\n");
		return EXIT_FAILURE;
	}

	std::vector<LASpoint> vecLASpoint;
	LASpoint laspoint;
    // 初始化点的格式及一些属性
	laspoint.init(&las_header_writer, las_header_writer.point_data_format, las_header_writer.point_data_record_length, 0);

 通过遍历给LASwriter对象写入点数据：

    double maxX = -DBL_MAX, maxY = -DBL_MAX, maxZ = -DBL_MAX;  
    int i_points= 0;      
    for (int i = 0; i < 100000; ++i)
    {
	    laspoint.set_X(pointX[i]);
	    laspoint.set_Y(pointY[i]);
	    laspoint.set_Z(pointZ[i]);

	    laspoint.set_R(pointR[i]);
	    laspoint.set_G(pointG[i]);
	    laspoint.set_B(pointB[i]);

	    // write the point
	    laswriter->write_point(&laspoint);

	    // add it to the inventory
	    laswriter->update_inventory(&laspoint);

	    //range
	    double x = pointX[i];
	    double y = pointY[i];
	    double z = pointZ[i];
	    if (x < minX) minX = x;
	    if (x > maxX) maxX = x;
	    if (y < minY) minY = y;
	    if (y > maxY) maxY = y;
	    if (z < minZ) minZ = z;
	    if (z > maxZ) maxZ = z;
	    i_points++;
    }
    // update the boundary
    las_header_writer.set_bounding_box(minX, minY, minZ, maxX, maxY, maxZ);

    las_header_writer.number_of_point_records = i_points;

    // update the header
    laswriter->update_header(&las_header_writer, TRUE, TRUE);

    // close the writer
    laswriter->close();
    delete laswriter;

以上记录了las文件的读写，以及las的投影信息的读写，如有错误望指正。