#include "stdafx.h"
#include "fangshibo.h"

#include <vector>


#include<math.h>
///gdal头文件
#include "..\include\gdal.h"
#include "..\include\gdal_priv.h"
#include "..\include\ogr_srs_api.h"
#include "..\include\cpl_string.h"
#include "..\include\cpl_conv.h"
#pragma comment (lib,"..\lib\gdal_i.lib")
/////////////////////////////////////////////////////////////////////////////


#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// The one and only application object

CWinApp theApp;
using namespace std;

///////////////////////
void  MaxNDVI(vector<CString> files,CString maxNDVfile)
{  
//    GDALAllRegister();
   vector<GDALDataset *> datasets;  
       //打开文件
       for(int i=0;i<files.size();i++)
       {  
           CString filename=files[i];          
           GDALDataset *tmpDataset = (GDALDataset *) GDALOpen( files[i], GA_ReadOnly);
           if(tmpDataset!=NULL)
           {              
               datasets.push_back(tmpDataset);              
           }
           else
           {
               fprintf( stderr, "GDALOpen failed - %dn%sn", CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
           }  
           tmpDataset=NULL;//以后再释放
       }
       ////读取数据咯，按行来读取
       //GDALDir
       if(datasets.size()==0)return;
       GDALDataset *tmpDataset=datasets[0];
       if(tmpDataset!=NULL)
       {  
           GDALDriver * driver=NULL;
           int index = maxNDVfile.ReverseFind('.');
           if (index < 0)
               return ;
           if (index == maxNDVfile.GetLength()-1)
               return ;
           CString suffix = maxNDVfile.Right(maxNDVfile.GetLength()-1-index);          
           suffix.MakeLower();
           if (suffix == "bmp")
               driver = GetGDALDriverManager()->GetDriverByName("BMP");
           else if (suffix == "jpg")
               driver = GetGDALDriverManager()->GetDriverByName("JPEG");
           else if (suffix == "tif")
               driver = GetGDALDriverManager()->GetDriverByName("GTiff");
           else if (suffix == "img")
               driver = GetGDALDriverManager()->GetDriverByName("HFA");
           else if (suffix == "bt")
               driver = GetGDALDriverManager()->GetDriverByName("BT");
           else if (suffix == "ecw")
               driver = GetGDALDriverManager()->GetDriverByName("ECW");
           else if (suffix == "fits")
               driver = GetGDALDriverManager()->GetDriverByName("FITS");
           else if (suffix == "gif")
               driver = GetGDALDriverManager()->GetDriverByName("GIF");
           else if (suffix == "hdf")
               driver = GetGDALDriverManager()->GetDriverByName("HDF4");
           else if (suffix == "hdr")
               driver = GetGDALDriverManager()->GetDriverByName("EHdr");
         
 
           int w=tmpDataset->GetRasterXSize();
           int h=tmpDataset->GetRasterYSize();          
           GDALDataset *maxNDV=driver->Create(maxNDVfile,w,h,1,GDT_Float32,NULL);      
           int xOff=0;
           int yOff=0;
           int width=w;
           int height=1;//一行一行地读取
           vector<float*> bufs;
           for(int i=0;i<datasets.size();i++)
           {
               float *buf=new float[width*height];
               bufs.push_back(buf);
               buf=NULL;
           }      
           float *newbuf=new float[width*height];
           GDALRasterBand * newpoband=maxNDV->GetRasterBand(1);      
           for(int j=0;j<h;j++)
           {              
               for(int i=0;i<datasets.size();i++)
               {
                   GDALDataset *tmpDt=datasets[i];
                   GDALRasterBand * poband=tmpDt->GetRasterBand(1);              
                   float *buf=bufs[i];
                   xOff=0;
                   yOff=j;
                   poband->RasterIO(GF_Read,xOff,yOff,width,height,buf,width,height,GDT_Float32,0,0);
                   buf=NULL;
                   tmpDt=NULL;//在后面再释放,因为很多个指针指向同一一段内存，不能随便释放内存哦
               }          
               //////////
               //在这可以开始运算咯
               for(int k=0;k<w;k++)//对于每个像元
               {
                   //下面这段代码以后要支持均值，极差，方差
                   for(int kk=0;kk<bufs.size();kk++)
                   {
                       float *buf=bufs[kk];
                       if(kk==0)//初始化
                       {
                           newbuf[k]=buf[k];
                       }
                       if(buf[k]>newbuf[k])//最大的,稍微改下就是最小值，
                       {
                           newbuf[k]=buf[k];
                       }                  
                   }
               }
               //////////写入数据
               newpoband->RasterIO(GF_Write,xOff,yOff,width,height,newbuf,width,height,GDT_Float32,0,0);
           }
           delete newbuf;
           newbuf=NULL;
           for( i=0;i<bufs.size();i++)
           {
               if(bufs[i]!=NULL)
                   delete bufs[i];          
           }
           bufs.resize(0);