C# GDAL将16位影像重采样2%线性拉伸并生成8位JPG影像

不多BB，直接上代码

private static void Resample()
        {
            Gdal.AllRegister();
            string fileName = @"D:\test.tif";
            string outFile = @"D:\result.jpg";
            Dataset ds = Gdal.Open(fileName, Access.GA_ReadOnly);
            OSGeo.GDAL.Driver drv = Gdal.GetDriverByName("JPEG");
            OSGeo.GDAL.Driver driver = Gdal.GetDriverByName("MEM");  // 内存类型，做数据中转
            int bandCount = ds.RasterCount;
            string strWkt = ds.GetProjectionRef();
            DataType dt = ds.GetRasterBand(1).DataType;
            double ratio = (double)ds.RasterXSize / ds.RasterYSize;
            double descX = 1024D / ds.RasterXSize;
            double descY = descX / ratio;
            int ySize = (int)(ds.RasterYSize * descY);
            Dataset newDs = driver.Create("mem", 1024, ySize, 3, DataType.GDT_Byte, null);
            Dataset newDsCopy = newDs;  // 复制一个Dataset,用来存储2%线性拉伸后的数据
            double[] geoTran = new double[6];
            ds.GetGeoTransform(geoTran);
            geoTran[1] = geoTran[1] / descX;
            geoTran[5] = geoTran[5] / descY;
            newDs.SetGeoTransform(geoTran);
            newDs.SetProjection(ds.GetProjection());
            newDsCopy.SetGeoTransform(geoTran);
            newDsCopy.SetProjection(ds.GetProjection());
            double[] scanLine = new double[1024 * ySize * Gdal.GetDataTypeSize(dt)];
            ds.GetRasterBand(1).ReadRaster(0, 0, ds.RasterXSize, ds.RasterYSize, scanLine, 1024, ySize, 0, 0);  //做数据重采样，转8位图并不需要重采样，我这里是因为项目需要
            double maxValue;
            double minValue;
            int hasValue;
            ds.GetRasterBand(1).GetMaximum(out maxValue, out hasValue);
            ds.GetRasterBand(1).GetMinimum(out minValue, out hasValue);
            for (int i = 0; i < scanLine.Length; i++)
            {
                if (scanLine[i] == 0) // 跳过无数据值，否则会出现负值
                {
                    continue;
                }
                scanLine[i] -= minValue;
                scanLine[i] /= maxValue;
                scanLine[i] *= 255;
            }

            newDs.GetRasterBand(1).WriteRaster(0, 0, 1024, ySize, scanLine, 1024, ySize, 0, 0);
            double[] scanLine2 = new double[1024 * ySize * Gdal.GetDataTypeSize(dt)];
            ds.GetRasterBand(2).ReadRaster(0, 0, ds.RasterXSize, ds.RasterYSize, scanLine2, 1024, ySize, 0, 0);
            ds.GetRasterBand(2).GetMaximum(out maxValue, out hasValue);
            ds.GetRasterBand(2).GetMinimum(out minValue, out hasValue);
            for (int i = 0; i < scanLine2.Length; i++)
            {
                if (scanLine2[i] == 0)
                {
                    continue;
                }
                scanLine2[i] -= minValue;
                scanLine2[i] /= maxValue;
                scanLine2[i] *= 255;
            }

            newDs.GetRasterBand(2).WriteRaster(0, 0, 1024, ySize, scanLine2, 1024, ySize, 0, 0);
            double[] scanLine3 = new double[1024 * ySize * Gdal.GetDataTypeSize(dt)];
            ds.GetRasterBand(3).ReadRaster(0, 0, ds.RasterXSize, ds.RasterYSize, scanLine3, 1024, ySize, 0, 0);
            ds.GetRasterBand(3).GetMaximum(out maxValue, out hasValue);
            ds.GetRasterBand(3).GetMinimum(out minValue, out hasValue);
            for (int i = 0; i < scanLine3.Length; i++)
            {
                if (scanLine3[i] == 0)
                {
                    continue;
                }
                scanLine3[i] -= minValue;
                scanLine3[i] /= maxValue;
                scanLine3[i] *= 255;
            }

            newDs.GetRasterBand(3).WriteRaster(0, 0, 1024, ySize, scanLine3, 1024, ySize, 0, 0);
            Dataset stretchSet = Liner2PercentStretch(newDs, newDsCopy, 0, 0, 1024, ySize, 1024, ySize);  // 进行2%线性拉伸
            drv.CreateCopy(outFile, stretchSet, 0, null, null, string.Empty);  // 将内存里数据保存到具体的文件里
            newDs.Dispose();
            stretchSet.Dispose();
            ds.Dispose();
            drv.Dispose();
            driver.Dispose();
            GC.Collect();
        }
