使用C#版本GDAL读取复数图像

GDAL的C#版本虽然在很多算法接口没有导出，但是在读写数据中的接口基本上都是完全导出了。使用ReadRaster和WriteRaster方法来进行读写，同时对这两个方法进行了重载，对于常用的数据类型可以不用指定数据类型直接进行读取即可。但是对于复数类型就有点复杂了。下面就针对GDAL如何来读取复数数据来进行一个简单的说明。

我们知道，在使用GDAL读取数据的时候使用的是ReadRaster这个函数，这个函数重载了6个，函数声明分别如下，以Dataset的ReadRaster为例，Band类中的ReadRaster与之类似。

OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Byte[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Int16[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Single[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Double[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.IntPtr, System.Int32, System.Int32, OSGeo.GDAL.DataType, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);

前五个函数的声明基本上完全一致，除了第五个参数，分别是Byte[]，Int16[]，Int32[]，Single[]，Double[]。分别对应的是将图像中的像素值使用8U，16U16S，32U32S，32F和64F。如果图像数据是普通的数据（这里指的是除复数图像之外），那么这五个函数完全可以满足所有的需求。但是如果图像是复数图像（复数图像就是图像的像素值是由复数组成的，比如通过傅立叶变换后的图像），这时，上面的五个函数就有点爱莫能助了，就需要上面的第六个函数来出场救急了。

第六个函数的声明和CC++的声明完全一致，使用CC++的同学对这个接口应该比较熟悉。这个函数可以替代上面五个函数来使用，但是比较麻烦的时，就是传入的参数会比较多，首先看一个例子。使用普通的方式来读取一个单波段8U的数据（Landsat5的一个数据）。从数据中读取100×100大小的数据，直接使用ReadRaster的Byte接口。

static void ReadRaster()
{
	string strFile = @"D:Datalandsatetp139r26_5t19900902p139r26_5t19900902_nn1.tif";

	OSGeo.GDAL.Gdal.AllRegister();
	OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(strFile, OSGeo.GDAL.Access.GA_ReadOnly);
   
	int []bandmap = new int [1];
	bandmap[0] = 1;

	//读取Byte类型
	Byte[] data = new Byte[100 * 100];
	OSGeo.GDAL.CPLErr err = ds.ReadRaster(3000, 2000, 100, 100, data, 100, 100, 1, bandmap, 0, 0, 0);
	
	//将读取的数据输出
	for (int i = 0; i < 100; i++ )
	{
		for (int j = 0; j < 100; j++)
		{
			Console.WriteLine("{0}", ddata[i*100+j]);
		}
	}
	Console.WriteLine("done
");
}

下面看看如何使用第六个接口来完成上面的功能，代码如下：

static void ReadRaster()
{
	string strFile = @"D:Datalandsatetp139r26_5t19900902p139r26_5t19900902_nn1.tif";

	OSGeo.GDAL.Gdal.AllRegister();
	OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(strFile, OSGeo.GDAL.Access.GA_ReadOnly);
   
	int []bandmap = new int [1];
	bandmap[0] = 1;

	//读取Byte类型
	Byte[] data = new Byte[100 * 100];
	
	//构造一个IntPtr，大小为100×100个Byte
	IntPtr ptr = Marshal.AllocCoTaskMem(data.Length);
	//需要指明读取的数据类型
	OSGeo.GDAL.CPLErr err = ds.ReadRaster(3000, 2000, 100, 100, ptr, 100, 100, OSGeo.GDAL.DataType.GDT_Byte, 1, bandmap, 0, 0, 0);
	//将读取到的数据拷贝到data中
	Marshal.Copy(ptr, data, 0, data.Length);

	//将读取的数据输出
	for (int i = 0; i < 100; i++ )
	{
		for (int j = 0; j < 100; j++)
		{
			Console.WriteLine("{0}", ddata[i*100+j]);
		}
	}
	Console.WriteLine("done
");
}

从上面的代码中可以看出，使用第六个函数的时候，需要一个IntPtr的类型，使用Marshal类中的AllocCoTaskMem函数来进行分配内存空间（不知道C#中的说法是啥，按照CC++中应该就是分配空间）。需要注意的是，分配的大小是按照字节（Byte）为单位的。接下来在调用ReadRaster时需要指定数据类型，调用之后，像素值应该都存储在这个IntPtr中了，然后再使用Marshal类中的Copy函数将IntPtr中的数据拷贝到Byte数组data中即可。

通过上面的代码，我们大致知道第六个ReadRaster的大致用法，接下来我们通过这个方法来读取一个复数类型的图像。没有复数类型的图像，可以使用ENVI随便打开一个图像，然后在Transform菜单下有个FFT，将打开的数据进行傅立叶变换，输出的结果就是一个复数图像。这里我将上面使用的Landsat的数据使用FFT转了下。使用下面的代码进行打开。

static void ReadRaster()
{
	string strFile = @"D:Datalandsatetp139r26_5t19900902p139r26_5t19900902_nn1_fft.tif";

	OSGeo.GDAL.Gdal.AllRegister();
	OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(strFile, OSGeo.GDAL.Access.GA_ReadOnly);
   
	int []bandmap = new int [1];
	bandmap[0] = 1;


	//读取复数类型，由于一个复数由两个数据组成，所以读取100×100的像元需要2倍普通图像的大小
	double[] data = new double[100 * 100 * 2];
	
	//构造一个IntPtr，大小为100×100×2个double，一个double为4个byte，所以下面的长度要乘以8
	IntPtr ptr = Marshal.AllocCoTaskMem(data.Length*8);
	//需要指明读取的数据类型
	OSGeo.GDAL.CPLErr err = ds.ReadRaster(0, 0, 100, 100, ptr, 100, 100, OSGeo.GDAL.DataType.GDT_CFloat64, 1, bandmap, 0, 0, 0);
	//将读取到的数据拷贝到data中
	Marshal.Copy(ptr, data, 0, data.Length*8);

	//将读取的数据输出
	int in1 = 0;
	for (int i = 0; i < 100; i++ )
	{
		for (int j = 0; j < 100; j++)
		{
			Console.WriteLine("({0}*i+{1})", ddata[in1++], ddata[in1++]);
		}
	}
	Console.WriteLine("done
");
}

在上面的代码中，一个复数由实部和虚部组成，所以一个像素值就对应于普通图像的2个像素值，故读取100×100大小的数据就需要分配普通图像2倍的大小。除此之外与读取普通图像一样，不过在最后获取像素值的时候，复数图像的像素值是顺序存储的，即实部,虚部,实部,虚部...这样的顺序来进行存储。上面的代码运行的前五个像素值如下图所示：

通过ENVI软件查看该图像的前五个像素值如下图所示（图中红色框中就是像素值），对比出来可以看到读取出来的实部和虚部完全一样。