GDAL 重采样？读取图像金字塔某分块？

Question

各位大虾，小弟现在用GDAL处理海量图像，学习了GDAL库后发现两个问题：
1、使用RasterIO不能进行采样设置，只能使用默认的nearest；
2、使用buildoverview可以设置，但是要读取其中的某一块却有困难。比如，现有一幅50000×50000的3波段图像，若缩放一倍后仍需分配25000×25000×3的空间存放数据集，超出内存范围会内存崩溃报错。我怎么读取其中的某一分块（比如前25000×1024部分)呢？
希望有良策的大虾们高手们前辈们不吝赐教，小弟感激不尽！

Answer 1

用RasterIO可以读取其中一部分内容啊。
我用过GDAL，但对它不熟，不少优秀的函数没用到。但重采样这个，即使RatserIO不支持，你也可以自己写啊。。我用GDAL只用到了读写，中间处理都是自己写……比如建立金字塔，重采样等。

Answer 2

想法不错，小弟需要支持多种格式的图像，弱弱问一句这自己建立金字塔和重采样是否和格式有关？还望不吝赐教哈，如果有什么实用的资料，可否分享一下？
对于GDAL我还有个困惑，帮助文档上貌似说颜色板不支持缩放，所以图像缩放颜色会有问题，孤独浮云老兄知否？

Answer 3

我处理的都是灰度图，或者在处理前将图片灰度化了。
我还没找到GDAL库调用调色板的函数。如果你能调用他的调色板读写函数，那问题就很简单了。
关于缩放，因为存在像素内插，对于灰度图其内插值仍然是灰度。对于彩色图，可以根据相邻像素3个分量梯度内插，会产生新的颜色。如果是二十四位等不带调色板的图，缩放就每任何问题。如果是带调色板的图，例如8位，其调色板只有256种颜色，新的颜色如果不是这其中的某种颜色，自然会出现问题。我觉得如果放大等等，可以把它最终存储为24位图片。这样就没任何问题了。如果强求8位的，1种方法是不用内插法，用最临近元（不产生新的颜色）。还有一种是存储前将24位真彩色转为8位彩色图（MSPAINT好像就是这么处理的）。但这样的效果反而不太好。如果GDAL提供了很好的处理函数，直接处理当然最好。实在不行，自己写也更加灵活。

Answer 4

