[pktools] 165/375: DWT in pkfilter automatically in spectral domain if dz>0
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 f9b38632cd97253dde85bb6b2d1c06bf67ce7f15
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Tue Dec 24 11:59:15 2013 +0100
DWT in pkfilter automatically in spectral domain if dz>0
---
ChangeLog | 3 +++
src/algorithms/Filter.cc | 34 +++++++++++++++++++------
src/algorithms/Filter.h | 37 +++++++++++++++++++++------
src/algorithms/Filter2d.h | 5 +---
src/apps/pkfilter.cc | 65 ++++++++++++++++++++++++-----------------------
src/apps/pkfilterascii.cc | 34 ++++++++++++-------------
6 files changed, 108 insertions(+), 70 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index c18f1dd..4403d82 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -211,7 +211,10 @@ version 2.4.3
version 2.5
- ImgReaderGdal and ImgWriterGdal
re-design of geotransform
+ - Filter2d.h
+ renamed mask to nodata
- pkinfo
hist: corrected nan when min=max
- pkfilter
corrected bug in 2d wavelet transform forward and inverse
+ renamed mask to nodata
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 6a6ffbb..6644631 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -30,14 +30,22 @@ filter::Filter::Filter(void)
filter::Filter::Filter(const vector<double> &taps)
- : m_taps(taps)
{
- assert(m_taps.size()%2);
+ setTaps(taps);
}
-void filter::Filter::setTaps(const vector<double> &taps)
+void filter::Filter::setTaps(const vector<double> &taps, bool normalize)
{
- m_taps=taps;
+ m_taps.resize(taps.size());
+ double norm=0;
+ for(int itap=0;itap<taps.size();++itap)
+ norm+=taps[itap];
+ if(norm){
+ for(int itap=0;itap<taps.size();++itap)
+ m_taps[itap]=taps[itap]/norm;
+ }
+ else
+ m_taps=taps;
assert(m_taps.size()%2);
}
@@ -197,8 +205,9 @@ void filter::Filter::dwtCut(std::vector<double>& data, const std::string& wavele
gsl_wavelet_workspace_free (work);
}
-void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, int offset)
+void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, int offset, short verbose)
{
+ bool bverbose=(verbose>1)? true:false;
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
const char* pszMessage;
@@ -213,7 +222,7 @@ void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& outpu
vector<double> pixelOutput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
- morphology(pixelInput,pixelOutput,method,dim,down,offset);
+ morphology(pixelInput,pixelOutput,method,dim,down,offset,bverbose);
for(int iband=0;iband<input.nrOfBand();++iband)
lineOutput[iband][x]=pixelOutput[iband];
}
@@ -230,7 +239,16 @@ void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& outpu
}
}
-void filter::Filter::doit(const ImgReaderGdal& input, ImgWriterGdal& output, short down, int offset)
+void filter::Filter::smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim, short down, int offset)
+{
+ assert(dim>0);
+ m_taps.resize(dim);
+ for(int itap=0;itap<dim;++itap)
+ m_taps[itap]=1.0/dim;
+ filter(input,output,down,offset);
+}
+
+void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output, short down, int offset)
{
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
@@ -246,7 +264,7 @@ void filter::Filter::doit(const ImgReaderGdal& input, ImgWriterGdal& output, sho
vector<double> pixelOutput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
- doit(pixelInput,pixelOutput,down,offset);
+ filter(pixelInput,pixelOutput,down,offset);
for(int iband=0;iband<input.nrOfBand();++iband)
lineOutput[iband][x]=pixelOutput[iband];
}
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 4af1420..d843ded 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -54,14 +54,17 @@ public:
initFilterMap(m_filterMap);
return m_filterMap[filterType];
};
- void setTaps(const std::vector<double> &taps);
+ void setTaps(const std::vector<double> &taps, bool normalize=true);
void pushClass(short theClass=1){m_class.push_back(theClass);};
void pushMask(short theMask=0){m_mask.push_back(theMask);};
- template<class T> void doit(const std::vector<T>& input, std::vector<T>& output, int down=1, int offset=0);
- template<class T> void doit(T* input, int inputSize, std::vector<T>& output, int down=1, int offset=0);
- template<class T> void morphology(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim, short down=1, int offset=0, bool verbose=0);
- void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=1, int offset=0);
- void doit(const ImgReaderGdal& input, ImgWriterGdal& output, short down=1, int offset=0);
+ template<class T> void filter(const std::vector<T>& input, std::vector<T>& output, int down=1, int offset=0);
+ template<class T> void smooth(const std::vector<T>& input, std::vector<T>& output, short dim, int down=1, int offset=0);
+ template<class T> void filter(T* input, int inputSize, std::vector<T>& output, int down=1, int offset=0);
+ template<class T> void smooth(T* input, int inputSize, std::vector<T>& output, short dim, int down=1, int offset=0);
+ template<class T> void morphology(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim, short down=1, int offset=0, bool verbose=false);
+ void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=1, int offset=0, short verbose=0);
+ void filter(const ImgReaderGdal& input, ImgWriterGdal& output, short down=1, int offset=0);
+ void smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim, short down=1, int offset=0);
double getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta=1.0, bool verbose=false);
template<class T> double applySrf(const std::vector<double> &wavelengthIn, const std::vector<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, T& output, double delta=1.0, bool normalize=false, bool verbose=false);
template<class T> double applySrf(const std::vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, std::vector<T>& output, double delta=1.0, bool normalize=false, int down=1, bool transposeInput=false, bool verbose=false);
@@ -397,7 +400,16 @@ template<class T> void Filter::applyFwhm(const std::vector<double> &wavelengthIn
}
}
-template<class T> void Filter::doit(const std::vector<T>& input, std::vector<T>& output, int down, int offset)
+ template<class T> void Filter::smooth(const std::vector<T>& input, std::vector<T>& output, short dim, int down, int offset)
+{
+ assert(dim>0);
+ m_taps.resize(dim);
+ for(int itap=0;itap<dim;++itap)
+ m_taps[itap]=1.0/dim;
+ filter(input,output,down,offset);
+ }
+
+template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T>& output, int down, int offset)
{
output.resize((input.size()-offset+down-1)/down);
int i=0;
@@ -587,7 +599,16 @@ template<class T> void Filter::morphology(const std::vector<T>& input, std::vect
}
}
-template<class T> void Filter::doit(T* input, int inputSize, std::vector<T>& output, int down, int offset)
+ template<class T> void Filter::smooth(T* input, int inputSize, std::vector<T>& output, short dim, int down, int offset)
+{
+ assert(dim>0);
+ m_taps.resize(dim);
+ for(int itap=0;itap<dim;++itap)
+ m_taps[itap]=1.0/dim;
+ filter(input,output,down,offset);
+ }
+
+template<class T> void Filter::filter(T* input, int inputSize, std::vector<T>& output, int down, int offset)
{
output.resize((inputSize-offset+down-1)/down);
int i=0;
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 873290e..6024ef5 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -53,7 +53,7 @@ extern "C" {
namespace filter2d
{
- 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 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};
enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };//bicubic not supported yet...
@@ -156,9 +156,6 @@ private:
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;
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 97824ee..7eaf741 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, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, majority voting (only for classes), smoothnodata (smooth nodata values only) values, threshold local filtering [...]
+ Optionpk<std::string> method_opt("f", "filter", "filter function (median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, majority voting (only for classes), smoothnodata (smooth nodata values only) values, threshold local filtering [...]
Optionpk<std::string> resample_opt("r", "resampling-method", "Resampling method for shifting operation (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
Optionpk<double> dimX_opt("dx", "dx", "filter kernel size in x, better use odd value to avoid image shift", 3);
Optionpk<double> dimY_opt("dy", "dy", "filter kernel size in y, better use odd value to avoid image shift", 3);
@@ -261,8 +261,14 @@ int main(int argc,char **argv) {
tapfile.close();
}
else if(tapz_opt.size()){
+ if(verbose_opt[0]){
+ std::cout << "taps: ";
+ for(int itap=0;itap<tapz_opt.size();++itap)
+ std::cout<< tapz_opt[itap] << " ";
+ std::cout<< std::endl;
+ }
filter1d.setTaps(tapz_opt);
- filter1d.doit(input,output,down_opt[0]);
+ filter1d.filter(input,output,down_opt[0]);
}
else if(fwhm_opt.size()){
if(verbose_opt[0])
@@ -369,7 +375,7 @@ int main(int argc,char **argv) {
if(dimZ_opt.size()){
if(verbose_opt[0])
std::cout<< "1-D filtering: dilate" << std::endl;
- filter1d.morphology(input,output,"dilate",dimZ_opt[0]);
+ filter1d.morphology(input,output,"dilate",dimZ_opt[0],1,0,verbose_opt[0]);
}
else{
filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
@@ -548,7 +554,7 @@ int main(int argc,char **argv) {
theTaps[2][1]=0.0;
theTaps[2][2]=1.0;
filter2d.setTaps(theTaps);
- filter2d.filter(input,output,true);
+ filter2d.filter(input,output,true);//absolute and normalize
break;
}
case(filter2d::sobely):{//Sobel edge detection in Y
@@ -567,7 +573,7 @@ int main(int argc,char **argv) {
theTaps[2][1]=-2.0;
theTaps[2][2]=-1.0;
filter2d.setTaps(theTaps);
- filter2d.filter(input,output,true);
+ filter2d.filter(input,output,true);//absolute and normalize
break;
}
case(filter2d::sobelxy):{//Sobel edge detection in XY
@@ -586,7 +592,7 @@ int main(int argc,char **argv) {
theTaps[2][1]=-1.0;
theTaps[2][2]=0.0;
filter2d.setTaps(theTaps);
- filter2d.filter(input,output,true);
+ filter2d.filter(input,output,true);//absolute and normalize
break;
}
case(filter2d::sobelyx):{//Sobel edge detection in XY
@@ -605,7 +611,7 @@ int main(int argc,char **argv) {
theTaps[2][1]=-1.0;
theTaps[2][2]=-2.0;
filter2d.setTaps(theTaps);
- filter2d.filter(input,output,true);
+ filter2d.filter(input,output,true);//absolute and normalize
break;
}
case(filter2d::smooth):{//Smoothing filter
@@ -613,7 +619,14 @@ int main(int argc,char **argv) {
std::cerr << "Error: down option not supported for this filter" << std::endl;
exit(1);
}
- filter2d.smooth(input,output,dimX_opt[0],dimY_opt[0]);
+ if(dimZ_opt.size()){
+ if(verbose_opt[0])
+ std::cout<< "1-D filtering: smooth" << std::endl;
+ filter1d.smooth(input,output,dimZ_opt[0]);
+ }
+ else{
+ filter2d.smooth(input,output,dimX_opt[0],dimY_opt[0]);
+ }
break;
}
case(filter2d::smoothnodata):{//Smoothing filter
@@ -624,47 +637,35 @@ int main(int argc,char **argv) {
filter2d.smoothNoData(input,output,dimX_opt[0],dimY_opt[0]);
break;
}
- case(filter2d::dwtz):
- if(down_opt[0]!=1){
- std::cerr << "Error: down option not supported for this filter" << std::endl;
- exit(1);
- }
- filter1d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
- break;
- case(filter2d::dwtzi):
- if(down_opt[0]!=1){
- std::cerr << "Error: down option not supported for this filter" << std::endl;
- exit(1);
- }
- filter1d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
- break;
- case(filter2d::dwtz_cut):
- if(down_opt[0]!=1){
- std::cerr << "Error: down option not supported for this filter" << std::endl;
- exit(1);
- }
- filter1d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
- break;
case(filter2d::dwt):
if(down_opt[0]!=1){
std::cerr << "Error: down option not supported for this filter" << std::endl;
exit(1);
}
- filter2d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
+ if(dimZ_opt.size())
+ filter1d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
+ else
+ filter2d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
break;
case(filter2d::dwti):
if(down_opt[0]!=1){
std::cerr << "Error: down option not supported for this filter" << std::endl;
exit(1);
}
- filter2d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
+ if(dimZ_opt.size())
+ filter1d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
+ else
+ filter2d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
break;
case(filter2d::dwt_cut):
if(down_opt[0]!=1){
std::cerr << "Error: down option not supported for this filter" << std::endl;
exit(1);
}
- filter2d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
+ if(dimZ_opt.size())
+ filter1d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
+ else
+ 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 6a6980c..d560e4b 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -174,24 +174,9 @@ int main(int argc,char **argv) {
}
if(method_opt.size()){
- switch(filter::Filter::getFilterType(method_opt[0])){
- case(filter::smooth):
- assert(dimZ_opt.size());
- if(tapZ_opt.empty()){
- tapZ_opt.resize(dimZ_opt[0]);
- for(int itap=0;itap<dimZ_opt[0];++itap)
- tapZ_opt[itap]=1.0/dimZ_opt[0];
- }
- filter1d.setTaps(tapZ_opt);
- case(filter::dwt_cut)://deliberate fall through
- wavelengthOut=wavelengthIn;
- break;
- }
+ wavelengthOut=wavelengthIn;
for(int icol=0;icol<inputCols_opt.size();++icol){
switch(filter::Filter::getFilterType(method_opt[0])){
- case(filter::smooth):
- filter1d.doit(inputData[icol],filteredData[icol]);
- break;
case(filter::dwt):
filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
break;
@@ -201,9 +186,22 @@ int main(int argc,char **argv) {
case(filter::dwt_cut):
filter1d.dwtCut(filteredData[icol],wavelet_type_opt[0],family_opt[0],threshold_opt[0]);
break;
+ case(filter::smooth):
+ if(tapZ_opt.size()){
+ filter1d.setTaps(tapZ_opt);
+ filter1d.filter(inputData[icol],filteredData[icol]);
+ }
+ else{
+ assert(dimZ_opt.size());
+ filter1d.smooth(inputData[icol],filteredData[icol],dimZ_opt[0]);
+ }
+ break;
default:
- if(verbose_opt[0])
- std::cout << "method to be implemented" << std::endl;
+ assert(tapZ_opt.size());
+ filter1d.filter(inputData[icol],filteredData[icol]);
+ // if(verbose_opt[0])
+ // std::cout << "method to be implemented" << std::endl;
+ // exit(1);
break;
}
}
--
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