pktools  2.6.7
Processing Kernel for geospatial data
pkfilter.cc
1 /**********************************************************************
2 pkfilter.cc: program to filter raster images: median, min/max, morphological, filtering
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 <fstream>
24 #include <math.h>
25 #include <sys/types.h>
26 #include <stdio.h>
27 #include "base/Optionpk.h"
28 #include "base/Vector2d.h"
29 #include "algorithms/Filter2d.h"
30 #include "algorithms/Filter.h"
31 #include "fileclasses/FileReaderAscii.h"
32 #include "imageclasses/ImgReaderGdal.h"
33 #include "imageclasses/ImgWriterGdal.h"
34 #include "algorithms/StatFactory.h"
35 
36 /******************************************************************************/
225 using namespace std;
226 /*------------------
227  Main procedure
228  ----------------*/
229 int main(int argc,char **argv) {
230  Optionpk<std::string> input_opt("i","input","input image file");
231  Optionpk<std::string> output_opt("o", "output", "Output image file");
232  // Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
233  Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
234  // Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation (North=0, East=90, South=180, West=270).");
235  Optionpk<std::string> method_opt("f", "filter", "filter function (nvalid, median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog (central pixel must be different than all other pixels within window), sauvola, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), density, countid, mode (majority voting), only for classes), smoothnodata (smooth nodata values only) values, ismin, ismax, order (rank pixels in order), stdev, mrf, dwt, dwti, dwt_cut, dwt_cut_from, scramble, shift, savgolay, percentile, proportion)");
236  Optionpk<std::string> resample_opt("r", "resampling-method", "Resampling method for shifting operation (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
237  Optionpk<double> dimX_opt("dx", "dx", "filter kernel size in x, use odd values only", 3);
238  Optionpk<double> dimY_opt("dy", "dy", "filter kernel size in y, use odd values only", 3);
239  Optionpk<int> dimZ_opt("dz", "dz", "filter kernel size in z (spectral/temporal dimension), must be odd (example: 3).. Set dz>0 if 1-D filter must be used in band domain");
240  Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
241  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);
242  Optionpk<int> savgolay_nl_opt("nl", "nl", "Number of leftward (past) data points used in Savitzky-Golay filter)", 2);
243  Optionpk<int> savgolay_nr_opt("nr", "nr", "Number of rightward (future) data points used in Savitzky-Golay filter)", 2);
244  Optionpk<int> savgolay_ld_opt("ld", "ld", "order of the derivative desired in Savitzky-Golay filter (e.g., ld=0 for smoothed function)", 0);
245  Optionpk<int> savgolay_m_opt("m", "m", "order of the smoothing polynomial in Savitzky-Golay filter, also equal to the highest conserved moment; usual values are m = 2 or m = 4)", 2);
246  Optionpk<short> class_opt("class", "class", "class value(s) to use for density, erosion, dilation, openening and closing, thresholding");
247  Optionpk<double> threshold_opt("t", "threshold", "threshold value(s) to use for threshold filter (one for each class), or threshold to cut for dwt_cut (use 0 to keep all) or dwt_cut_from, or sigma for shift", 0);
248  Optionpk<double> nodata_opt("nodata", "nodata", "nodata value(s) (used for smoothnodata filter)");
249  Optionpk<std::string> tap_opt("tap", "tap", "text file containing taps used for spatial filtering (from ul to lr). Use dimX and dimY to specify tap dimensions in x and y. Leave empty for not using taps");
250  Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering");
251  Optionpk<string> padding_opt("pad","pad", "Padding method for filtering (how to handle edge effects). Choose between: symmetric, replicate, circular, zero (pad with 0).", "symmetric");
252  Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
253  Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing spectral response functions (two columns: wavelength response)");
254  Optionpk<double> wavelengthIn_opt("win", "wavelengthIn", "list of wavelengths in input spectrum (-win band1 -win band2 ...)");
255  Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-wout band1 -wout band2 ...)");
256  Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima");
257  Optionpk<std::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","");
258  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff");
259  Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid). Use none to omit color table");
260  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
261  Optionpk<short> down_opt("d", "down", "down sampling factor. Use value 1 for no downsampling). Use value n>1 for downsampling (aggregation)", 1);
262  Optionpk<string> beta_opt("beta", "beta", "ASCII file with beta for each class transition in Markov Random Field");
263  // Optionpk<double> eps_opt("eps","eps", "error marging for linear feature",0);
264  // Optionpk<bool> l1_opt("l1","l1", "obtain longest object length for linear feature",false);
265  // Optionpk<bool> l2_opt("l2","l2", "obtain shortest object length for linear feature",false,2);
266  // Optionpk<bool> a1_opt("a1","a1", "obtain angle found for longest object length for linear feature",false);
267  // Optionpk<bool> a2_opt("a2","a2", "obtain angle found for shortest object length for linear feature",false);
268  Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
269 
270  resample_opt.setHide(1);
271  option_opt.setHide(1);
272  wavelet_type_opt.setHide(1);
273  family_opt.setHide(1);
274  savgolay_nl_opt.setHide(1);
275  savgolay_nr_opt.setHide(1);
276  savgolay_ld_opt.setHide(1);
277  savgolay_m_opt.setHide(1);
278  class_opt.setHide(1);
279  threshold_opt.setHide(1);
280  tap_opt.setHide(1);
281  tapz_opt.setHide(1);
282  padding_opt.setHide(1);
283  wavelengthIn_opt.setHide(1);
284  wavelengthOut_opt.setHide(1);
285  down_opt.setHide(1);
286  beta_opt.setHide(1);
287  // eps_opt.setHide(1);
288  // l1_opt.setHide(1);
289  // l2_opt.setHide(1);
290  // a1_opt.setHide(1);
291  // a2_opt.setHide(1);
292  interpolationType_opt.setHide(1);
293  otype_opt.setHide(1);
294  oformat_opt.setHide(1);
295  colorTable_opt.setHide(1);
296  disc_opt.setHide(1);
297 
298  bool doProcess;//stop process when program was invoked with help option (-h --help)
299  try{
300  doProcess=input_opt.retrieveOption(argc,argv);
301  output_opt.retrieveOption(argc,argv);
302  // tmpdir_opt.retrieveOption(argc,argv);
303  // angle_opt.retrieveOption(argc,argv);
304  method_opt.retrieveOption(argc,argv);
305  srf_opt.retrieveOption(argc,argv);
306  fwhm_opt.retrieveOption(argc,argv);
307  dimX_opt.retrieveOption(argc,argv);
308  dimY_opt.retrieveOption(argc,argv);
309  dimZ_opt.retrieveOption(argc,argv);
310  nodata_opt.retrieveOption(argc,argv);
311  resample_opt.retrieveOption(argc,argv);
312  option_opt.retrieveOption(argc,argv);
313  wavelet_type_opt.retrieveOption(argc,argv);
314  family_opt.retrieveOption(argc,argv);
315  savgolay_nl_opt.retrieveOption(argc,argv);
316  savgolay_nr_opt.retrieveOption(argc,argv);
317  savgolay_ld_opt.retrieveOption(argc,argv);
318  savgolay_m_opt.retrieveOption(argc,argv);
319  class_opt.retrieveOption(argc,argv);
320  threshold_opt.retrieveOption(argc,argv);
321  tap_opt.retrieveOption(argc,argv);
322  tapz_opt.retrieveOption(argc,argv);
323  padding_opt.retrieveOption(argc,argv);
324  wavelengthIn_opt.retrieveOption(argc,argv);
325  wavelengthOut_opt.retrieveOption(argc,argv);
326  down_opt.retrieveOption(argc,argv);
327  beta_opt.retrieveOption(argc,argv);
328  // eps_opt.retrieveOption(argc,argv);
329  // l1_opt.retrieveOption(argc,argv);
330  // l2_opt.retrieveOption(argc,argv);
331  // a1_opt.retrieveOption(argc,argv);
332  // a2_opt.retrieveOption(argc,argv);
333  interpolationType_opt.retrieveOption(argc,argv);
334  otype_opt.retrieveOption(argc,argv);
335  oformat_opt.retrieveOption(argc,argv);
336  colorTable_opt.retrieveOption(argc,argv);
337  disc_opt.retrieveOption(argc,argv);
338  verbose_opt.retrieveOption(argc,argv);
339  }
340  catch(string predefinedString){
341  std::cout << predefinedString << std::endl;
342  exit(0);
343  }
344  if(!doProcess){
345  cout << endl;
346  cout << "Usage: pkfilter -i input -o output [-f filter | -perc value | -srf file [-srf file]* -win wavelength [-win wavelength]* | -wout wavelength -fwhm value [-wout wavelength -fwhm value]* -win wavelength [-win wavelength]*]" << endl;
347  cout << endl;
348  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
349  exit(0);//help was invoked, stop processing
350  }
351 
352  //not implemented yet, must debug first...
353  vector<double> angle_opt;
354 
355  ImgReaderGdal input;
356  ImgWriterGdal output;
357  if(input_opt.empty()){
358  cerr << "Error: no input file selected, use option -i" << endl;
359  exit(1);
360  }
361  if(output_opt.empty()){
362  cerr << "Error: no output file selected, use option -o" << endl;
363  exit(1);
364  }
365  input.open(input_opt[0]);
366  GDALDataType theType=GDT_Unknown;
367  if(verbose_opt[0])
368  cout << "possible output data types: ";
369  for(int iType = 0; iType < GDT_TypeCount; ++iType){
370  if(verbose_opt[0])
371  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
372  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
373  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
374  otype_opt[0].c_str()))
375  theType=(GDALDataType) iType;
376  }
377  if(theType==GDT_Unknown)
378  theType=input.getDataType();
379 
380  if(verbose_opt[0])
381  std::cout << std::endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
382 
383  string imageType;//=input.getImageType();
384  if(oformat_opt.size())
385  imageType=oformat_opt[0];
386 
387  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
388  string theInterleave="INTERLEAVE=";
389  theInterleave+=input.getInterleave();
390  option_opt.push_back(theInterleave);
391  }
392 
393  try{
394  string errorString;
395  int nband=input.nrOfBand();
396 
397  if(fwhm_opt.size())
398  nband=fwhm_opt.size();
399  else if(srf_opt.size())
400  nband=srf_opt.size();
401  else if(tap_opt.size()||tapz_opt.size())
402  nband=input.nrOfBand();
403  else{
404  if(method_opt.empty()){
405  errorString="Error: no filter selected, use option -f";
406  throw(errorString);
407  }
408  switch(filter2d::Filter2d::getFilterType(method_opt[0])){
409  case(filter2d::dilate):
410  case(filter2d::erode):
411  case(filter2d::close):
412  case(filter2d::open):
413  case(filter2d::smooth):
414  //implemented in spectral/temporal domain (dimZ>1) and spatial domain
415  if(dimZ_opt.size())
416  assert(dimZ_opt[0]>1);
417  nband=input.nrOfBand();
418  break;
419  case(filter2d::dwt):
420  case(filter2d::dwti):
421  case(filter2d::dwt_cut):
422  case(filter2d::smoothnodata):
423  //implemented in spectral/temporal/spatial domain and nband always input.nrOfBand()
424  nband=input.nrOfBand();
425  break;
426  case(filter2d::savgolay):
427  nband=input.nrOfBand();
428  if(dimZ_opt.empty())
429  dimZ_opt.push_back(1);
430  case(filter2d::dwt_cut_from):
431  //only implemented in spectral/temporal domain
432  if(dimZ_opt.size()){
433  nband=input.nrOfBand();
434  assert(threshold_opt.size());
435  }
436  else{
437  errorString="filter not implemented in spatial domain";
438  throw(errorString);
439  }
440  break;
441  case(filter2d::mrf)://deliberate fall through
442  assert(class_opt.size()>1);
443  if(verbose_opt[0])
444  std::cout << "opening output image " << output_opt[0] << std::endl;
445  nband=class_opt.size();
446  case(filter2d::ismin):
447  case(filter2d::ismax):
448  case(filter2d::shift):
449  case(filter2d::scramble):
450  case(filter2d::mode):
451  case(filter2d::sobelx):
452  case(filter2d::sobely):
453  case(filter2d::sobelxy):
454  case(filter2d::countid):
455  case(filter2d::order):
456  case(filter2d::density):
457  case(filter2d::homog):
458  case(filter2d::heterog):
459  case(filter2d::sauvola):
460  //only implemented in spatial domain
461  if(dimZ_opt.size()){
462  errorString="filter not implemented in spectral/temporal domain";
463  throw(errorString);
464  }
465  break;
466  // case(filter2d::percentile):
467  // //implemented in spectral/temporal/spatial domain and nband 1 if dimZ>0
468  // if(dimZ_opt.size()){
469  // dimZ_opt[0]=1;
470  // nband=1;
471  // }
472  // else
473  // nband=input.nrOfBand();
474  // break;
475  case(filter2d::sum):
476  case(filter2d::mean):
477  case(filter2d::min):
478  case(filter2d::max):
479  case(filter2d::var):
480  case(filter2d::stdev):
481  case(filter2d::nvalid):
482  case(filter2d::median):
483  case(filter2d::percentile):
484  case(filter2d::proportion):
485  //implemented in spectral/temporal/spatial domain and nband 1 if dimZ==1
486  if(dimZ_opt.size()==1)
487  if(dimZ_opt[0]==1)
488  nband=1;
489  else
490  nband=input.nrOfBand();
491  break;
492  default:
493  errorString="filter not implemented";
494  throw(errorString);
495  break;
496  }
497  }
498  std::cout << "opening output image " << output_opt[0] << " with " << nband << " bands" << std::endl;
499  output.open(output_opt[0],(input.nrOfCol()+down_opt[0]-1)/down_opt[0],(input.nrOfRow()+down_opt[0]-1)/down_opt[0],nband,theType,imageType,option_opt);
500  }
501  catch(string errorstring){
502  cerr << errorstring << endl;
503  exit(1);
504  }
505  output.setProjection(input.getProjection());
506  double gt[6];
507  input.getGeoTransform(gt);
508  gt[1]*=down_opt[0];//dx
509  gt[5]*=down_opt[0];//dy
510  output.setGeoTransform(gt);
511 
512  if(colorTable_opt.size()){
513  if(colorTable_opt[0]!="none"){
514  if(verbose_opt[0])
515  cout << "set colortable " << colorTable_opt[0] << endl;
516  assert(output.getDataType()==GDT_Byte);
517  output.setColorTable(colorTable_opt[0]);
518  }
519  }
520  else if(input.getColorTable()!=NULL)
521  output.setColorTable(input.getColorTable());
522 
523  if(nodata_opt.size()){
524  for(int iband=0;iband<output.nrOfBand();++iband)
525  output.GDALSetNoDataValue(nodata_opt[0],iband);
526  }
527 
529  filter::Filter filter1d;
530  if(verbose_opt[0])
531  cout << "Set padding to " << padding_opt[0] << endl;
532  filter1d.setPadding(padding_opt[0]);
533  if(class_opt.size()){
534  if(verbose_opt[0])
535  std::cout<< "class values: ";
536  for(int iclass=0;iclass<class_opt.size();++iclass){
537  if(!dimZ_opt.size())
538  filter2d.pushClass(class_opt[iclass]);
539  else
540  filter1d.pushClass(class_opt[iclass]);
541  if(verbose_opt[0])
542  std::cout<< class_opt[iclass] << " ";
543  }
544  if(verbose_opt[0])
545  std::cout<< std::endl;
546  }
547 
548  if(nodata_opt.size()){
549  if(verbose_opt[0])
550  std::cout<< "mask values: ";
551  for(int imask=0;imask<nodata_opt.size();++imask){
552  if(verbose_opt[0])
553  std::cout<< nodata_opt[imask] << " ";
554  filter1d.pushNoDataValue(nodata_opt[imask]);
555  filter2d.pushNoDataValue(nodata_opt[imask]);
556  }
557  if(verbose_opt[0])
558  std::cout<< std::endl;
559  }
560  filter1d.setThresholds(threshold_opt);
561  filter2d.setThresholds(threshold_opt);
562 
563  if(tap_opt.size()){
564  ifstream tapfile(tap_opt[0].c_str());
565  assert(tapfile);
566  Vector2d<double> taps(dimY_opt[0],dimX_opt[0]);
567 
568  for(int j=0;j<dimY_opt[0];++j){
569  for(int i=0;i<dimX_opt[0];++i){
570  tapfile >> taps[j][i];
571  }
572  }
573  if(verbose_opt[0]){
574  std::cout << "taps: ";
575  for(int j=0;j<dimY_opt[0];++j){
576  for(int i=0;i<dimX_opt[0];++i){
577  std::cout<< taps[j][i] << " ";
578  }
579  std::cout<< std::endl;
580  }
581  }
582  filter2d.setTaps(taps);
583  try{
584  filter2d.filter(input,output);
585  }
586  catch(string errorstring){
587  cerr << errorstring << endl;
588  }
589  tapfile.close();
590  }
591  else if(tapz_opt.size()){
592  if(verbose_opt[0]){
593  std::cout << "taps: ";
594  for(int itap=0;itap<tapz_opt.size();++itap)
595  std::cout<< tapz_opt[itap] << " ";
596  std::cout<< std::endl;
597  }
598  filter1d.setTaps(tapz_opt);
599  filter1d.filter(input,output);
600  }
601  else if(fwhm_opt.size()){
602  if(verbose_opt[0])
603  std::cout << "spectral filtering to " << fwhm_opt.size() << " bands with provided fwhm " << std::endl;
604  assert(wavelengthOut_opt.size()==fwhm_opt.size());
605  assert(wavelengthIn_opt.size());
606 
607  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
608  Vector2d<double> lineOutput(wavelengthOut_opt.size(),input.nrOfCol());
609  const char* pszMessage;
610  void* pProgressArg=NULL;
611  GDALProgressFunc pfnProgress=GDALTermProgress;
612  double progress=0;
613  pfnProgress(progress,pszMessage,pProgressArg);
614  for(int y=0;y<input.nrOfRow();++y){
615  if((y+1+down_opt[0]/2)%down_opt[0])
616  continue;
617  for(int iband=0;iband<input.nrOfBand();++iband)
618  input.readData(lineInput[iband],y,iband);
619  filter1d.applyFwhm<double>(wavelengthIn_opt,lineInput,wavelengthOut_opt,fwhm_opt, interpolationType_opt[0], lineOutput, down_opt[0], verbose_opt[0]);
620  for(int iband=0;iband<output.nrOfBand();++iband){
621  try{
622  output.writeData(lineOutput[iband],y/down_opt[0],iband);
623  }
624  catch(string errorstring){
625  cerr << errorstring << "in band " << iband << ", line " << y << endl;
626  }
627  }
628  progress=(1.0+y)/output.nrOfRow();
629  pfnProgress(progress,pszMessage,pProgressArg);
630  }
631  }
632  else if(srf_opt.size()){
633  if(verbose_opt[0])
634  std::cout << "spectral filtering to " << srf_opt.size() << " bands with provided SRF " << std::endl;
635  assert(wavelengthIn_opt.size());
636  vector< Vector2d<double> > srf(srf_opt.size());//[0] srf_nr, [1]: wavelength, [2]: response
637  ifstream srfFile;
638  for(int isrf=0;isrf<srf_opt.size();++isrf){
639  srf[isrf].resize(2);
640  srfFile.open(srf_opt[isrf].c_str());
641  double v;
642  //add 0 to make sure srf is 0 at boundaries after interpolation step
643  srf[isrf][0].push_back(0);
644  srf[isrf][1].push_back(0);
645  srf[isrf][0].push_back(1);
646  srf[isrf][1].push_back(0);
647  while(srfFile >> v){
648  srf[isrf][0].push_back(v);
649  srfFile >> v;
650  srf[isrf][1].push_back(v);
651  }
652  srfFile.close();
653  //add 0 to make sure srf[isrf] is 0 at boundaries after interpolation step
654  srf[isrf][0].push_back(srf[isrf][0].back()+1);
655  srf[isrf][1].push_back(0);
656  srf[isrf][0].push_back(srf[isrf][0].back()+1);
657  srf[isrf][1].push_back(0);
658  if(verbose_opt[0])
659  cout << "srf file details: " << srf[isrf][0].size() << " wavelengths defined" << endl;
660  }
661  assert(output.nrOfBand()==srf.size());
662  double centreWavelength=0;
663  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
664  const char* pszMessage;
665  void* pProgressArg=NULL;
666  GDALProgressFunc pfnProgress=GDALTermProgress;
667  double progress=0;
668  pfnProgress(progress,pszMessage,pProgressArg);
669  for(int y=0;y<input.nrOfRow();++y){
670  if((y+1+down_opt[0]/2)%down_opt[0])
671  continue;
672  for(int iband=0;iband<input.nrOfBand();++iband)
673  input.readData(lineInput[iband],y,iband);
674  for(int isrf=0;isrf<srf.size();++isrf){
675  vector<double> lineOutput(output.nrOfCol());
676  double delta=1.0;
677  bool normalize=true;
678  centreWavelength=filter1d.applySrf<double>(wavelengthIn_opt,lineInput,srf[isrf], interpolationType_opt[0], lineOutput, delta, normalize);
679  if(verbose_opt[0])
680  std::cout << "centre wavelength srf " << isrf << ": " << centreWavelength << std::endl;
681  try{
682  output.writeData(lineOutput,GDT_Float64,y/down_opt[0],isrf);
683  }
684  catch(string errorstring){
685  cerr << errorstring << "in srf " << srf_opt[isrf] << ", line " << y << endl;
686  }
687 
688  }
689  progress=(1.0+y)/output.nrOfRow();
690  pfnProgress(progress,pszMessage,pProgressArg);
691  }
692 
693  }
694  else{
695  switch(filter2d::Filter2d::getFilterType(method_opt[0])){
696  case(filter2d::dilate):
697  if(down_opt[0]!=1){
698  std::cerr << "Error: down option not supported for morphological operator" << std::endl;
699  exit(1);
700  }
701  try{
702  if(dimZ_opt.size()){
703  if(verbose_opt[0])
704  std::cout<< "1-D filtering: dilate" << std::endl;
705  filter1d.morphology(input,output,"dilate",dimZ_opt[0],verbose_opt[0]);
706  }
707  else
708  filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
709  }
710  catch(string errorstring){
711  cerr << errorstring << endl;
712  }
713  break;
714  case(filter2d::erode):
715  if(down_opt[0]!=1){
716  std::cerr << "Error: down option not supported for morphological operator" << std::endl;
717  exit(1);
718  }
719  try{
720  if(dimZ_opt.size()>0){
721  if(verbose_opt[0])
722  std::cout<< "1-D filtering: dilate" << std::endl;
723  filter1d.morphology(input,output,"erode",dimZ_opt[0]);
724  }
725  else{
726  filter2d.morphology(input,output,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
727  }
728  }
729  catch(string errorstring){
730  cerr << errorstring << endl;
731  }
732  break;
733  case(filter2d::close):{//closing
734  if(down_opt[0]!=1){
735  std::cerr << "Error: down option not supported for morphological operator" << std::endl;
736  exit(1);
737  }
738 
739  ImgWriterGdal tmpout;
740  tmpout.open("/vsimem/dilation.tif",input.nrOfCol(),input.nrOfRow(),input.nrOfBand(),input.getDataType(),input.getImageType());
741  try{
742  if(dimZ_opt.size()){
743  filter1d.morphology(input,tmpout,"dilate",dimZ_opt[0]);
744  }
745  else{
746  filter2d.morphology(input,tmpout,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
747  }
748  }
749  catch(std::string errorString){
750  std::cout<< errorString;
751  exit(1);
752  }
753  tmpout.close();
754  ImgReaderGdal tmpin;
755  tmpin.open("/vsimem/dilation.tif");
756  try{
757  if(dimZ_opt.size()){
758  filter1d.morphology(tmpin,output,"erode",dimZ_opt[0]);
759  }
760  else{
761  filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
762  }
763  }
764  catch(string errorstring){
765  cerr << errorstring << endl;
766  }
767  tmpin.close();
768  break;
769  }
770  case(filter2d::open):{//opening
771  if(down_opt[0]!=1){
772  std::cerr << "Error: down option not supported for morphological operator" << std::endl;
773  exit(1);
774  }
775  ImgWriterGdal tmpout;
776  tmpout.open("/vsimem/erosion.tif",input.nrOfCol(),input.nrOfRow(),input.nrOfBand(),input.getDataType(),input.getImageType());
777  try{
778  if(dimZ_opt.size()){
779  filter1d.morphology(input,tmpout,"erode",dimZ_opt[0]);
780  }
781  else{
782  filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
783  }
784  }
785  catch(std::string errorString){
786  std::cout<< errorString;
787  exit(1);
788  }
789  tmpout.close();
790  ImgReaderGdal tmpin;
791  try{
792  tmpin.open("/vsimem/erosion.tif");
793  if(dimZ_opt.size()){
794  filter1d.morphology(tmpin,output,"dilate",dimZ_opt[0]);
795  }
796  else{
797  filter2d.morphology(tmpin,output,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
798  }
799  tmpin.close();
800  tmpout.close();
801  }
802  catch(string errorstring){
803  cerr << errorstring << endl;
804  }
805  break;
806  }
807  case(filter2d::homog):{//spatially homogeneous
808  try{
809  filter2d.doit(input,output,"homog",dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
810  }
811  catch(string errorstring){
812  cerr << errorstring << endl;
813  }
814  break;
815  }
816  case(filter2d::heterog):{//spatially heterogeneous
817  try{
818  filter2d.doit(input,output,"heterog",dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
819  }
820  catch(string errorstring){
821  cerr << errorstring << endl;
822  }
823  break;
824  }
825  case(filter2d::sauvola):{//Implements Sauvola's thresholding method (http://fiji.sc/Auto_Local_Threshold)
826 
827  //test
828  Vector2d<unsigned short> inBuffer;
829  for(int iband=0;iband<input.nrOfBand();++iband){
830  input.readDataBlock(inBuffer,0,input.nrOfCol()-1,0,input.nrOfRow()-1,iband);
831  }
832  try{
833  filter2d.doit(input,output,"sauvola",dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
834  }
835  catch(string errorstring){
836  cerr << errorstring << endl;
837  }
838  break;
839  }
840  case(filter2d::shift):{//shift
841  if(down_opt[0]!=1){
842  std::cerr << "Error: down option not supported for shift operator" << std::endl;
843  exit(1);
844  }
845  assert(input.nrOfBand());
846  assert(input.nrOfCol());
847  assert(input.nrOfRow());
848  try{
849  filter2d.shift(input,output,dimX_opt[0],dimY_opt[0],threshold_opt[0],filter2d::Filter2d::getResampleType(resample_opt[0]));
850  }
851  catch(string errorstring){
852  cerr << errorstring << endl;
853  }
854  break;
855  }
856  // case(filter2d::linearfeature):{
857  // if(down_opt[0]!=1){
858  // std::cerr << "Error: down option not supported for linear feature" << std::endl;
859  // exit(1);
860  // }
861  // assert(input.nrOfBand());
862  // assert(input.nrOfCol());
863  // assert(input.nrOfRow());
864  // float theAngle=361;
865  // if(angle_opt.size())
866  // theAngle=angle_opt[0];
867  // if(verbose_opt[0])
868  // std::cout << "using angle " << theAngle << std::endl;
869  // try{
870  // //using an angle step of 5 degrees and no maximum distance
871  // filter2d.linearFeature(input,output,theAngle,5,0,eps_opt[0],l1_opt[0],a1_opt[0],l2_opt[0],a2_opt[0],0,verbose_opt[0]);
872  // }
873  // catch(string errorstring){
874  // cerr << errorstring << endl;
875  // }
876  // break;
877  // }
878  case(filter2d::mrf):{//Markov Random Field
879  if(verbose_opt[0])
880  std::cout << "Markov Random Field filtering" << std::endl;
881  try{
882  if(beta_opt.size()){
883  //in file: classFrom classTo
884  //in variable: beta[classTo][classFrom]
885  FileReaderAscii betaReader(beta_opt[0]);
886  Vector2d<double> beta(class_opt.size(),class_opt.size());
887  vector<int> cols(class_opt.size());
888  for(int iclass=0;iclass<class_opt.size();++iclass)
889  cols[iclass]=iclass;
890  betaReader.readData(beta,cols);
891  if(verbose_opt[0]){
892  std::cout << "using values for beta:" << std::endl;
893  for(int iclass1=0;iclass1<class_opt.size();++iclass1)
894  std::cout << " " << iclass1 << " (" << class_opt[iclass1] << ")";
895  std::cout << std::endl;
896  for(int iclass1=0;iclass1<class_opt.size();++iclass1){
897  std::cout << iclass1 << " (" << class_opt[iclass1] << ")";
898  for(int iclass2=0;iclass2<class_opt.size();++iclass2)
899  std::cout << " " << beta[iclass2][iclass1] << " (" << class_opt[iclass2] << ")";
900  std::cout << std::endl;
901  }
902  }
903  filter2d.mrf(input, output, dimX_opt[0], dimY_opt[0], beta, true, down_opt[0], verbose_opt[0]);
904  }
905  else
906  filter2d.mrf(input, output, dimX_opt[0], dimY_opt[0], 1, true, down_opt[0], verbose_opt[0]);
907  }
908  catch(string errorstring){
909  cerr << errorstring << endl;
910  }
911  break;
912  }
913  case(filter2d::sobelx):{//Sobel edge detection in X
914  if(down_opt[0]!=1){
915  std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
916  exit(1);
917  }
918  Vector2d<double> theTaps(3,3);
919  theTaps[0][0]=-1.0;
920  theTaps[0][1]=0.0;
921  theTaps[0][2]=1.0;
922  theTaps[1][0]=-2.0;
923  theTaps[1][1]=0.0;
924  theTaps[1][2]=2.0;
925  theTaps[2][0]=-1.0;
926  theTaps[2][1]=0.0;
927  theTaps[2][2]=1.0;
928  filter2d.setTaps(theTaps);
929  try{
930  filter2d.filter(input,output,true,true);//absolute and normalize
931  }
932  catch(string errorstring){
933  cerr << errorstring << endl;
934  }
935  break;
936  }
937  case(filter2d::sobely):{//Sobel edge detection in Y
938  if(down_opt[0]!=1){
939  std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
940  exit(1);
941  }
942  Vector2d<double> theTaps(3,3);
943  theTaps[0][0]=1.0;
944  theTaps[0][1]=2.0;
945  theTaps[0][2]=1.0;
946  theTaps[1][0]=0.0;
947  theTaps[1][1]=0.0;
948  theTaps[1][2]=0.0;
949  theTaps[2][0]=-1.0;
950  theTaps[2][1]=-2.0;
951  theTaps[2][2]=-1.0;
952  filter2d.setTaps(theTaps);
953  try{
954  filter2d.filter(input,output,true,true);//absolute and normalize
955  }
956  catch(string errorstring){
957  cerr << errorstring << endl;
958  }
959  break;
960  }
961  case(filter2d::sobelxy):{//Sobel edge detection in XY
962  if(down_opt[0]!=1){
963  std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
964  exit(1);
965  }
966  Vector2d<double> theTaps(3,3);
967  theTaps[0][0]=0.0;
968  theTaps[0][1]=1.0;
969  theTaps[0][2]=2.0;
970  theTaps[1][0]=-1.0;
971  theTaps[1][1]=0.0;
972  theTaps[1][2]=1.0;
973  theTaps[2][0]=-2.0;
974  theTaps[2][1]=-1.0;
975  theTaps[2][2]=0.0;
976  filter2d.setTaps(theTaps);
977  try{
978  filter2d.filter(input,output,true,true);//absolute and normalize
979  }
980  catch(string errorstring){
981  cerr << errorstring << endl;
982  }
983  break;
984  }
985  case(filter2d::sobelyx):{//Sobel edge detection in XY
986  if(down_opt[0]!=1){
987  std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
988  exit(1);
989  }
990  Vector2d<double> theTaps(3,3);
991  theTaps[0][0]=2.0;
992  theTaps[0][1]=1.0;
993  theTaps[0][2]=0.0;
994  theTaps[1][0]=1.0;
995  theTaps[1][1]=0.0;
996  theTaps[1][2]=-1.0;
997  theTaps[2][0]=0.0;
998  theTaps[2][1]=-1.0;
999  theTaps[2][2]=-2.0;
1000  filter2d.setTaps(theTaps);
1001  try{
1002  filter2d.filter(input,output,true,true);//absolute and normalize
1003  }
1004  catch(string errorstring){
1005  cerr << errorstring << endl;
1006  }
1007  break;
1008  }
1009  case(filter2d::smooth):{//Smoothing filter
1010  if(down_opt[0]!=1){
1011  std::cerr << "Error: down option not supported for this filter" << std::endl;
1012  exit(1);
1013  }
1014  try{
1015  if(dimZ_opt.size()){
1016  if(verbose_opt[0])
1017  std::cout<< "1-D filtering: smooth" << std::endl;
1018  filter1d.smooth(input,output,dimZ_opt[0]);
1019  }
1020  else{
1021  filter2d.smooth(input,output,dimX_opt[0],dimY_opt[0]);
1022  }
1023  }
1024  catch(string errorstring){
1025  cerr << errorstring << endl;
1026  }
1027  break;
1028  }
1029  case(filter2d::smoothnodata):{//Smoothing filter
1030  if(down_opt[0]!=1){
1031  std::cerr << "Error: down option not supported for this filter" << std::endl;
1032  exit(1);
1033  }
1034  try{
1035  if(dimZ_opt.size()){
1036  if(verbose_opt[0])
1037  std::cout<< "1-D filtering: smooth" << std::endl;
1038  filter1d.smoothNoData(input,interpolationType_opt[0],output);
1039  }
1040  else{
1041  if(verbose_opt[0])
1042  std::cout<< "2-D filtering: smooth" << std::endl;
1043  filter2d.smoothNoData(input,output,dimX_opt[0],dimY_opt[0]);
1044  }
1045  }
1046  catch(string errorstring){
1047  cerr << errorstring << endl;
1048  }
1049  break;
1050  }
1051  case(filter2d::dwt):
1052  if(down_opt[0]!=1){
1053  std::cerr << "Error: down option not supported for this filter" << std::endl;
1054  exit(1);
1055  }
1056  try{
1057  if(dimZ_opt.size()){
1058  if(verbose_opt[0])
1059  std::cout<< "DWT in spectral domain" << std::endl;
1060  filter1d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
1061  }
1062  else
1063  filter2d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
1064  }
1065  catch(string errorstring){
1066  cerr << errorstring << endl;
1067  }
1068  break;
1069  case(filter2d::dwti):
1070  if(down_opt[0]!=1){
1071  std::cerr << "Error: down option not supported for this filter" << std::endl;
1072  exit(1);
1073  }
1074  try{
1075  if(dimZ_opt.size()){
1076  if(verbose_opt[0])
1077  std::cout<< "inverse DWT in spectral domain" << std::endl;
1078  filter1d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
1079  }
1080  else
1081  filter2d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
1082  }
1083  catch(string errorstring){
1084  cerr << errorstring << endl;
1085  }
1086  break;
1087  case(filter2d::dwt_cut):
1088  if(down_opt[0]!=1){
1089  std::cerr << "Error: down option not supported for this filter" << std::endl;
1090  exit(1);
1091  }
1092  if(dimZ_opt.size()){
1093  if(verbose_opt[0])
1094  std::cout<< "DWT approximation in spectral domain" << std::endl;
1095  filter1d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
1096  }
1097  else
1098  filter2d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
1099  break;
1100  case(filter2d::dwt_cut_from):
1101  if(down_opt[0]!=1){
1102  std::cerr << "Error: down option not supported for this filter" << std::endl;
1103  exit(1);
1104  }
1105  try{
1106  if(dimZ_opt.size()){
1107  if(verbose_opt[0])
1108  std::cout<< "DWT approximation in spectral domain" << std::endl;
1109  filter1d.dwtCutFrom(input, output, wavelet_type_opt[0], family_opt[0], static_cast<int>(threshold_opt[0]));
1110  }
1111  else{
1112  string errorString="Error: this filter is not supported in 2D";
1113  throw(errorString);
1114  }
1115  }
1116  catch(string errorstring){
1117  cerr << errorstring << endl;
1118  }
1119  break;
1120  case(filter2d::savgolay):{
1121  assert(savgolay_nl_opt.size());
1122  assert(savgolay_nr_opt.size());
1123  assert(savgolay_ld_opt.size());
1124  assert(savgolay_m_opt.size());
1125  if(verbose_opt[0])
1126  std::cout << "Calculating Savitzky-Golay coefficients: " << endl;
1127  filter1d.getSavGolayCoefficients(tapz_opt, input.nrOfBand(), savgolay_nl_opt[0], savgolay_nr_opt[0], savgolay_ld_opt[0], savgolay_m_opt[0]);
1128  if(verbose_opt[0]){
1129  std::cout << "taps (size is " << tapz_opt.size() << "): ";
1130  for(int itap=0;itap<tapz_opt.size();++itap)
1131  std::cout<< tapz_opt[itap] << " ";
1132  std::cout<< std::endl;
1133  }
1134  filter1d.setTaps(tapz_opt);
1135  filter1d.filter(input,output);
1136  break;
1137  }
1138  case(filter2d::percentile)://deliberate fall through
1139  case(filter2d::threshold)://deliberate fall through
1140  assert(threshold_opt.size());
1141  if(dimZ_opt.size())
1142  filter1d.setThresholds(threshold_opt);
1143  else
1144  filter2d.setThresholds(threshold_opt);
1145  case(filter2d::density)://deliberate fall through
1146  filter2d.setClasses(class_opt);
1147  if(verbose_opt[0])
1148  std::cout << "classes set" << std::endl;
1149  default:
1150  try{
1151  if(dimZ_opt.size()){
1152  if(dimZ_opt[0]==1)
1153  filter1d.stat(input,output,method_opt[0]);
1154  else{
1155  assert(down_opt[0]==1);//not implemented yet...
1156  filter1d.filter(input,output,method_opt[0],dimZ_opt[0]);
1157  }
1158  }
1159  else
1160  filter2d.doit(input,output,method_opt[0],dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
1161  }
1162  catch(string errorstring){
1163  cerr << errorstring << endl;
1164  }
1165  break;
1166  }
1167  }
1168  input.close();
1169  output.close();
1170  return 0;
1171 }
STL namespace.
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
void setColorTable(const std::string &filename, int band=0)
Set the color table using an (ASCII) file with 5 columns (value R G B alpha)
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
int nrOfRow(void) const
Get the number of rows of this dataset.
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
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
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 getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
void close(void)
Close the image.
void readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
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 nrOfBand(void) const
Get the number of bands of this dataset.
std::string getInterleave() const
Get the band coding (interleave)