使用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，16U\16S，32U\32S，32F和64F。如果圖像數據是普通的數據（這裡指的是除復數圖像之外），那麼這五個函數完全可以滿足所有的需求。但是如果圖像是復數圖像（復數圖像就是圖像的像素值是由復數組成的，比如通過傅立葉變換後的圖像），這時，上面的五個函數就有點愛莫能助了，就需要上面的第六個函數來出場救急了。

第六個函數的聲明和C\C++的聲明完全一致，使用C\C++的同學對這個接口應該比較熟悉。這個函數可以替代上面五個函數來使用，但是比較麻煩的時，就是傳入的參數會比較多，首先看一個例子。使用普通的方式來讀取一個單波段8U的數據（Landsat5的一個數據）。從數據中讀取100×100大小的數據，直接使用ReadRaster的Byte接口。

static void ReadRaster()
{
	string strFile = @"D:\Data\landsat\etp139r26_5t19900902\p139r26_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\n");
}

下面看看如何使用第六個接口來完成上面的功能，代碼如下：

static void ReadRaster()
{
	string strFile = @"D:\Data\landsat\etp139r26_5t19900902\p139r26_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\n");
}

從上面的代碼中可以看出，使用第六個函數的時候，需要一個IntPtr的類型，使用Marshal類中的AllocCoTaskMem函數來進行分配內存空間（不知道C#中的說法是啥，按照C\C++中應該就是分配空間）。需要注意的是，分配的大小是按照字節（Byte）為單位的。接下來在調用ReadRaster時需要指定數據類型，調用之後，像素值應該都存儲在這個IntPtr中了，然後再使用Marshal類中的Copy函數將IntPtr中的數據拷貝到Byte數組data中即可。

通過上面的代碼，我們大致知道第六個ReadRaster的大致用法，接下來我們通過這個方法來讀取一個復數類型的圖像。沒有復數類型的圖像，可以使用ENVI隨便打開一個圖像，然後在Transform菜單下有個FFT，將打開的數據進行傅立葉變換，輸出的結果就是一個復數圖像。這裡我將上面使用的Landsat的數據使用FFT轉了下。使用下面的代碼進行打開。

static void ReadRaster()
{
	string strFile = @"D:\Data\landsat\etp139r26_5t19900902\p139r26_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\n");
}

在上面的代碼中，一個復數由實部和虛部組成，所以一個像素值就對應於普通圖像的2個像素值，故讀取100×100大小的數據就需要分配普通圖像2倍的大小。除此之外與讀取普通圖像一樣，不過在最後獲取像素值的時候，復數圖像的像素值是順序存儲的，即實部,虛部,實部,虛部...這樣的順序來進行存儲。上面的代碼運行的前五個像素值如下圖所示：

通過ENVI軟件查看該圖像的前五個像素值如下圖所示（圖中紅色框中就是像素值），對比出來可以看到讀取出來的實部和虛部完全一樣。<喎?http://www.Bkjia.com/kf/ware/vc/" target="_blank" class="keylink">vcD48cD48aW1nIHNyYz0="http://www.2cto.com/uploadfile/Collfiles/20140219/20140219093641418.jpg" alt="" />