贴3段代码，第一段是用GDAL读取数据，并转化为位图。
我处理的都是灰度图。所以对于8位的位图，GDAL读取数据信息虽然是调色板序号，但序号就是其灰度值。
（对于8位彩色图自然不是）其中str是打开路径，对于波段为1的图像，按照灰度图读取。对于波段》=3的图像，读取其中3个波段进行灰度化。bool CPicMatchDoc::CreatImage(int LorR, CString str)
{

GDALDataset *poDataset; //数据集对象指针
GDALAllRegister();  //注册驱动 USES_CONVERSION; // 在前面加上这句
char *pz =T2A(str); poDataset = (GDALDataset *) GDALOpen(pz, GA_ReadOnly); if (poDataset == NULL) /*检查是否正常打开文件*/
{
if(LorR==0)
AfxMessageBox(_T("打开左片影像文件失败！"),0,0);
else
   AfxMessageBox(_T("打开右片影像文件失败！"),0,0);
return 0;
} BeginWaitCursor(); //输入图像的信息
int nBands = poDataset->GetRasterCount(); GDALRasterBand ** poBand; //各波段信息指针 poBand = new GDALRasterBand *[nBands]; if (poBand == NULL)
{
AfxMessageBox(_T("创建各波段数据集失败！"), MB_ICONWARNING,0);
return 0;
} for (int i=0; i<nBands; i++)
{
poBand[i] = poDataset->GetRasterBand(i+1);
if (poBand[i] == NULL)
{
CString str;
str.Format(_T("创建第%d波段数据集失败！"),i+1);
AfxMessageBox(str, MB_ICONWARNING);
return 0;
}
} CImageInfoDlg dlg;
dlg.Totalband=nBands;
if(LorR==0)
dlg.dlgType=1;
else
dlg.dlgType=2;
dlg.DoModal(); int nXsize,nYsize; //图像高宽
nXsize = poBand[0]->GetXSize();
nYsize = poBand[0]->GetYSize();
int maxSize=max(nXsize,nYsize); if (LorR==0) //左片
{ CPicTreeView *pView=(CPicTreeView*)(((CMainFrame*)AfxGetMainWnd())-> m_wndSplitter.GetPane(0,0));
CTreeCtrl& CtlTree=(CTreeCtrl&)pView->GetTreeCtrl (); //根节点
HTREEITEM htemroot=CtlTree.GetRootItem();
//左片
HTREEITEM htemleft=CtlTree.GetChildItem(htemroot);
HTREEITEM htemdele=CtlTree.GetChildItem(htemleft);
// htemdele=CtlTree.GetNextSiblingItem(htemdele); while (htemdele!=NULL)
{
HTREEITEM htemNext=CtlTree.GetNextItem(htemdele,TVGN_NEXT);
CtlTree.DeleteItem(htemdele);
htemdele=htemNext;
}

if(600<=maxSize&&maxSize<=1200)
nImgL=2;
else
{
if (maxSize<=2500)
nImgL=3;
else
{
if (maxSize<=5000)
nImgL=4;
else
nImgL=5;
} } pBitmapInfoL=new BITMAPINFO *[nImgL];
pDataL=new LPBYTE[nImgL]; pBitmapInfoL[0] = (BITMAPINFO*)new char[sizeof(BITMAPINFO) + 
sizeof(RGBQUAD) * (256)];  pBitmapInfoL[0]->bmiHeader.biClrUsed = 0;
pBitmapInfoL[0]->bmiHeader.biBitCount = 8;
pBitmapInfoL[0]->bmiHeader.biClrImportant = 0;
pBitmapInfoL[0]->bmiHeader.biCompression = BI_RGB;
pBitmapInfoL[0]->bmiHeader.biWidth = nXsize;
pBitmapInfoL[0]->bmiHeader.biHeight = nYsize;
pBitmapInfoL[0]->bmiHeader.biPlanes = 1;
pBitmapInfoL[0]->bmiHeader.biSize = sizeof(BITMAPINFOHEADER);//sizeof(BITMAPINFO);
pBitmapInfoL[0]->bmiHeader.biSizeImage = (nXsize*8+31)/32*4*nYsize; 
pBitmapInfoL[0]->bmiHeader.biXPelsPerMeter = 0;
pBitmapInfoL[0]->bmiHeader.biYPelsPerMeter = 0; //为颜色表赋值
for (int i=0; i<256; i++)
{
pBitmapInfoL[0]->bmiColors[i].rgbBlue = i;
pBitmapInfoL[0]->bmiColors[i].rgbGreen = i;
pBitmapInfoL[0]->bmiColors[i].rgbRed = i;
pBitmapInfoL[0]->bmiColors[i].rgbReserved = 0;
}
}
else
{ CPicTreeView *pView=(CPicTreeView*)(((CMainFrame*)AfxGetMainWnd())-> m_wndSplitter.GetPane(0,0)); CTreeCtrl& CtlTree=(CTreeCtrl&)pView->GetTreeCtrl (); //根节点
HTREEITEM htemroot=CtlTree.GetRootItem();
//左片
HTREEITEM htemleft=CtlTree.GetChildItem(htemroot);
//右片
htemleft=CtlTree.GetNextSiblingItem(htemleft);
HTREEITEM htemdele=CtlTree.GetChildItem(htemleft);
// htemdele=CtlTree.GetNextSiblingItem(htemdele); while (htemdele!=NULL)
{
HTREEITEM htemNext=CtlTree.GetNextItem(htemdele,TVGN_NEXT);
CtlTree.DeleteItem(htemdele);
htemdele=htemNext;
}

if(600<=maxSize&&maxSize<=1200)
nImgR=2;
else
{
if (maxSize<=2500)
nImgR=3;
else
{
if (maxSize<=5000)
nImgR=4;
else
nImgR=5;
} } pBitmapInfoR=new BITMAPINFO *[nImgR];
pDataR=new LPBYTE[nImgR]; pBitmapInfoR[0] = (BITMAPINFO*)new char[sizeof(BITMAPINFO) + 
sizeof(RGBQUAD) * (256)];  pBitmapInfoR[0]->bmiHeader.biClrUsed = 0;
pBitmapInfoR[0]->bmiHeader.biBitCount = 8;
pBitmapInfoR[0]->bmiHeader.biClrImportant = 0;
pBitmapInfoR[0]->bmiHeader.biCompression = BI_RGB;
pBitmapInfoR[0]->bmiHeader.biWidth = nXsize;
pBitmapInfoR[0]->bmiHeader.biHeight = nYsize;
pBitmapInfoR[0]->bmiHeader.biPlanes = 1;
pBitmapInfoR[0]->bmiHeader.biSize = sizeof(BITMAPINFOHEADER);//sizeof(BITMAPINFO);
pBitmapInfoR[0]->bmiHeader.biSizeImage = (nXsize*8+31)/32*4*nYsize; 
pBitmapInfoR[0]->bmiHeader.biXPelsPerMeter = 0;
pBitmapInfoR[0]->bmiHeader.biYPelsPerMeter = 0; //为颜色表赋值
for (int i=0; i<256; i++)
{
pBitmapInfoR[0]->bmiColors[i].rgbBlue = i;
pBitmapInfoR[0]->bmiColors[i].rgbGreen = i;
pBitmapInfoR[0]->bmiColors[i].rgbRed = i;
pBitmapInfoR[0]->bmiColors[i].rgbReserved = 0;
}
}
BYTE *poBandBlock[3]; //存储像素数据 if (dlg.nstyle == 1) //单一波段打开灰度影像
{
poBandBlock[0]= (BYTE*)CPLMalloc(sizeof(BYTE)*nXsize);
//建立数据区
LONG LineBytes = (nXsize*8+31)/32*4;
if (LorR==0)
pDataL[0]= (LPBYTE)new char[LineBytes*nYsize];
else
pDataR[0]= (LPBYTE)new char[LineBytes*nYsize]; for (int i=0; i<nYsize; i++)
{

//一次性读取
poBand[dlg.Gray]->RasterIO(GF_Read, 0, i, nXsize, 
1, poBandBlock[0],nXsize, 1,GDT_Byte, 0, 0);

if (LorR==0) //left
{
//为各像素赋值
for (int j=0; j<nXsize; j++)
pDataL[0][(nYsize-i-1)*LineBytes + j] = poBandBlock[0][j]; }
else
{
//为各像素赋值
for (int j=0; j<nXsize; j++)
pDataR[0][(nYsize-i-1)*LineBytes + j] = poBandBlock[0][j];
}
} CPLFree(poBandBlock[0]); //释放GDAL读取数据

} if (dlg.nstyle == 0) //多波段灰度化打开影像
{
//一次性读取
poBandBlock[0]= (BYTE*)CPLMalloc(sizeof(BYTE)*nXsize);
poBandBlock[1]= (BYTE*)CPLMalloc(sizeof(BYTE)*nXsize);
poBandBlock[2]= (BYTE*)CPLMalloc(sizeof(BYTE)*nXsize); //建立数据区
LONG LineBytes = (nXsize*8+31)/32*4;
if (LorR==0)
pDataL[0]= (LPBYTE)new char[LineBytes*nYsize];
else
pDataR[0]= (LPBYTE)new char[LineBytes*nYsize]; for (int i=0; i<nYsize; i++)
{
poBand[dlg.Red]->RasterIO(GF_Read, 0, i, nXsize, 
1, poBandBlock[0], nXsize, 1, GDT_Byte, 0, 0);
poBand[dlg.Green]->RasterIO(GF_Read, 0, i, nXsize, 
1, poBandBlock[1], nXsize, 1, GDT_Byte, 0, 0);
poBand[dlg.Blue]->RasterIO(GF_Read, 0, i, nXsize, 
1, poBandBlock[2], nXsize, 1, GDT_Byte, 0, 0); if (LorR==0) //left
{
//为各像素赋值
for (int j=0; j<nXsize; j++)
pDataL[0][(nYsize-i-1)*LineBytes + j] = (unsigned char)(0.299*poBandBlock[0][j]
+0.587*poBandBlock[1][j] + 0.114*poBandBlock[2][j]+0.5);
}
else
{
for (int j=0; j<nXsize; j++)
//为各像素赋值
pDataR[0][(nYsize-i-1)*LineBytes + j] = (unsigned char)(0.299*poBandBlock[0][j]
+0.587*poBandBlock[1][j] + 0.114*poBandBlock[2][j]+0.5);
}
} CPLFree(poBandBlock[0]); //释放GDAL读取数据
CPLFree(poBandBlock[1]); //释放GDAL读取数据
CPLFree(poBandBlock[2]); //释放GDAL读取数据 }
EndWaitCursor(); return 1;
}

Answer 5

第二段代码是创建金字塔影像。事实上我是定义了一个位图文件头数组，数据区数组来存放金字塔影像信息。（当然由于对于不同分辨率图像，金字塔分层不同。所以数组是动态的，用指针实现）金字塔影像创建一般为先低通滤波，然后重采样。
低通滤波一般选用奇数大小窗口的模板。
我这段代码其实是为了计算方便，选用的是2X2的平均低通模板。由于用双线性内插法重采样时，若缩小的倍数是整数倍（我每层是按照1：2缩小），则双线性内插此时就是最临近元了（因为重采样像素值是整数）。把这2步结合起来考虑，也就是可以直接利用上一层4个像素值用2X2平均就得到下一层影像对应像素值。
这样比一般创建步骤少一个循环。这只是种特殊情况。一般还是要先低通滤波卷积，利用低通后图像的像素信息重采样。// 创建金字塔影像
bool CPicMatchDoc::CreatPyramid(int LorR)
{
CPicTreeView *pView=(CPicTreeView*)(((CMainFrame*)AfxGetMainWnd())-> m_wndSplitter.GetPane(0,0)); int i,j,k; //循环变量
if (LorR==0) //左片金字塔影像
{
for(i=1;i<nImgL;i++)
{
//影像高宽
int nXsize=pBitmapInfoL[i-1]->bmiHeader.biWidth/2;
int nYsize=pBitmapInfoL[i-1]->bmiHeader.biHeight/2; pBitmapInfoL[i] = (BITMAPINFO*)new char[sizeof(BITMAPINFO) + 
sizeof(RGBQUAD) * (256)];  pBitmapInfoL[i]->bmiHeader.biClrUsed = 0;
pBitmapInfoL[i]->bmiHeader.biBitCount = 8;
pBitmapInfoL[i]->bmiHeader.biClrImportant = 0;
pBitmapInfoL[i]->bmiHeader.biCompression = BI_RGB;
pBitmapInfoL[i]->bmiHeader.biWidth = nXsize;
pBitmapInfoL[i]->bmiHeader.biHeight = nYsize;
pBitmapInfoL[i]->bmiHeader.biPlanes = 1;
pBitmapInfoL[i]->bmiHeader.biSize = sizeof(BITMAPINFOHEADER);//sizeof(BITMAPINFO);
pBitmapInfoL[i]->bmiHeader.biSizeImage = (nXsize*8+31)/32*4*nYsize; 
pBitmapInfoL[i]->bmiHeader.biXPelsPerMeter = 0;
pBitmapInfoL[i]->bmiHeader.biYPelsPerMeter = 0; //为颜色表赋值
for (j=0; j<256; j++)
{
pBitmapInfoL[i]->bmiColors[j].rgbBlue = j;
pBitmapInfoL[i]->bmiColors[j].rgbGreen = j;
pBitmapInfoL[i]->bmiColors[j].rgbRed = j;
pBitmapInfoL[i]->bmiColors[j].rgbReserved = 0;
} //建立数据区
LONG LineBytes = (nXsize*8+31)/32*4;
pDataL[i]= (LPBYTE)new char[LineBytes*nYsize]; //上一层图片数据
int lastx,lasty;
lastx=pBitmapInfoL[i-1]->bmiHeader.biWidth;
lasty=pBitmapInfoL[i-1]->bmiHeader.biHeight;
LONG lastlinebyte =(lastx*8+31)/32*4; for(j=0;j<nYsize;j++)
{
for (k=0;k<nXsize;k++)
pDataL[i][(nYsize-j-1)*LineBytes + k]=(unsigned char)(pDataL[i-1][(lasty-2*j-1)*lastlinebyte+2*k]*0.25+
pDataL[i-1][(lasty-2*j-2)*lastlinebyte+2*k-1]*0.25+
pDataL[i-1][(lasty-2*j-2)*lastlinebyte+2*k]*0.25+
pDataL[i-1][(lasty-2*j-1)*lastlinebyte+2*k-1]*0.25+0.5);
} } for (i=0;i<nImgL+1;i++)
{
CTreeCtrl& CtlTree=(CTreeCtrl&)pView->GetTreeCtrl (); TV_INSERTSTRUCT tvThree;//树叶 HTREEITEM htemroot=CtlTree.GetRootItem();
HTREEITEM htemleft=CtlTree.GetChildItem(htemroot);
//CString stemp=CtlTree.GetItemText(htemleft);
CString strTV; if (i==0)
strTV.Format(_T("特征点列表"));
else
strTV.Format(_T("第 %d 层"),i);
tvThree.hParent=htemleft;
tvThree.item.pszText=(LPTSTR)(LPCTSTR)strTV;
tvThree.item.mask=TVIF_TEXT|TVIF_IMAGE|TVIF_SELECTEDIMAGE;
tvThree.item.iImage=1;
tvThree.item.iSelectedImage=2;
CtlTree.InsertItem (&tvThree);

}
}
else
{
for(i=1;i<nImgR;i++)
{
//影像高宽
int nXsize=pBitmapInfoR[i-1]->bmiHeader.biWidth/2;
int nYsize=pBitmapInfoR[i-1]->bmiHeader.biHeight/2; pBitmapInfoR[i] = (BITMAPINFO*)new char[sizeof(BITMAPINFO) + 
sizeof(RGBQUAD) * (256)];  pBitmapInfoR[i]->bmiHeader.biClrUsed = 0;
pBitmapInfoR[i]->bmiHeader.biBitCount = 8;
pBitmapInfoR[i]->bmiHeader.biClrImportant = 0;
pBitmapInfoR[i]->bmiHeader.biCompression = BI_RGB;
pBitmapInfoR[i]->bmiHeader.biWidth = nXsize;
pBitmapInfoR[i]->bmiHeader.biHeight = nYsize;
pBitmapInfoR[i]->bmiHeader.biPlanes = 1;
pBitmapInfoR[i]->bmiHeader.biSize = sizeof(BITMAPINFOHEADER);//sizeof(BITMAPINFO);
pBitmapInfoR[i]->bmiHeader.biSizeImage = (nXsize*8+31)/32*4*nYsize; 
pBitmapInfoR[i]->bmiHeader.biXPelsPerMeter = 0;
pBitmapInfoR[i]->bmiHeader.biYPelsPerMeter = 0; //为颜色表赋值
for (j=0; j<256; j++)
{
pBitmapInfoR[i]->bmiColors[j].rgbBlue = j;
pBitmapInfoR[i]->bmiColors[j].rgbGreen = j;
pBitmapInfoR[i]->bmiColors[j].rgbRed = j;
pBitmapInfoR[i]->bmiColors[j].rgbReserved = 0;
} //建立数据区
LONG LineBytes = (nXsize*8+31)/32*4;
pDataR[i]= (LPBYTE)new char[LineBytes*nYsize]; //上一层图片数据
int lastx,lasty;
lastx=pBitmapInfoR[i-1]->bmiHeader.biWidth;
lasty=pBitmapInfoR[i-1]->bmiHeader.biHeight;
LONG lastlinebyte =(lastx*8+31)/32*4; for(j=0;j<nYsize;j++)
{
for (k=0;k<nXsize;k++)
pDataR[i][(nYsize-j-1)*LineBytes + k]=(unsigned char)(pDataR[i-1][(lasty-2*j-1)*lastlinebyte+2*k]*0.25+
pDataR[i-1][(lasty-2*j-2)*lastlinebyte+2*k-1]*0.25+
pDataR[i-1][(lasty-2*j-2)*lastlinebyte+2*k]*0.25+
pDataR[i-1][(lasty-2*j-1)*lastlinebyte+2*k-1]*0.25+0.5);
} }
for (i=0;i<nImgR+1;i++)
{
CTreeCtrl& CtlTree=(CTreeCtrl&)pView->GetTreeCtrl (); TV_INSERTSTRUCT tvThree;//树叶 HTREEITEM htemroot=CtlTree.GetRootItem();
HTREEITEM htemleft=CtlTree.GetChildItem(htemroot);
HTREEITEM htemRight=CtlTree.GetNextSiblingItem(htemleft); CString strTV; if (i==0)
strTV.Format(_T("特征点列表"));
else
strTV.Format(_T("第 %d 层"),i);
tvThree.hParent=htemRight;
tvThree.item.pszText=(LPTSTR)(LPCTSTR)strTV;
tvThree.item.mask=TVIF_TEXT|TVIF_IMAGE|TVIF_SELECTEDIMAGE;
tvThree.item.iImage=1;
tvThree.item.iSelectedImage=2;
CtlTree.InsertItem (&tvThree); } }

return 1;
}

Answer 6

非常感谢孤独浮云兄如此热心回答，但是我发现代码都是位图的，我需要支持BMP JPG PNG IMG TIF GIF等等各种格式，一个个解析其格式不太可能，所以只能依靠开源库了。利用GDAL可以读取其波段数据，重采样和建立金子塔若和图像格式无关的话可以自己写，若有关的话也变的不现实了。
我再找找资料
再次表示感谢，结帖给分，有关于GDAL深入研究的资料的话希望能发我一份哈，邮箱[email protected]。

Answer 7

并没有一个一个解析其图片格式啊。
原理是将所有格式图片数据用GDAL读取，构造一副位图，对位图处理。存储的时候也可以选多种方式。这样调色板之类的问题都容易解决了。
你可以看看这个
http://download.csdn.net/source/2380538
我的这个里面有建立金字塔以及支持你所需要的格式的读、写。