[pktools] 162/375: worked on pkfilter, changed some filter name to dwtz[i|_cut]
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:10 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit 3ee7575278c912664693271b2e7f6ea4809e3c3f
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Sat Dec 21 23:54:24 2013 +0100
worked on pkfilter, changed some filter name to dwtz[i|_cut]
---
src/algorithms/Filter.cc | 130 +++++++++++++++++++++++++++++++++++++++++++++
src/algorithms/Filter.h | 13 +++--
src/algorithms/Filter2d.cc | 4 +-
src/algorithms/Filter2d.h | 37 ++++++-------
src/apps/pkfilter.cc | 25 +++++----
src/apps/pkfilterascii.cc | 54 ++++++-------------
6 files changed, 189 insertions(+), 74 deletions(-)
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 69812c7..40357ac 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -39,10 +39,105 @@ void filter::Filter::setTaps(const vector<double> &taps)
assert(m_taps.size()%2);
}
+void filter::Filter::dwtForward(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+ Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
+ for(int y=0;y<input.nrOfRow();++y){
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ input.readData(lineInput[iband],GDT_Float64,y,iband);
+ vector<double> pixelInput(input.nrOfBand());
+ for(int x=0;x<input.nrOfCol();++x){
+ pixelInput=lineInput.selectCol(x);
+ dwtForward(pixelInput,wavelet_type,family);
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ lineOutput[iband][x]=pixelInput[iband];
+ }
+ for(int iband=0;iband<input.nrOfBand();++iband){
+ try{
+ output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in band " << iband << ", line " << y << endl;
+ }
+ }
+ progress=(1.0+y)/output.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+}
+
+void filter::Filter::dwtInverse(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+ Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
+ for(int y=0;y<input.nrOfRow();++y){
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ input.readData(lineInput[iband],GDT_Float64,y,iband);
+ vector<double> pixelInput(input.nrOfBand());
+ for(int x=0;x<input.nrOfCol();++x){
+ pixelInput=lineInput.selectCol(x);
+ dwtInverse(pixelInput,wavelet_type,family);
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ lineOutput[iband][x]=pixelInput[iband];
+ }
+ for(int iband=0;iband<input.nrOfBand();++iband){
+ try{
+ output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in band " << iband << ", line " << y << endl;
+ }
+ }
+ progress=(1.0+y)/output.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+}
+
+void filter::Filter::dwtCut(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+ Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
+ for(int y=0;y<input.nrOfRow();++y){
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ input.readData(lineInput[iband],GDT_Float64,y,iband);
+ vector<double> pixelInput(input.nrOfBand());
+ for(int x=0;x<input.nrOfCol();++x){
+ pixelInput=lineInput.selectCol(x);
+ dwtCut(pixelInput,wavelet_type,family,cut);
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ lineOutput[iband][x]=pixelInput[iband];
+ }
+ for(int iband=0;iband<input.nrOfBand();++iband){
+ try{
+ output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in band " << iband << ", line " << y << endl;
+ }
+ }
+ progress=(1.0+y)/output.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+}
+
void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family){
+ int origsize=data.size();
//make sure data size if power of 2
while(data.size()&(data.size()-1))
data.push_back(data.back());
+
int nsize=data.size();
gsl_wavelet *w;
gsl_wavelet_workspace *work;
@@ -50,19 +145,54 @@ void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wa
w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
work=gsl_wavelet_workspace_alloc(nsize);
gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
+ data.erase(data.begin()+origsize,data.end());
+ gsl_wavelet_free (w);
+ gsl_wavelet_workspace_free (work);
}
void filter::Filter::dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family){
+ int origsize=data.size();
//make sure data size if power of 2
while(data.size()&(data.size()-1))
data.push_back(data.back());
int nsize=data.size();
+ gsl_wavelet *w;
+ gsl_wavelet_workspace *work;
assert(nsize);
+ w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
+ work=gsl_wavelet_workspace_alloc(nsize);
+ gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
+ data.erase(data.begin()+origsize,data.end());
+ gsl_wavelet_free (w);
+ gsl_wavelet_workspace_free (work);
+}
+
+void filter::Filter::dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut){
+ int origsize=data.size();
+ //make sure data size if power of 2
+ while(data.size()&(data.size()-1))
+ data.push_back(data.back());
+ int nsize=data.size();
gsl_wavelet *w;
gsl_wavelet_workspace *work;
+ assert(nsize);
w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
work=gsl_wavelet_workspace_alloc(nsize);
+ gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
+ std::vector<double> abscoeff(data.size());
+ size_t* p=new size_t[data.size()];
+ for(int index=0;index<data.size();++index){
+ abscoeff[index]=fabs(data[index]);
+ }
+ int nc=(100-cut)/100.0*nsize;
+ gsl_sort_index(p,&(abscoeff[0]),1,nsize);
+ for(int i=0;(i+nc)<nsize;i++)
+ data[p[i]]=0;
gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
+ data.erase(data.begin()+origsize,data.end());
+ delete[] p;
+ gsl_wavelet_free (w);
+ gsl_wavelet_workspace_free (work);
}
void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, int offset)
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 73bf1ce..ecca488 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -23,6 +23,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <vector>
#include <iostream>
extern "C" {
+#include <gsl/gsl_sort.h>
#include <gsl/gsl_wavelet.h>
}
#include "StatFactory.h"
@@ -32,7 +33,7 @@ extern "C" {
namespace filter
{
- enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwtForward=26, dwtInverse=27, dwtQuantize=28};
+ enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwt=26, dwti=27, dwt_cut=28};
class Filter
{
@@ -67,16 +68,20 @@ public:
template<class T> void applyFwhm(const std::vector<double> &wavelengthIn, const std::vector<T>& input, const std::vector<double> &wavelengthOut, const std::vector<double> &fwhm, const std::string& interpolationType, std::vector<T>& output, bool verbose=false);
template<class T> void applyFwhm(const std::vector<double> &wavelengthIn, const Vector2d<T>& input, const std::vector<double> &wavelengthOut, const std::vector<double> &fwhm, const std::string& interpolationType, Vector2d<T>& output, int down=1, bool verbose=false);
+ void dwtForward(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
+ void dwtInverse(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
+ void dwtCut(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut);
void dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family);
void dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family);
+ void dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut);
private:
static void initFilterMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
//initialize Map
- m_filterMap["dwtForward"]=filter::dwtForward;
- m_filterMap["dwtInverse"]=filter::dwtInverse;
- m_filterMap["dwtQuantize"]=filter::dwtQuantize;
+ m_filterMap["dwt"]=filter::dwt;
+ m_filterMap["dwti"]=filter::dwti;
+ m_filterMap["dwt_cut"]=filter::dwt_cut;
m_filterMap["stdev"]=filter::stdev;
m_filterMap["var"]=filter::var;
m_filterMap["min"]=filter::min;
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index 91d2695..26faaad 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -1113,12 +1113,12 @@ void filter2d::Filter2d::dwtInverse(const ImgReaderGdal& input, ImgWriterGdal& o
}
}
-void filter2d::Filter2d::dwtQuantize(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double quantize, bool verbose){
+void filter2d::Filter2d::dwtCut(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut, bool verbose){
Vector2d<float> theBuffer;
for(int iband=0;iband<input.nrOfBand();++iband){
input.readDataBlock(theBuffer, GDT_Float32, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
std::cout << "filtering band " << iband << std::endl << std::flush;
- dwtQuantize(theBuffer, wavelet_type, family, quantize);
+ dwtCut(theBuffer, wavelet_type, family, cut);
output.writeDataBlock(theBuffer,GDT_Float32,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
}
}
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 1718fda..066c546 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -53,7 +53,7 @@ extern "C" {
namespace filter2d
{
- enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, mrf=26, dwtForward=27, dwtInverse=28, dwtQuantize=29, scramble=30, shift=31, linearfeature=32};
+ enum FILTER_TYPE { median=100, var=101 , min=102, max=103, sum=104, mean=105, minmax=106, dilate=107, erode=108, close=109, open=110, homog=111, sobelx=112, sobely=113, sobelxy=114, sobelyx=115, smooth=116, density=117, majority=118, mixed=119, threshold=120, ismin=121, ismax=122, heterog=123, order=124, stdev=125, mrf=126, dwt=127, dwti=128, dwt_cut=129, scramble=130, shift=131, linearfeature=132, smoothnodata=133, dwtz=134, dwtzi=135,dwtz_cut=136};
enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };//bicubic not supported yet...
@@ -96,10 +96,10 @@ public:
template<class T1, class T2> void smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dimX, int dimY);
void dwtForward(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
void dwtInverse(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
- void dwtQuantize(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double quantize, bool verbose=false);
+ void dwtCut(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut, bool verbose=false);
template<class T> void dwtForward(Vector2d<T>& data, const std::string& wavelet_type, int family);
template<class T> void dwtInverse(Vector2d<T>& data, const std::string& wavelet_type, int family);
- template<class T> void dwtQuantize(Vector2d<T>& data, const std::string& wavelet_type, int family, double quantize);
+ template<class T> void dwtCut(Vector2d<T>& data, const std::string& wavelet_type, int family, double cut);
void majorVoting(const std::string& inputFilename, const std::string& outputFilename,int dim=0,const std::vector<int> &prior=std::vector<int>());
/* void homogeneousSpatial(const std::string& inputFilename, const std::string& outputFilename, int dim, bool disc=false, int noValue=0); */
void doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=2, bool disc=false);
@@ -150,12 +150,15 @@ private:
m_filterMap["order"]=filter2d::order;
m_filterMap["stdev"]=filter2d::stdev;
m_filterMap["mrf"]=filter2d::mrf;
- m_filterMap["dwtForward"]=filter2d::dwtForward;
- m_filterMap["dwtInverse"]=filter2d::dwtInverse;
- m_filterMap["dwtQuantize"]=filter2d::dwtQuantize;
+ m_filterMap["dwt"]=filter2d::dwt;
+ m_filterMap["dwti"]=filter2d::dwti;
+ m_filterMap["dwt_cut"]=filter2d::dwt_cut;
m_filterMap["scramble"]=filter2d::scramble;
m_filterMap["shift"]=filter2d::shift;
m_filterMap["linearfeature"]=filter2d::linearfeature;
+ m_filterMap["dwtz"]=filter2d::dwtz;
+ m_filterMap["dwtzi"]=filter2d::dwtzi;
+ m_filterMap["dwtz_cut"]=filter2d::dwtz_cut;
}
Vector2d<double> m_taps;
@@ -781,11 +784,11 @@ template<class T> void Filter2d::dwtForward(Vector2d<T>& theBuffer, const std::s
double progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
- //make sure data size if power of 2
int nRow=theBuffer.size();
assert(nRow);
int nCol=theBuffer[0].size();
assert(nCol);
+ //make sure data size if power of 2
while(theBuffer.size()&(theBuffer.size()-1))
theBuffer.push_back(theBuffer.back());
for(int irow=0;irow<theBuffer.size();++irow)
@@ -827,11 +830,11 @@ template<class T> void Filter2d::dwtInverse(Vector2d<T>& theBuffer, const std::s
double progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
- //make sure data size if power of 2
int nRow=theBuffer.size();
assert(nRow);
int nCol=theBuffer[0].size();
assert(nCol);
+ //make sure data size if power of 2
while(theBuffer.size()&(theBuffer.size()-1))
theBuffer.push_back(theBuffer.back());
for(int irow=0;irow<theBuffer.size();++irow)
@@ -866,18 +869,18 @@ template<class T> void Filter2d::dwtInverse(Vector2d<T>& theBuffer, const std::s
gsl_wavelet_workspace_free (work);
}
-template<class T> void Filter2d::dwtQuantize(Vector2d<T>& theBuffer, const std::string& wavelet_type, int family, double quantize){
+template<class T> void Filter2d::dwtCut(Vector2d<T>& theBuffer, const std::string& wavelet_type, int family, double cut){
const char* pszMessage;
void* pProgressArg=NULL;
GDALProgressFunc pfnProgress=GDALTermProgress;
double progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
- //make sure data size if power of 2
int nRow=theBuffer.size();
assert(nRow);
int nCol=theBuffer[0].size();
assert(nCol);
+ //make sure data size if power of 2
while(theBuffer.size()&(theBuffer.size()-1))
theBuffer.push_back(theBuffer.back());
for(int irow=0;irow<theBuffer.size();++irow)
@@ -904,18 +907,12 @@ template<class T> void Filter2d::dwtQuantize(Vector2d<T>& theBuffer, const std::
for(int icol=0;icol<theBuffer[0].size();++icol){
int index=irow*theBuffer[0].size()+icol;
abscoeff[index]=fabs(data[index]);
- if(quantize<0){//absolute threshold
- if(abscoeff[index]<-quantize)
- data[index]=0;
- }
}
}
- if(quantize>0){//percentual threshold
- int nc=quantize/100.0*nsize;
- gsl_sort_index(p,abscoeff,1,nsize);
- for(int i=0;(i+nc)<nsize;i++)
- data[p[i]]=0;
- }
+ int nc=(100-cut)/100.0*nsize;
+ gsl_sort_index(p,abscoeff,1,nsize);
+ for(int i=0;(i+nc)<nsize;i++)
+ data[p[i]]=0;
gsl_wavelet2d_nstransform_inverse (w, data, theBuffer.size(), theBuffer[0].size(),theBuffer[0].size(), work);
for(int irow=0;irow<theBuffer.size();++irow){
for(int icol=0;icol<theBuffer[irow].size();++icol){
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index fa1a94c..2ae4261 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -41,7 +41,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> output_opt("o", "output", "Output image file");
Optionpk<bool> disc_opt("c", "circular", "circular disc kernel for dilation and erosion", false);
Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation (North=0, East=90, South=180, West=270).");
- Optionpk<std::string> method_opt("f", "filter", "filter function (median,var,min,max,sum,mean,minmax,dilate,erode,close,open,spatially homogeneous (central pixel must be identical to all other pixels within window),SobelX edge detection in X,SobelY edge detection in Y,SobelXY,SobelYX,smooth,density,majority voting (only for classes),forest aggregation (mixed),smooth no data (mask) values,threshold local filtering,ismin,ismax,heterogeneous (central pixel must be different than all other [...]
+ Optionpk<std::string> method_opt("f", "filter", "filter function (median,var,min,max,sum,mean,minmax,dilate,erode,close,open,spatially homogeneous (central pixel must be identical to all other pixels within window),SobelX edge detection in X,Sobely edge detection in Y,Sobelxy,Sobelyx,smooth,density,majority voting (only for classes),forest aggregation (mixed),smooth no data (mask) values,threshold local filtering,ismin,ismax,heterogeneous (central pixel must be different than all other [...]
Optionpk<std::string> resample_opt("r", "resampling-method", "Resampling method for shifting operation (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
Optionpk<int> dimX_opt("dx", "dx", "filter kernel size in x, better use odd value to avoid image shift", 3);
Optionpk<int> dimY_opt("dy", "dy", "filter kernel size in y, better use odd value to avoid image shift", 3);
@@ -49,7 +49,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
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);
Optionpk<short> class_opt("class", "class", "class value(s) to use for density, erosion, dilation, openening and closing, thresholding");
- Optionpk<double> threshold_opt("t", "threshold", "threshold value(s) to use for threshold filter (one for each class), or quantization for dwtQuantize, or sigma for shift", 0);
+ 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 sigma for shift", 0);
Optionpk<short> mask_opt("m", "mask", "mask value(s) ");
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");
Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering");
@@ -57,7 +57,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing spectral response functions (two columns: wavelength response)");
Optionpk<double> wavelengthIn_opt("win", "wavelengthIn", "list of wavelengths in input spectrum (-win band1 -win band2 ...)");
Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-wout band1 -wout band2 ...)");
- 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");
+ 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",1);
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","");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid). Use none to ommit color table");
@@ -563,16 +563,23 @@ int main(int argc,char **argv) {
filter2d.smoothNoData(input,output,dimX_opt[0],dimY_opt[0]);
break;
}
- case(filter2d::dwtForward):
+ case(filter2d::dwtz):
+ filter1d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
+ break;
+ case(filter2d::dwtzi):
+ filter1d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
+ break;
+ case(filter2d::dwtz_cut):
+ filter1d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
+ break;
+ case(filter2d::dwt):
filter2d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
break;
- case(filter2d::dwtInverse):
+ case(filter2d::dwti):
filter2d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
break;
- case(filter2d::dwtQuantize):
- if(verbose_opt[0])
- std::cout << "Quantization filtering" << std::endl;
- filter2d.dwtQuantize(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
+ case(filter2d::dwt_cut):
+ filter2d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
break;
case(filter2d::threshold):
filter2d.setThresholds(threshold_opt);
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
index 2052da2..6ef0a82 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -40,10 +40,10 @@ int main(int argc,char **argv) {
Optionpk<std::string> input_opt("i","input","input ASCII file");
Optionpk<std::string> output_opt("o", "output", "Output ASCII file");
Optionpk<int> inputCols_opt("ic", "inputCols", "input columns (e.g., for three dimensional input data in first three columns use: -ic 0 -ic 1 -ic 2");
- Optionpk<std::string> method_opt("f", "filter", "filter function (to be implemented: dwtForward, dwtInverse,dwtQuantize)");
+ Optionpk<std::string> method_opt("f", "filter", "filter function (to be implemented: dwt, dwti,dwt_cut)");
Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
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);
- Optionpk<double> threshold_opt("qt", "threshold", "Quantize threshold value", 0);
+ Optionpk<double> threshold_opt("cut", "cut", "threshold to cut dwt coefficients. Use 0 to keep all.", 0);
Optionpk<int> dimZ_opt("dz", "dz", "filter kernel size in z (band or spectral dimension), must be odd (example: 3).. Set dz>0 if 1-D filter must be used in band domain");
Optionpk<double> tapZ_opt("tapz", "tapz", "taps used for spectral filtering");
Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
@@ -96,7 +96,8 @@ int main(int argc,char **argv) {
std::cout << "wavelengthIn.size(): " << wavelengthIn.size() << std::endl;
std::cout << "inputData[0].size(): " << inputData[0].size() << std::endl;
}
- assert(wavelengthIn.size()==inputData[0].size());
+ if(wavelengthIn.size())
+ assert(wavelengthIn.size()==inputData[0].size());
asciiReader.close();
filter::Filter filter1d;
if(fwhm_opt.size()){
@@ -180,7 +181,7 @@ int main(int argc,char **argv) {
tapZ_opt[itap]=1.0/dimZ_opt[0];
}
filter1d.setTaps(tapZ_opt);
- case(filter::dwtQuantize)://deliberate fall through
+ case(filter::dwt_cut)://deliberate fall through
wavelengthOut=wavelengthIn;
break;
}
@@ -189,38 +190,15 @@ int main(int argc,char **argv) {
case(filter::smooth):
filter1d.doit(inputData[icol],filteredData[icol]);
break;
- case(filter::dwtForward):
+ case(filter::dwt):
filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
break;
- case(filter::dwtInverse):
+ case(filter::dwti):
filter1d.dwtInverse(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
break;
- case(filter::dwtQuantize):{
- int origSize=filteredData[icol].size();
- filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
- std::vector<double> abscoeff(filteredData[icol].size());
- for(int iband=0;iband<filteredData[icol].size();++iband){
- abscoeff[iband]=fabs(filteredData[icol][iband]);
- if(threshold_opt[0]<0){//absolute threshold
- if(abscoeff[iband]<-threshold_opt[0])
- filteredData[icol][iband]=0;
- }
- // if(fabs(filteredData[icol][iband])<threshold_opt[0])
- // filteredData[icol][iband]=0;
- }
- if(threshold_opt[0]>0){//percentual threshold
- int nsize=abscoeff.size();
- size_t* p=new size_t[nsize];
- int nc=threshold_opt[0]/100.0*nsize;
- gsl_sort_index(p,&(abscoeff[0]),1,nsize);
- for(int i=0;(i+nc)<nsize;i++)
- filteredData[icol][p[i]]=0;
- }
- filter1d.dwtInverse(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
- //remove extended samples
- filteredData[icol].erase(filteredData[icol].begin()+origSize,filteredData[icol].end());
+ case(filter::dwt_cut):
+ filter1d.dwtCut(filteredData[icol],wavelet_type_opt[0],family_opt[0],threshold_opt[0]);
break;
- }
default:
if(verbose_opt[0])
std::cout << "method to be implemented" << std::endl;
@@ -259,16 +237,14 @@ int main(int argc,char **argv) {
int nband=0;
if(method_opt.size()){
switch(filter::Filter::getFilterType(method_opt[0])){
- case(filter::dwtForward):
+ case(filter::dwt):
nband=filteredData[0].size();
break;
- case(filter::dwtInverse):
+ case(filter::dwti):
nband=filteredData[0].size();
break;
- case(filter::dwtQuantize):
- nband=wavelengthOut.size();
- assert(wavelengthOut.size()==nband);
- assert(filteredData[0].size()==nband);
+ case(filter::dwt_cut):
+ nband=filteredData[0].size();
break;
default:
nband=wavelengthOut.size();
@@ -285,13 +261,13 @@ int main(int argc,char **argv) {
if(!output_opt.empty()){
if(wavelengthOut.size())
outputStream << wavelengthOut[iband] << " ";
- else
+ else if(wavelengthIn_opt.size())
outputStream << iband << " ";
}
else{
if(wavelengthOut.size())
std::cout << wavelengthOut[iband] << " ";
- else
+ else if(wavelengthIn_opt.size())
std::cout << iband << " ";
}
for(int icol=0;icol<inputCols_opt.size();++icol){
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list