GDAL Qt 开发

       最近把投影转换终于弄清楚了，在gdal中，只是支持部分投影的转换，这些投影是比较常见的几种类型，如果需要其他gdal不直接提供的，可以自己根据proj.4编写。下一步的任务是用gdal读取文件数据，用qt提供的界面进行显示。需要解决的关键问题是大数据（如n G的数据），怎么在文件中进行显示的问题。

      解决的问题是初步想法：

     一：通过在qt界面上获取显示的范围，如果显示范围小于图像的范围，则只读取界面范围的图像。具体1：用gdal读取图像像素信息，2：根据像素信息构造qimage，3：将qimage显示在qpixmap中，4：研究金字塔图像）。

     二：建一个qtGUI工程，准备一个做几何纠正的框架。具体

         1：在工程中，首先需要打开文件，显示文件信息对话框，用树状图显示。

         2：在对话框显示图像。首先显示金字塔影像全图，先在geometry类中把qimage构造出来，然后在对应的类中进行显示。需要使用多视图来对影像进行显示，以便后续的选点工作，建立的新的widget（必须基于Qmainwindows来建，不然在后面无法将qscroolArea设置为主widget）。

         3：在缩略图中显示红色的范围框，首先弄清楚金字塔影像产生的规则，即可以得到红线框所在区域影像的原始范围，通过这个范围，用影像分块读取法，读取影像像素，继而构造qimage：根据金字塔的层数，可以得到现在的像素在原图像上的位置。如：第n层，则对应原图像的像素坐标为2的n次方（需要判断有没有溢出）。

遥感影像过大，而且数据类型多样，因此显示是个问题，现有解决方案：

       利用金字塔加速，这种方案应用还算比较用的，如果不用这种方案，在缩放的时候就会速度很慢，但是会比较占硬盘，特别是1个多G的数据的时候。一般是2倍缩放，但是这种方案不知道是取平均还是直接取4个像素中的某个，按理说第二种是可行的，如果仅仅是显示的话，因为具体选点的时候，如果仅在2倍大小选点的话，还是有0.5个像素的选点误差。在实际处理的时候，我还是倾向3倍。

      利用帧缓存技术，这种方案使得拖动较为流畅，但是还是有缝的，即拖动的时候会有黑块出现，这样处理的方案是再拖动的时候，计算出需要显示的东东，存储在某个对象中，显示时交换即可。改进方案，因为屏幕一般不过1280×1024，如果我以显示的部分为中心，读取9倍大小的影像块(内存要的也不过10M)，这样在拖动的时候，你不管怎么拖动都在我的影像范围内，这样拖动就会显得无缝，在拖动时，还是需要一个缓存，来存储要需用的区域，拖动完的时候进行交换。

     多种数据格式的话，可以在生成金字塔影像的时候，对影像的各个波段进行计算，获取其灰度分布直方图，然后显示的时候进行实时计算，将原始格式转换为8位位图。

         4：当在缩略图中显示确定显示范围时，将此范围换算为原图像的范围，并把该范围交给下一个处理大图的ui，缩略图所在ui类中，并没有获得主界面ui的gdaldataset指针，但是在下一个显示原图的范围的过程中，我们需要获得原图像数据信息，所以，在构造下一个显示界面时，我们必须把相应的数据指针传进去，还要把相应的overview层数传进去。：在缩略图中显示选框大小时，应该把金字塔层数带进去计算，如果选框所表示的原始图像过大，将耗损内存严重，这一步放在优化阶段进行。

         5：首先根据overview层数，起始点与终止点数，计算出原图像的起始点，终止点。通过这两个点的位置，读取原数据图像。图像读取完毕，问题是显示范围太大，在主视图中调整。设计在主界面中设计每次显示原始图像的500*500.不用scroolarea。这样界面更清新。通过金字塔层数，算出500*500在缩略图中的位置，从而确定线框的大小。实验证明，这种限制死了原图像显示的方法不好，因为在图像很大时，或者金字塔层数较多时，这样使得在缩略图的显示界面上，矩形框变得非常小。既不能把矩形框限制死，也不能把缩略图限制死，而是在中间寻找一个平衡，即使得缩略图中的矩形框大小适中，也要让图形显示大小适中，实验：只记录缩略图中的线框起始点，对于原图像界面只显示固定大小。而缩略图，总是显示金字塔的最上层。

        6：初步实现界面显示工作，但是在检查过程中发现：当显示的图像精度太高，而只显示缩略图框左上角500个像素点的位置根本不够用，所以还是要想办法把50*50的缩略图中的内容全部显示出来，解决方案：用graphicsview来显示每个block块的图像。研究graphicsview的实现方式。现在确定用layout可以解决这个问题。把所有的图像label全部排列在gridlayout中，根绝图像block的大小，确定每一个label的大小。当试图变换时，显示新区域的label，同时把以前的image删除掉。从而实现不占用太多内存的目的。

 实验数据表明：每个界面要显示8*8个qimge，但是label.setpixel（）非常耗时，估计也非常耗费内存。不然不会那么慢，现在只有通过 graphicview来显示这大量的qimage了。先实验看效果。实验表明，用分块读取影像，然后放在qgraphics中，显示效率是相当高的，而且，如果是第二个显示界面是subclassing和reimplementation 那个graphicsview类的话，它还能支持海量图像的自动管理，省去了每次缩放都调用金字塔影像的麻烦。

       7：。感觉每次显示的时候，两个显示框的位置好像不大对劲，先检查一下，看有没有问题，在计算金字塔的scalefactor时出问题了。把范围算小了，这样一来，当移动缩略图中的范围框时，对应区域的刷新就相当慢了，解决这个问题，用多线程计算来解决，qthread实验，由于gdal中的rater不能支持qt中的并行（？）。所以不用并行了，寻找其他方案，思路为减少缩略图中的框的大小，这个大小应该根据影像的范围的动态的确定。下一步实现这个。先实验，看看主屏在显示多大的时不卡。（图像横纵为1500时比较合适）。

       8：对于warp1中的缩略图显示，采用的是继承mainwindow的形式，对比可知，warp2显示效果较好，对warp1进行修改，让它也继承QGraphicsView来进行显示，这样对于十字叉过小的问题就可以解决了，实验：这样rect和image不在一个图层，解决方法：设置z值，失败，暂时解决方案，在realeasemouse中移动方框，可以保证两个item的相对位置不变。在移动方框的时候，setpos出现了bug，现在认定是qt出现的一个bug，因为在方框在左上角的时候，操作时没有问题的。而且

 QPointF px;
 px=pRect+QPointF(dx,dy);
 items->setPos(px);
 px=items->pos()+QPointF(dx,dy);
 items->setPos(px);与下面这句

 QPointF px;
 px=items->pos()+QPointF(dx,dy);
 items->setPos(px);当在第一次移动的时候，竟然效果一样，这完全不能解释，当第二次拖动的时候，两者就不一样了（本来第一次就应该不一样呀，还请高人指教）。

