pktools  2.6.7
Processing Kernel for geospatial data
pkenhance.cc
1 /**********************************************************************
2 pkenhance.cc: program to enhance raster images: histogram matching
3 Copyright (C) 2008-2014 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #include <assert.h>
21 #include <iostream>
22 #include <string>
23 #include "base/Optionpk.h"
24 #include "imageclasses/ImgReaderGdal.h"
25 #include "imageclasses/ImgWriterGdal.h"
26 
27 using namespace std;
28 
29 int main(int argc,char **argv) {
30  Optionpk<std::string> input_opt("i","input","input image file");
31  Optionpk<string> reference_opt("r", "reference", "Reference image file");
32  Optionpk<std::string> output_opt("o", "output", "Output image file");
33  Optionpk<double> minRef_opt("minref", "minref", "Sets minimum for histogram of reference image");
34  Optionpk<double> maxRef_opt("maxref", "maxref", "Sets maximum for histogram of reference image");
35  Optionpk<double> minInput_opt("mininput", "mininput", "Sets minimum for histogram of input image");
36  Optionpk<double> maxInput_opt("maxinput", "maxinput", "Sets maximum for histogram of input image");
37  Optionpk<double> nodata_opt("nodata", "nodata", "Sets no data value(s) for calculations (nodata values in input image)");
38  Optionpk<std::string> method_opt("m", "method", "enhancement method (histmatch)", "histmatch");
39  // Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
40  // Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing moment, see also http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
41  // Optionpk<double> quantize_opt("q", "quantize", "Quantize threshold",0);
42  Optionpk<short> nbin_opt("nbin", "nbin", "Number of bins used in histogram. Use 0 for all input values as integers",0);
43  Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
44  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
45  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
46  Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
47 
48  bool doProcess;//stop process when program was invoked with help option (-h --help)
49  try{
50  doProcess=input_opt.retrieveOption(argc,argv);
51  reference_opt.retrieveOption(argc,argv);
52  output_opt.retrieveOption(argc,argv);
53  minRef_opt.retrieveOption(argc,argv);
54  maxRef_opt.retrieveOption(argc,argv);
55  minInput_opt.retrieveOption(argc,argv);
56  maxInput_opt.retrieveOption(argc,argv);
57  nodata_opt.retrieveOption(argc,argv);
58  method_opt.retrieveOption(argc,argv);
59  nbin_opt.retrieveOption(argc,argv);
60  otype_opt.retrieveOption(argc,argv);
61  oformat_opt.retrieveOption(argc,argv);
62  option_opt.retrieveOption(argc,argv);
63  verbose_opt.retrieveOption(argc,argv);
64  }
65  catch(string predefinedString){
66  std::cout << predefinedString << std::endl;
67  exit(0);
68  }
69  if(!doProcess){
70  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
71  exit(0);//help was invoked, stop processing
72  }
73 
74  ImgReaderGdal inputImg;
75  ImgReaderGdal referenceImg;
76  ImgWriterGdal outputImg;
77  assert(input_opt.size());
78  inputImg.open(input_opt[0]);
79  for(int inodata=0;inodata<nodata_opt.size();++inodata){
80  if(!inodata)
81  inputImg.GDALSetNoDataValue(nodata_opt[0],0);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
82  inputImg.pushNoDataValue(nodata_opt[inodata]);
83  }
84 
85  int nband=inputImg.nrOfBand();
86  GDALDataType theType=GDT_Unknown;
87  if(verbose_opt[0])
88  cout << "possible output data types: ";
89  for(int iType = 0; iType < GDT_TypeCount; ++iType){
90  if(verbose_opt[0])
91  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
92  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
93  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
94  otype_opt[0].c_str()))
95  theType=(GDALDataType) iType;
96  }
97  if(theType==GDT_Unknown)
98  theType=inputImg.getDataType();
99 
100  if(verbose_opt[0])
101  std::cout << std::endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
102 
103  string imageType=inputImg.getImageType();
104  if(oformat_opt.size())
105  imageType=oformat_opt[0];
106 
107  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
108  string theInterleave="INTERLEAVE=";
109  theInterleave+=inputImg.getInterleave();
110  option_opt.push_back(theInterleave);
111  }
112  try{
113  assert(output_opt.size());
114  if(verbose_opt[0])
115  std::cout << "opening output image " << output_opt[0] << std::endl;
116  outputImg.open(output_opt[0],inputImg.nrOfCol(),inputImg.nrOfRow(),inputImg.nrOfBand(),theType,imageType,option_opt);
117  }
118  catch(string errorstring){
119  cout << errorstring << endl;
120  exit(1);
121  }
122 
123  if(method_opt[0]=="histmatch"){
124  assert(reference_opt.size());
125  referenceImg.open(reference_opt[0]);
126  assert(nband==referenceImg.nrOfBand());
127  for(int inodata=0;inodata<nodata_opt.size();++inodata){
128  if(!inodata)
129  referenceImg.GDALSetNoDataValue(nodata_opt[0],0);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
130  referenceImg.pushNoDataValue(nodata_opt[inodata]);
131  }
132  const char* pszMessage;
133  void* pProgressArg=NULL;
134  GDALProgressFunc pfnProgress=GDALTermProgress;
135  double progress=0;
136  pfnProgress(progress,pszMessage,pProgressArg);
137  for(int iband=0;iband<nband;++iband){
138  //calculate histograms
139  unsigned int nbinRef=nbin_opt[0];
140  unsigned int nbinInput=nbin_opt[0];
141  std::vector<double> histRef(nbinRef);
142  std::vector<double> histInput(nbinInput);
143  double minValueRef=0;
144  double maxValueRef=0;
145  double minValueInput=0;
146  double maxValueInput=0;
147  if(minRef_opt.size())
148  minValueRef=minRef_opt[0];
149  if(maxRef_opt.size())
150  maxValueRef=maxRef_opt[0];
151  if(minInput_opt.size())
152  minValueInput=minInput_opt[0];
153  if(maxInput_opt.size())
154  maxValueInput=maxInput_opt[0];
155  unsigned long int nsampleRef=referenceImg.getHistogram(histRef,minValueRef,maxValueRef,nbinRef,iband);
156  unsigned long int nsampleInput=inputImg.getHistogram(histInput,minValueInput,maxValueInput,nbinInput,iband);
157  //create cumulative historgrams
158  for(unsigned int bin=0;bin<nbinRef;++bin){
159  histRef[bin]+=100.0*static_cast<double>(histRef[bin])/static_cast<double>(nsampleRef);
160  }
161  for(unsigned int bin=0;bin<nbinInput;++bin)
162  histInput[bin]+=100.0*static_cast<double>(histInput[bin])/static_cast<double>(nsampleInput);
163  //match histograms
164  vector<double> lineBuffer(inputImg.nrOfCol());
165  for(int irow=0;irow<inputImg.nrOfRow();++irow){
166  inputImg.readData(lineBuffer,GDT_Float64, irow, iband);
167  for(int icol=0;icol<inputImg.nrOfCol();++icol){
168  //find bin in input image histogram
169  int inputBin=static_cast<int>(static_cast<double>(lineBuffer[icol]-minValueInput)/(maxValueInput-minValueInput)*(histInput.size()-1));
170  //find corresponding bin in reference image histogram
171  //todo: optimize with lower_bound?
172  // std::vector<unsigned long int>::const_iterator hit=histRef.begin();
173  int ibin=0;
174  for(ibin=0;ibin<histRef.size();++ibin){
175  if(histRef[ibin]>histInput[inputBin])
176  break;
177  }
178  if(ibin)
179  --ibin;
180  lineBuffer[icol]=(maxValueRef-minValueRef)/(histRef.size()-1)*(ibin)+minValueRef;
181  // std::vector<unsigned long int>::const_iterator it=std::lower_bound(histRef.begin(),histRef.end(),culInput);
182  }
183  outputImg.writeData(lineBuffer,GDT_Float64,irow,iband);
184  progress=(1.0+irow);
185  progress+=(outputImg.nrOfRow()*iband);
186  progress/=outputImg.nrOfBand()*outputImg.nrOfRow();
187  assert(progress>=0);
188  assert(progress<=1);
189  pfnProgress(progress,pszMessage,pProgressArg);
190  }
191  }
192  }
193 
194  inputImg.close();
195  if(method_opt[0]=="histmatch")
196  referenceImg.close();
197  outputImg.close();
198  return 0;
199 }
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
double getHistogram(std::vector< double > &histvector, double &min, double &max, unsigned int &nbin, int theBand=0, bool kde=false)
Calculate the image histogram for a specific band using a defined number of bins and constrained by a...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
std::string getImageType() const
Get the image type (implemented as the driver description)
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
std::string getInterleave() const
Get the band coding (interleave)
void close(void)
Close the image.
CPLErr GDALSetNoDataValue(double noDataValue, int band=0)
Set the GDAL (internal) no data value for this data set. Only a single no data value per band is supp...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98