/// <summary>
        /// 2%线性拉伸.
        /// </summary>
        /// <param name="srcSet">原始的Dataset数据.</param>
        /// <param name="desSet">拉伸后存储数据的Dataset.</param>
        /// <param name="pylY">ReadRaster函数的X偏移量.</param>
        /// <param name="pylX">ReadRaster函数的Y偏移量.</param>
        /// <param name="buffX">读取后数据的X尺寸.</param>
        /// <param name="buffY">读取后数据的Y尺寸.</param>
        /// <param name="xSize">读取前数据的X尺寸.</param>
        /// <param name="ySize">读取前数据的Y尺寸.</param>
        /// <returns>拉伸后的Dataset.</returns>
        private static Dataset Liner2PercentStretch(Dataset srcSet, Dataset desSet, int pylY, int pylX, int buffX, int buffY, int xSize, int ySize)
        {
            for (int i = 1; i < srcSet.RasterCount + 1; i++)
            {
                Band srcBand = srcSet.GetRasterBand(i);
                double max = 0;
                double min = 0;
                int hasValue = 0;
                srcBand.GetMaximum(out max, out hasValue);
                srcBand.GetMinimum(out min, out hasValue);
                Band desBand = desSet.GetRasterBand(i);
                byte[] buffer = new byte[buffX * buffY];
                int[] panHistogram = new int[256];
                double[] ratioHistogram = new double[256];
                srcBand.GetHistogram(-0.5, 255.5, 256, panHistogram, 1, 0, null, null);
                int totalCount = panHistogram.Sum();
                ratioHistogram[0] = (double)panHistogram[0] / totalCount;
                for (int j = 1; j <= 255; j++)
                {
                    ratioHistogram[j] = ratioHistogram[j - 1] + ((double)panHistogram[j] / totalCount);
                }
                srcBand.ReadRaster(pylY, pylX, xSize, ySize, buffer, buffX, buffY, 0, 0);
                int percent2Value = 0;  // 2%处灰度值
                int percent98Value = 0; // 98%处灰度值
                for (int m = 0; m <= 255; m++)
                {
                    if (ratioHistogram[m] <= 0.02)
                    {
                        percent2Value = m;
                    }
                    if (ratioHistogram[m] <= 0.98)
                    {
                        percent98Value = m;
                    }
                }
                for (int n = 0; n < buffer.Length; n++)
                {
                    if (buffer[n] <= percent2Value)
                    {
                        buffer[n] = 0;
                    }
                    else if (buffer[n] >= percent98Value)
                    {
                        buffer[n] = 255;
                    }
                    else
                    {
                        buffer[n] = (byte)(255 / (percent98Value - percent2Value) * (buffer[n] - percent2Value));
                    }
                }
                desBand.WriteRaster(0, 0, buffX, buffY, buffer, buffX, buffY, 0, 0);
            }
            return desSet;
        }

其中下面三行代码是转8位图的关键

scanLine2[i] -= minValue;
scanLine2[i] /= maxValue;
scanLine2[i] *= 255;

可以看到我在最后做了一个2%线性拉伸，这是因为只是单纯的转8位图，最后出来的结果很暗，效果非常不好，所以我做了一下2%线性拉伸。

不做线性拉伸结果：

拉伸后的结果，跟原图简直一毛一样：