初步打算就把方框初始位置放在左上角算了。本来这也不影响最终结果的。这个问题纠结了一上午加上大半个下午，真不值得，以后干嘛还是得先参考帮助文档，按照帮助文档上面标准的做法来做，因为所有的东西不可能做到尽善尽美，你不可能做到，人家也不可能，弹性可能没有你想象的那么大，所以还是按照标准来，不然指不定什么时候就出问题了，在这种问题上耗时，非常不划算。

       8：subclassing and reimplementation：重载qgraphicsview的功能， 在graphicsview中有两套坐标，一套是view的坐标。一套是scene中的坐标。由于在放置scene中，我们已经指定了单个scene的位置。所以，每次获取坐标时候需要把它用maptoscene转换到scene下。这个坐标加上该图本来的左上角图像坐标，就可以得出任意一点，在原图像上的坐标了。

     9：研究envi的几何纠正过程。对RST（resampling，scaling and translation）进行研究。对于几何纠正的理解，下图表示：

 

 

 

 

 

 根据控制点，图像上面的点转换到地理坐标范围，这个过程是translate，通过地理坐标范围，确定新的图像范围，并对空白处进行插值，这是resampling，根据要求的像素，对图像进行缩放，这是scaling。

 10：完成程序中base的界面（和warp中一样），完成选点界面及其相关控制的设计。界面之间的通信在qt中是怎么完成的，推测信号槽机制应该可以完成。设置全局变量四个，分别记录base和warp中的像素点坐标（pointcentral加上左上角坐标即可得出）。全局变量的使用比较麻烦，而且不符合面向对象的规则，最好不用，在这里我们使用qt的信号槽机制解决这一问题。把坐标位置一层一层传到住界面类中。

 11：设计控制点选取界面，一个表格table，三个putton；实现选点后自动显示

12：研究gdal图像translate工具的实现，在研究的过程中，发现以前的关于投影变换的东西需要重新改进。现在暂时不打乱计划，先做几何纠正。重点研究gcp部分的实现情况。几何纠正流程（参考gdal中的变换）
a：获取基准影像的对应点的地理坐标，获取方式为仿射变换，结合纠正影像的影像坐标位置，就可以得到纠正影像的GCPs位置了，写入gdal_gcp链表中。
      Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
      Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
b：用gdal打开纠正影像GDALOpen，获取其波段数nBandCount；获取影像的长和宽nRasterXSize、nRasterYSize；
c：获取输出格式（与纠正影像格式应该保持一致）driver:hDriver；
d：申明一个virtual dataset：poVDS；
e：计算仿射变换参数，并设置为poVDS的仿射变换参数（错误，在这里不能设置放射变换参数，而是应该用控制点模式）。获取地理参考坐标系，设置poVDS控制点。
f： 设置poVDS元数据。
g：获取波段数据类型，申明vrt raster band变量，逐波段处理纠正影像数据，poVDS中增加一个波段-将源波段添加到vrt波段中去。
h：creatcopy将数据创建目标文件。释放过程中申明的变量。

