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 <kempe...@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
_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
