遙感影像16位轉8位
阿新 • • 發佈:2019-01-06
現在常用衛星影像基本上都是16位影像,如GF,ZY3,Landsat8,WV等,有時我們需要將16位影像降到8位影像,這樣不僅減小了資料量,也便於後期處理。
通常的軟體在處理降位時會存在一些問題,如曝光,出現空值等。因此,我自己根據常用的降位方法,如最簡單的線性拉伸,分段拉伸以及對數變換和指數變換,結合常用影像的特點,使用百分比截斷和指數(冪為0.7)變換將影像從16位降到8位,在抑制高光的同時保證了影像的對比度。
int imageprocessing::stretch_percent_16to8(const char *inFilename, const char *dstFilename) { GDALAllRegister(); //為了支援中文路徑 CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO"); int src_height = 0; int src_width = 0; GDALDataset *poIn = (GDALDataset *)GDALOpen(inFilename,GA_ReadOnly); //開啟影像 //獲取影像大小 src_width = poIn ->GetRasterXSize(); src_height = poIn ->GetRasterYSize(); //獲取影像波段數 int InBands = poIn ->GetRasterCount(); //獲取影像格式 GDALDataType eDataType = poIn -> GetRasterBand(1) -> GetRasterDataType(); //定義儲存影像的空間參考陣列 double adfInGeoTransform[6] = {0}; const char *pszWKT = NULL; //獲取影像空間參考 poIn ->GetGeoTransform(adfInGeoTransform); pszWKT = poIn ->GetProjectionRef(); //建立檔案 GDALDriver *poDriver = (GDALDriver *)GDALGetDriverByName("GTiff"); GDALDataset *poOutputDS = poDriver -> Create(dstFilename,src_width,src_height,InBands,GDT_Byte,NULL); //設定拉伸後圖像的空間參考以及地理座標 poOutputDS -> SetGeoTransform(adfInGeoTransform); poOutputDS -> SetProjection(pszWKT); //讀取影像 cout<<"16位影像降到8位影像處理..."<<endl; for(int iBand = 0; iBand < InBands; iBand++ ) { cout<<"正在處理第 "<<iBand + 1<<" 波段"<<endl; //讀取影像 uint16_t *srcData = (uint16_t *)malloc(sizeof(uint16_t) *src_width * src_height *1); memset(srcData,0,sizeof(uint16_t ) * 1 *src_width * src_height); int src_max = 0, src_min = 65500; //讀取多光譜影像到快取 poIn ->GetRasterBand( iBand + 1) -> RasterIO( GF_Read, 0, 0, src_width, src_height , srcData + 0 * src_width * src_height,src_width,src_height, GDT_UInt16, 0, 0 ); //} //統計最大值 for (int src_row = 0; src_row < src_height; src_row ++) { for (int src_col = 0; src_col < src_width; src_col++) { uint16_t src_temVal = *(srcData + src_row * src_width + src_col); if (src_temVal > src_max) src_max = src_temVal; if(src_temVal < src_min ) src_min = src_temVal; } } double *numb_pix = (double *)malloc(sizeof(double)*(src_max+1)); //存畫素值直方圖,即每個畫素值的個數 memset(numb_pix,0,sizeof(double) * (src_max+1)); // ------- 統計畫素值直方圖 ------------ // for (int src_row = 0; src_row < src_height; src_row ++) { for (int src_col = 0; src_col < src_width; src_col++) { uint16_t src_temVal = *(srcData + src_row * src_width + src_col); *(numb_pix + src_temVal) += 1; } } double *frequency_val = (double *)malloc(sizeof(double)*(src_max+1)); //畫素值出現的頻率 memset(frequency_val,0.0,sizeof(double)*(src_max+1)); for (int val_i = 0; val_i <= src_max; val_i++) { *(frequency_val + val_i) = *(numb_pix + val_i) / double(src_width * src_height); } double *accumlt_frequency_val = (double*)malloc(sizeof(double)*(src_max+1)); //畫素出現的累計頻率 memset(accumlt_frequency_val, 0.0,sizeof(double)*(src_max+1)); for (int val_i = 0; val_i <= src_max; val_i ++) { for (int val_j = 0; val_j < val_i; val_j ++ ) { *(accumlt_frequency_val + val_i) += *(frequency_val + val_j); } } //統計畫素兩端截斷值 int minVal = 0, maxVal = 0; for (int val_i = 1; val_i < src_max; val_i++) { double acc_fre_temVal0 = *(frequency_val + 0); double acc_fre_temVal = *(accumlt_frequency_val + val_i); if((acc_fre_temVal - acc_fre_temVal0) > 0.0015 ) { minVal = val_i; break; } } for (int val_i = src_max-1; val_i > 0; val_i--) { double acc_fre_temVal0 = *(accumlt_frequency_val + src_max); double acc_fre_temVal = *(accumlt_frequency_val + val_i); if(acc_fre_temVal < (acc_fre_temVal0 - 0.00012) ) { maxVal = val_i; break; } } for (int src_row = 0; src_row < src_height; src_row ++) { uint8_t *dstData = (uint8_t*)malloc(sizeof(uint8_t)*src_width); memset(dstData, 0, sizeof(uint8_t)*src_width); for (int src_col = 0; src_col < src_width; src_col++) { uint16_t src_temVal = *(srcData + src_row * src_width + src_col); double stre_temVal = (src_temVal - minVal) / double(maxVal - minVal) ; if(src_temVal < minVal) { *(dstData + src_col) = (src_temVal) *(20.0/double(minVal)) ; } else if(src_temVal > maxVal) { stre_temVal = (src_temVal - src_min) / double(src_max - src_min); *(dstData + src_col) = 254; } else *(dstData + src_col) = pow(stre_temVal,0.7) * 250; } poOutputDS->GetRasterBand(iBand + 1)->RasterIO(GF_Write, 0,src_row,src_width,1,dstData,src_width,1,GDT_Byte,0,0); free(dstData); } free(numb_pix); free(frequency_val); free(accumlt_frequency_val); free(srcData); } GDALClose(poIn); GDALClose(poOutputDS); return 0; }