//至此，创建了虚拟的dataset，与原始文件相比，文件中增加了新的控制点文件，其他的不变。有两个思路去做几何纠正，第一个：直接创建中间文件。把中间文件带入到
   gdalwarp中去，形成新的最终文件。第二：利用虚拟的dataset，通过对它进行warp变换，从而达到纠正的目的。实验表明，第二种方法不易实现，所以采用第一种方法。

i：将中间文件看做源文件，创建目标文件，从源文件向目标文件转换。删除源文件。

13：考虑一个新问题，当base有投影，warp也有投影的时候，是采用base下得投影还是采用warp下得投影呢。应该是base下得投影，因为如果采用warp，在选点计算实际地理坐标时，不能利用在base中采集到的像素点，加上warp中的仿射变换参数来计算。只能用base中的像素点和base的仿射变换参数，计算地理坐标，这个地理坐标同时也就是warp中的像素点对应的地理坐标，这样来说，warp中的点都应该被变换到base投影类型下，所以在新文件设置投影的时候，应该以base中的投影为准，而波段信息，当然是采用warp中的。

14：调试程序，问题一：base2中方框移动问题，解决，    问题二：参数传递问题，解决，   问题三：单位问题，有的仿射变换单位为经纬度，不是meter，怎么办，直接按照degree进行计算即可。都是degree坐标，统一即可，解决。     问题四：读取格式时，没有预先读取波段格式。而是只用了byte，先实验转换部分，这个以后完善。问题五：poVDS无法读取数据，看gdal中提供tranlate中的实现，问题是忘记给povds加上波段，解决。      问题五：转换之后的图像，是以gcp的形式描述的，首先要为其加上投影信息，然后再研究如何把gcp形式转化到左上角点形式。实验：在vrt的控制点中增加投影信息，得到的投影信息是是控制点的，现在研究gdalwarp中的变换方式，把控制点附带的投影信息，运用到全局上，这一问题回溯到几何纠正的流程上，解决。    问题六：拖动之后，放大界面的框不见了，解决。  问题七:当图像范围较小时，缩略显示框有问题，解决。   

15：到目前为止，程序基本上可以完成所要求的任务。下面进行程序优化。     一：主界面改造，解决；    二：投影信息重新获取。有难度，osg的几个函数不知道怎么应用，暂时放下。双击投影信息获取wkt格式投影，解决。             三：像素值超过255的压缩处理处理，gdal里面有一个压缩的函数，translate中出现过，解决                 四：显示彩色影像，完成了波段的选择。下一步将波段参数传入到缩略图和放大图中进行显示，解决。                 五：对于不同的波段数据类型，都予以支持，在base1中，读取时用的rasterio来读取的，可以自动进行类型转换，但是base2中，readblock并不支持直接的类型转换，所以必须先用void指针类型来读取数据，然后根据不同的数据类型对void指针进行强制转换，程序中支持8,16,32，64位简单数据类型，对于complex类型的由于不常见暂不支持，如果到时候有需要，可以自行添加支持。 解决      六：界面美化，几个窗口的位置摆放完成。 七：nblocksize过下或者过大的问题。解决  。 七：波段名字显示问题。投影数据显示问题    八：提示建立金字塔文件 ，解决。 九：错误检测，研究c++异常处理，决定使用return和assert来完成程序中出现的错误， 十：对于先前转换错误，失败原因是gdaltranslate和warp不能支持envi格式，即不能支持virtual dataset，所以以后还是都写成tif格式比较好。解决。 十一：内存释放。界面关闭。

今天下午调试了一下午，最后的问题出在调用gdaltranslate和warp函数上面，原因是对于堆栈的问题不敏感，导致反复调试也不知道怎么回事儿，下面贴出一段代码

  char **p=new char *[7];
 for (int i=0;i<7;i++)
 {
  p[i]=new char[512];
 }
 sprintf(p[0],"%s","gdalwarp");
 sprintf(p[1],"%s","-s_srs");
 sprintf(p[2],"%s",pszGCPProjection);
 sprintf(p[3],"%s","-t_srs");
 sprintf(p[4],"%s",pszGCPProjection);
 sprintf(p[5],"%s",tempfilename);
 sprintf(p[6],"%s",destFileName);
 Mywarp(7,p);

 

16：用真实数据测试。可以满足正确性，有空了再把ui程序员应该注意的地方看一下

17：把影像切割的软件重新做一下。一：打开影像，首先看有没有金字塔影像，如果没有，则建立金字塔影像。二：显示图像。对缩放等内容进行优化。三：存储。

18：把重投影再做一下。重点研究ogr、virtual dataset; vrd可以做的工作有。获取virtual daset，通过virtual dataset获得源文件的projection。创建新的projection，更改projection的一些参数（投影，地理参考，仿射变换参数。）。首先更改地理参考类型、更改投影类型。设置virtual dataset投影类型。获取源文件仿射变换参数，转换到新文件下得仿射变换参数，然后设置virtual dataset的放射变换参数。转换失败后提示错误信息。用utm实验验证，本程序可以得出正确的结果（与envi作比较）。

19：修改影像切割工具，修改：在显示影像之前，选择显示波段，同时选择目标文件的波段。这样可以减去那些无用波段。

 ps：到目前为止，暑假的所有任务完成。congratulations。