pktools 2.6.7
Processing Kernel for geospatial data
pkfilter.cc
1/**********************************************************************
2pkfilter.cc: program to filter raster images: median, min/max, morphological, filtering
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along 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/******************************************************************************/
225using namespace std;
226/*------------------
227 Main procedure
228 ----------------*/
229int 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
528 filter2d::Filter2d filter2d;
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
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}
int nrOfRow(void) const
Get the number of rows of this dataset.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
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.
int nrOfBand(void) const
Get the number of bands of this dataset.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
std::string getInterleave() const
Get the band coding (interleave)
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...
std::string getImageType() const
Get the image type (implemented as the driver description)
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 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
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an 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...
void close(void)
Close the image.
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...
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)
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