[pktools] 141/375: introduced line features in Filter2d and pkfilter
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:08 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 6fd42fae0df8d3011650013eb1ee3b9fabea5ed9
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Thu Nov 21 08:22:37 2013 +0100
introduced line features in Filter2d and pkfilter
---
src/algorithms/Filter2d.cc | 130 +++++++++++++++++++++++++++++++++++++++++++++
src/algorithms/Filter2d.h | 11 +++-
src/apps/pkfilter.cc | 43 ++++++++++++++-
3 files changed, 181 insertions(+), 3 deletions(-)
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index 17df45f..e11c021 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -1111,3 +1111,133 @@ void filter2d::Filter2d::dwtQuantize(const ImgReaderGdal& input, ImgWriterGdal&
output.writeDataBlock(theBuffer,GDT_Float32,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
}
}
+
+void filter2d::Filter2d::linearFeature(const ImgReaderGdal& input, ImgWriterGdal& output, float angle, float angleStep, float maxDistance, float eps, bool l1, bool a1, bool l2, bool a2, int band, bool verbose){
+ Vector2d<float> inputBuffer;
+ vector< Vector2d<float> > outputBuffer;
+ input.readDataBlock(inputBuffer, GDT_Float32, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, band);
+ if(maxDistance<=0)
+ maxDistance=sqrt(input.nrOfCol()*input.nrOfRow());
+ linearFeature(inputBuffer,outputBuffer,angle,angleStep,maxDistance,eps, l1, a1, l2, a2,verbose);
+ for(int iband=0;iband<outputBuffer.size();++iband)
+ output.writeDataBlock(outputBuffer[iband],GDT_Float32,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
+}
+
+void filter2d::Filter2d::linearFeature(const Vector2d<float>& input, vector< Vector2d<float> >& output, float angle, float angleStep, float maxDistance, float eps, bool l1, bool a1, bool l2, bool a2, bool verbose)
+{
+ output.clear();
+ int nband=0;
+ if(l1)
+ ++nband;
+ if(a1)
+ ++nband;
+ if(l2)
+ ++nband;
+ if(a2)
+ ++nband;
+ output.resize(nband);
+ for(int iband=0;iband<output.size();++iband)
+ output[iband].resize(input.nRows(),input.nCols());
+ if(maxDistance<=0)
+ maxDistance=sqrt(input.nRows()*input.nCols());
+ int indexI=0;
+ int indexJ=0;
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ for(int y=0;y<input.nRows();++y){
+ for(int x=0;x<input.nCols();++x){
+ float currentValue=input[y][x];
+ //find values equal to current value with some error margin
+ //todo: add distance for two opposite directions
+ float lineDistance1=0;//longest line of object
+ float lineDistance2=maxDistance;//shortest line of object
+ float lineAngle1=0;//angle to longest line (North=0)
+ float lineAngle2=0;//angle to shortest line (North=0)
+ float northAngle=0;//rotating angle
+ for(northAngle=0;northAngle<180;northAngle+=angleStep){
+ if(angle<=360&&angle>=0&&angle!=northAngle)
+ continue;
+ //test
+ if(verbose)
+ std::cout << "northAngle: " << northAngle << std::endl;
+ float currentDistance=0;
+ float theDir=0;
+ for(short side=0;side<=1;side+=1){
+ theDir=PI/2.0-DEG2RAD(northAngle)+side*PI;//in radians
+ //test
+ if(verbose)
+ std::cout << "theDir in deg: " << RAD2DEG(theDir) << std::endl;
+ if(theDir<0)
+ theDir+=2*PI;
+ //test
+ if(verbose)
+ std::cout << "theDir in deg: " << RAD2DEG(theDir) << std::endl;
+ float nextValue=currentValue;
+ for(float currentRay=1;currentRay<maxDistance;++currentRay){
+ indexI=x+currentRay*cos(theDir);
+ indexJ=y-currentRay*sin(theDir);
+ if(indexJ<0||indexJ>=input.size())
+ break;
+ if(indexI<0||indexI>=input[indexJ].size())
+ break;
+ nextValue=input[indexJ][indexI];
+ if(verbose){
+ std::cout << "x: " << x << std::endl;
+ std::cout << "y: " << y << std::endl;
+ std::cout << "currentValue: " << currentValue << std::endl;
+ std::cout << "theDir in degrees: " << RAD2DEG(theDir) << std::endl;
+ std::cout << "cos(theDir): " << cos(theDir) << std::endl;
+ std::cout << "sin(theDir): " << sin(theDir) << std::endl;
+ std::cout << "currentRay: " << currentRay << std::endl;
+ std::cout << "currentDistance: " << currentDistance << std::endl;
+ std::cout << "indexI: " << indexI << std::endl;
+ std::cout << "indexJ: " << indexJ << std::endl;
+ std::cout << "nextValue: " << nextValue << std::endl;
+ }
+ if(fabs(currentValue-nextValue)<=eps){
+ ++currentDistance;
+ //test
+ if(verbose)
+ std::cout << "currentDistance: " << currentDistance << ", continue" << std::endl;
+ }
+ else{
+ if(verbose)
+ std::cout << "currentDistance: " << currentDistance << ", break" << std::endl;
+ break;
+ }
+ }
+ }
+ if(lineDistance1<currentDistance){
+ lineDistance1=currentDistance;
+ lineAngle1=northAngle;
+ }
+ if(lineDistance2>currentDistance){
+ lineDistance2=currentDistance;
+ lineAngle2=northAngle;
+ }
+ if(verbose){
+ std::cout << "lineDistance1: " << lineDistance1 << std::endl;
+ std::cout << "lineAngle1: " << lineAngle1 << std::endl;
+ std::cout << "lineDistance2: " << lineDistance2 << std::endl;
+ std::cout << "lineAngle2: " << lineAngle2 << std::endl;
+ }
+ }
+ int iband=0;
+ if(l1)
+ output[iband++][y][x]=lineDistance1;
+ if(a1)
+ output[iband++][y][x]=lineAngle1;
+ if(l2)
+ output[iband++][y][x]=lineDistance2;
+ if(a2)
+ output[iband++][y][x]=lineAngle2;
+ assert(iband==nband);
+ }
+ progress=(1.0+y);
+ progress/=input.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+}
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index ea5d690..905f596 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -28,6 +28,10 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#define DEG2RAD(DEG) (DEG/180.0*PI)
#endif
+#ifndef RAD2DEG
+#define RAD2DEG(RAD) (RAD/PI*180)
+#endif
+
#include <assert.h>
#include <math.h>
#include <limits>
@@ -49,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};
+ 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 RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };//bicubic not supported yet...
@@ -111,7 +115,8 @@ public:
void dwt_texture(const std::string& inputFilename, const std::string& outputFilename,int dim, int scale, int down=1, int iband=0, bool verbose=false);
void shift(const ImgReaderGdal& input, ImgWriterGdal& output, int offsetX=0, int offsetY=0, double randomSigma=0, RESAMPLE resample=BILINEAR, bool verbose=false);
template<class T> void shift(const Vector2d<T>& input, Vector2d<T>& output, int offsetX=0, int offsetY=0, double randomSigma=0, RESAMPLE resample=0, bool verbose=false);
-
+ void linearFeature(const Vector2d<float>& input, vector< Vector2d<float> >& output, float angle=361, float angleStep=1, float maxDistance=0, float eps=0, bool l1=true, bool a1=true, bool l2=true, bool a2=true, bool verbose=false);
+ void linearFeature(const ImgReaderGdal& input, ImgWriterGdal& output, float angle=361, float angleStep=1, float maxDistance=0, float eps=0, bool l1=true, bool a1=true, bool l2=true, bool a2=true, int band=0, bool verbose=false);
private:
static void initMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
@@ -147,6 +152,7 @@ private:
m_filterMap["dwtQuantize"]=filter2d::dwtQuantize;
m_filterMap["scramble"]=filter2d::scramble;
m_filterMap["shift"]=filter2d::shift;
+ m_filterMap["linearfeature"]=filter2d::linearfeature;
}
Vector2d<double> m_taps;
@@ -156,6 +162,7 @@ private:
std::vector<double> m_threshold;
};
+
template<class T1, class T2> void Filter2d::smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dim)
{
smooth(inputVector,outputVector,dim,dim);
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 50350a4..6ef1e33 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);
@@ -64,6 +64,11 @@ int main(int argc,char **argv) {
Optionpk<std::string> option_opt("co", "co", "options: NAME=VALUE [-co COMPRESS=LZW] [-co INTERLEAVE=BAND]");
Optionpk<short> down_opt("d", "down", "down sampling factor. Use value 1 for no downsampling). Use value n>1 for downsampling (aggregation)", 1);
Optionpk<string> beta_opt("beta", "beta", "ASCII file with beta for each class transition in Markov Random Field");
+ Optionpk<double> eps_opt("eps","eps", "error marging for linear feature",0);
+ Optionpk<bool> l1_opt("l1","l1", "obtain longest object length for linear feature",false);
+ Optionpk<bool> l2_opt("l2","l2", "obtain shortest object length for linear feature",false);
+ Optionpk<bool> a1_opt("a1","a1", "obtain angle found for longest object length for linear feature",false);
+ Optionpk<bool> a2_opt("a2","a2", "obtain angle found for shortest object length for linear feature",false);
Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
bool doProcess;//stop process when program was invoked with help option (-h --help)
@@ -91,6 +96,11 @@ int main(int argc,char **argv) {
wavelengthOut_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
beta_opt.retrieveOption(argc,argv);
+ eps_opt.retrieveOption(argc,argv);
+ l1_opt.retrieveOption(argc,argv);
+ l2_opt.retrieveOption(argc,argv);
+ a1_opt.retrieveOption(argc,argv);
+ a2_opt.retrieveOption(argc,argv);
interpolationType_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
@@ -145,6 +155,20 @@ int main(int argc,char **argv) {
std::cout << "opening output image " << output_opt[0] << std::endl;
output.open(output_opt[0],(input.nrOfCol()+down_opt[0]-1)/down_opt[0],(input.nrOfRow()+down_opt[0]-1)/down_opt[0],class_opt.size(),theType,imageType,option_opt);
}
+ else if(filter2d::Filter2d::getFilterType(method_opt[0])==filter2d::linearfeature){
+ if(verbose_opt[0])
+ std::cout << "opening output image " << output_opt[0] << std::endl;
+ int nband=0;
+ if(l1_opt[0])
+ ++nband;
+ if(a1_opt[0])
+ ++nband;
+ if(l2_opt[0])
+ ++nband;
+ if(a2_opt[0])
+ ++nband;
+ output.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),nband,theType,imageType,option_opt);
+ }
else if(fwhm_opt.size()||srf_opt.size()){
//todo: support down and offset
int nband=fwhm_opt.size()? fwhm_opt.size():srf_opt.size();
@@ -424,6 +448,23 @@ int main(int argc,char **argv) {
}
break;
}
+ case(filter2d::linearfeature):{
+ assert(input.nrOfBand());
+ assert(input.nrOfCol());
+ assert(input.nrOfRow());
+ float theAngle=361;
+ if(angle_opt.size())
+ theAngle=angle_opt[0];
+ if(verbose_opt[0])
+ std::cout << "using angle " << theAngle << std::endl;
+ try{
+ 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]);
+ }
+ catch(string errorstring){
+ cerr << errorstring << endl;
+ }
+ break;
+ }
case(filter2d::mrf):{//Markov Random Field
if(verbose_opt[0])
std::cout << "Markov Random Field filtering" << std::endl;
--
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