#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <sys/stat.h>
#include <math.h>
#include <algorithm>
#include <string.h>
#include <stdint.h>
#include <iomanip>
#include <gdal_priv.h>
#include <cpl_conv.h>
#include <ogr_spatialref.h>
using namespace std;

int main(int argc, char* argv[])
{
if(argc != 4){cout<<"wrong arguments. need file"<<endl; exit(1);}
string file = argv[1];
int zone = atoi (argv[2]);
string id = argv[3];

GDALAllRegister();
GDALDataset  *INGDAL;
GDALRasterBand  *INBAND;
GDALDriver *OUTDRIVER;
GDALDataset *OUTGDAL;
GDALRasterBand *OUTBAND;
OGRSpatialReference oSRS;
char *OUTPRJ = NULL;
char **papszOptions = NULL;
OUTDRIVER = GetGDALDriverManager()->GetDriverByName("GTiff"); if( OUTDRIVER == NULL ) {cout << "no driver" << endl; exit( 1 );};
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.SetUTM( zone, TRUE );
oSRS.exportToWkt( &OUTPRJ );
papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "LZW");
double GeoTransform[6];

int i, k, j, x, y, xloc, yloc;
double sum, count;
int mean, std, newmin, newmax, val;
//arrays

int xsize=(15*6+3)*20;
int ysize =(15*6+3)*20;

string filename = "thumbs_monthly_red_10_s/"+id+"/"+file+"_temp.tif";
cout<<filename<<endl;
INGDAL = (GDALDataset *) GDALOpen( filename.c_str(), GA_ReadOnly ); if(INGDAL==NULL){cout << file << " not found" << endl; exit(1);}
INGDAL->GetGeoTransform(GeoTransform);
//int xsize=INBAND->GetXSize(); //if(xsize != (15+3)*10){cout<<(int)xsize<<"wrong size"<<(15+3)*10<<endl;exit(1);}
//int ysize=INBAND->GetYSize(); //if(ysize != (15+3)*10){cout<<ysize<<"wrong size"<<(15+3)*10<<endl;exit(1);}
uint8_t bands[3][ysize][xsize];
uint8_t box[ysize][xsize];
INBAND = INGDAL->GetRasterBand(1); if(INBAND==NULL){cout << "band not found" << endl; exit(1);}
INBAND->RasterIO(GF_Read, 0, 0, xsize, ysize, bands[0], xsize, ysize, GDT_Byte, 0, 0); 
INBAND = INGDAL->GetRasterBand(2); if(INBAND==NULL){cout << "band not found" << endl; exit(1);}
INBAND->RasterIO(GF_Read, 0, 0, xsize, ysize, bands[1], xsize, ysize, GDT_Byte, 0, 0); 
INBAND = INGDAL->GetRasterBand(3); if(INBAND==NULL){cout << "band not found" << endl; exit(1);}
INBAND->RasterIO(GF_Read, 0, 0, xsize, ysize, bands[2], xsize, ysize, GDT_Byte, 0, 0); 
GDALClose(INGDAL);

//cout<<xsize<<','<<ysize<<endl;
//cout<<15*10-1<<endl;
//cout<<(15+3)*10<<endl;
for(y=0; y<ysize; y++) {for(x=0; x<xsize; x++) {
if (((y==(15*60)-1 or y==15*60 or y==15*60+1 or y==(15)*60+60-1 or y==(15)*60+60 or y==(15)*60+60+1) and (x>=15*60-1 and x<=(15)*60+60+1)) or ((y>=15*60-1 and y<=15*60+3*20+1) and (x==15*60-1 or x==15*60+3*20 or x==15*60 or x==15*60+3*20-1 or x==15*60+3*20+1 or x==15*60+1))){
bands[0][y][x]=255;
bands[1][y][x]=0;
bands[2][y][x]=0;
}
//if (((y==(15*60)-1 or y==(15)*60+60 or y==15*60 or y==(15)*60+60-1) and (x>=15*60-1 and x<=(15)*60+60) and x%3==0) or ((y>=15*60-1 and y<=15*60+3*20 and y%3==0) and (x==15*60-1 or x==15*60+3*20 or x==15*60 or x==15*60+3*20-1))){
//bands[0][y][x]=255;
//bands[1][y][x]=0;
//bands[2][y][x]=0;
//}
if((y>=(15)*60-1 and y<=(15)*60+3*20 and y%3==0) and (x==15*60+1*20-1 or x==15*60+2*20 or x==15*60+1*20 or x==15*60+2*20-1)){
  bands[0][y][x]=255;
  bands[1][y][x]=0;
  bands[2][y][x]=0;
}
if((x>=(15)*60-1 and x<=(15)*60+3*20 and x%3==0) and (y==15*60+1*20-1 or y==15*60+2*20 or y==15*60+1*20 or y==15*60+2*20-1)){
  bands[0][y][x]=255;
  bands[1][y][x]=0;
  bands[2][y][x]=0;
}
}}

filename = "thumbs_monthly_red_10_s/"+id+"/"+file+"_temp3.tif";
OUTGDAL = OUTDRIVER->Create(filename.c_str(), xsize, ysize, 3, GDT_Byte, papszOptions );
OUTGDAL->SetGeoTransform(GeoTransform); OUTGDAL->SetProjection(OUTPRJ); 
OUTBAND = OUTGDAL->GetRasterBand(1); OUTBAND->RasterIO( GF_Write, 0, 0, xsize, ysize, bands[0], xsize, ysize, GDT_Byte, 0, 0 ); 
OUTBAND = OUTGDAL->GetRasterBand(2); OUTBAND->RasterIO( GF_Write, 0, 0, xsize, ysize, bands[1], xsize, ysize, GDT_Byte, 0, 0 ); 
OUTBAND = OUTGDAL->GetRasterBand(3); OUTBAND->RasterIO( GF_Write, 0, 0, xsize, ysize, bands[2], xsize, ysize, GDT_Byte, 0, 0 ); 
GDALClose((GDALDatasetH)OUTGDAL);
return(0);
}