[pktools] 100/375: geotransfer problem solved in pkfilter, directional morphology improved in Filter2d

Bas Couwenberg sebastic at xs4all.nl
Wed Dec 3 21:54:03 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 b2c9b2e61fa7501cf7b935c966d62c263eb894cd
Author: user <user at osgeolive.(none)>
Date:   Fri May 3 22:50:04 2013 +0200

    geotransfer problem solved in pkfilter, directional morphology improved in Filter2d
---
 src/algorithms/Filter2d.cc | 74 ++++++++++++++++++++++------------------------
 src/algorithms/Filter2d.h  |  2 +-
 src/apps/pkeditogr.cc      | 22 ++++++++++++++
 src/apps/pkfilter.cc       | 46 ++++++++--------------------
 src/apps/pklas2img.cc      |  4 +++
 src/apps/pksetmask.cc      |  2 +-
 6 files changed, 75 insertions(+), 75 deletions(-)

diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index ee6ecde..11a2d88 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -488,7 +488,6 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
 	    }
 	    if(!masked){
               vector<short>::const_iterator vit=m_class.begin();
-              //todo: test if this works (only add occurrence if within defined classes)!
               if(!m_class.size())
                 ++occurrence[inBuffer[indexJ][indexI]];
               else{
@@ -937,7 +936,7 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
 //   output.close();
 // }
 
-void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, bool disc, double angle)
+void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, const vector<double> &angle, bool disc)
 {
   assert(dimX);
   assert(dimY);
@@ -983,7 +982,6 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
 	outBuffer[x]=currentValue;
 	vector<double> statBuffer;
 	bool currentMasked=false;
-        double rse=0;
 	for(int imask=0;imask<m_mask.size();++imask){
 	  if(currentValue==m_mask[imask]){
 	    currentMasked=true;
@@ -999,42 +997,40 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
 	      if(disc&&(i*i+j*j>(dimX/2)*(dimY/2)))
 		continue;
               bool masked=false;
-	      if(angle>=-180){
-	      	// double theta;
-	      	// if(i>0)
-	      	//   theta=atan(static_cast<double>(j)/(static_cast<double>(i)));
-	      	// else if(i<0)
-	      	//   theta=PI+atan(static_cast<double>(j)/(static_cast<double>(i)));
-	      	// else if(j>0)
-	      	//   theta=PI/2.0;
-	      	// else if(j<0)
-	      	//   theta=3.0*PI/2.0;
-                // rse=sqrt((theta-DEG2RAD(angle))*(theta-DEG2RAD(angle))/theta/theta);
-                // //test
-                if(angle<45||angle>315)
-                  if((j!=0)||(i>0))//RIGHT
-                    continue;
-                if(angle>135&&angle<225)
-                  if(j!=0||i<0)//LEFT
-                    continue;
-                if(angle>45&&angle<135)
-                  if(j<0||i!=0)//UP
-                    continue;
-                if(angle>225&&angle<315)
-                  if(j>0||i!=0)//DOWN
-                    continue;
-                if(angle>270&&angle<360)
-                  if(j>0||i>0)//LOWER RIGHT
-                    continue;
-                if(angle>0&&angle<90)
-                  if(j<0||i>0)//UPPER RIGHT
-                    continue;
-                if(angle>180&&angle<270)
-                  if(j>0||i<0)//LOWER LEFT
-                    continue;
-                if(angle>90&&angle<180)
-                  if(j<0||i<0)//UPPER LEFT
-                    continue;
+	      if(angle.size()){
+	      	double theta;
+		//use polar coordinates in radians
+	      	if(i>0){
+		  if(j<0)
+		    theta=atan(static_cast<double>(-j)/(static_cast<double>(i)));
+		  else
+		    theta=-atan(static_cast<double>(j)/(static_cast<double>(i)));
+		}
+	      	else if(i<0){
+		  if(j<0)
+		    theta=PI-atan(static_cast<double>(-j)/(static_cast<double>(-i)));
+		  else
+		    theta=PI+atan(static_cast<double>(j)/(static_cast<double>(-i)));
+		}
+	      	else if(j<0)
+	      	  theta=PI/2.0;
+	      	else if(j>0)
+	      	  theta=3.0*PI/2.0;
+		//convert to North (0), East (90), South (180), West (270) in degrees
+		theta=360-(theta/PI*180)+90;
+		if(theta<0)
+		  theta+=360;
+		while(theta>360)
+		  theta-=360;
+		bool alligned=false;
+		for(int iangle=0;iangle<angle.size();++iangle){
+		  if(sqrt((theta-angle[iangle])*(theta-angle[iangle]))<10){
+		    alligned=true;
+		    break;
+		  }
+		}
+		if(!alligned)
+		  continue;
 	      }
 	      indexI=x+i;
 	      //check if out of bounds
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 0f9d592..fb2d3ea 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -81,7 +81,7 @@ public:
   template<class T1, class T2> void doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, const std::string& method, int dimX, int dimY, short down=1, bool disc=false);
   void median(const string& inputFilename, const string& outputFilename, int dim, bool disc=false);
   void var(const string& inputFilename, const string& outputFilename, int dim, bool disc=false);
-  void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, bool disc=false, double angle=-190);
+  void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, const vector<double> &angle, bool disc=false);
   template<class T> unsigned long int morphology(const Vector2d<T>& input, Vector2d<T>& output, const std::string& method, int dimX, int dimY, bool disc=false, double hThreshold=0);
   template<class T> void shadowDsm(const Vector2d<T>& input, Vector2d<T>& output, double sza, double saa, double pixelSize, short shadowFlag=1);
   void shadowDsm(const ImgReaderGdal& input, ImgWriterGdal& output, double sza, double saa, double pixelSize, short shadowFlag=1);
diff --git a/src/apps/pkeditogr.cc b/src/apps/pkeditogr.cc
index 9fd62f2..f764d51 100644
--- a/src/apps/pkeditogr.cc
+++ b/src/apps/pkeditogr.cc
@@ -29,6 +29,8 @@ int main(int argc, char *argv[])
 {
   Optionpk<string> input_opt("i", "input", "Input image");
   Optionpk<string> output_opt("o", "output", "Output mask file");
+  Optionpk<string> selectField_opt("select", "select", "select field (combined with like opt)");
+  Optionpk<string> like_opt("like", "like", "substring(s) to be found in select field (if multiple substrings are provided, feature will be selected if one of them is found)");
   Optionpk<string> field_opt("f", "field", "output field names (number must exactly match input fields)");
   Optionpk<long int> setfeature_opt("sf", "sf", "id of feature(s) to set (start from 0)");
   Optionpk<string> setname_opt("sn", "sn", "name(s) of field(s) to set");
@@ -42,6 +44,8 @@ int main(int argc, char *argv[])
   try{
     doProcess=input_opt.retrieveOption(argc,argv);
     output_opt.retrieveOption(argc,argv);
+    selectField_opt.retrieveOption(argc,argv);
+    like_opt.retrieveOption(argc,argv);
     field_opt.retrieveOption(argc,argv);
     addname_opt.retrieveOption(argc,argv);
     addtype_opt.retrieveOption(argc,argv);
@@ -134,6 +138,24 @@ int main(int argc, char *argv[])
   OGRFeature *poFeature;
   while((poFeature = ogrReader.getLayer()->GetNextFeature()) != NULL ){
     ++ifeature;
+    bool doSelect;
+    if(like_opt.empty())
+      doSelect=true;
+    else{
+      doSelect=false;
+      int fieldIndex=poFeature->GetFieldIndex(selectField_opt[0].c_str());
+      string fieldValue=poFeature->GetFieldAsString(fieldIndex);
+      for(int ilike=0;ilike<like_opt.size();++ilike){
+	if(fieldValue.find(like_opt[ilike])!=std::string::npos){
+	  if(verbose_opt[0])
+	    std::cout << "found " << like_opt[ilike] << " in " << fieldValue << std::endl;
+	  doSelect=true;
+	  break;
+	}
+      }
+    }
+    if(!doSelect)
+      continue;
     OGRFeature *poDstFeature = NULL;
     poDstFeature=ogrWriter.createFeature();
     if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 213fbe1..13e1d6a 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -40,8 +40,8 @@ int main(int argc,char **argv) {
   Optionpk<std::string> input_opt("i","input","input image file");
   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.");
-  Optionpk<std::string> method_opt("f", "filter", "filter function (median,variance,min,max,sum,mean,minmax,dilation,erosion,closing,opening,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  [...]
+  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,variance,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  [...]
   Optionpk<int> dimX_opt("dx", "dx", "filter kernel size in x, must be odd", 3);
   Optionpk<int> dimY_opt("dy", "dy", "filter kernel size in y, must be odd", 3);
   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");
@@ -166,11 +166,6 @@ int main(int argc,char **argv) {
   else if(input.getColorTable()!=NULL)
     output.setColorTable(input.getColorTable());
 
-  if(input.isGeoRef()){
-    output.setProjection(input.getProjection());
-    output.copyGeoTransform(input);
-  }
-
   filter2d::Filter2d filter2d;
   filter::Filter filter1d;
   if(class_opt.size()){
@@ -187,8 +182,9 @@ int main(int argc,char **argv) {
     if(verbose_opt[0])
       std::cout<< std::endl;
   }
-  if(mask_opt.size()&&verbose_opt[0]){
-    std::cout<< "mask values: ";
+  if(mask_opt.size()){
+    if(verbose_opt[0])
+      std::cout<< "mask values: ";
     for(int imask=0;imask<mask_opt.size();++imask){
       if(verbose_opt[0])
         std::cout<< mask_opt[imask] << " ";
@@ -330,23 +326,17 @@ int main(int argc,char **argv) {
         filter1d.morphology(input,output,"dilate",dimZ_opt[0]);
       }
       else{
-        if(angle_opt.size())
-          filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
-        else
-          filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0]);
+	filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
       }
       break;
     case(filter2d::erode):
-      if(dimZ_opt[0]>0){
+      if(dimZ_opt.size()>0){
         if(verbose_opt[0])
           std::cout<< "1-D filtering: dilate" << std::endl;
         filter1d.morphology(input,output,"erode",dimZ_opt[0]);
       }
       else{
-        if(angle_opt.size())
-          filter2d.morphology(input,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
-        else
-          filter2d.morphology(input,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0]);
+	filter2d.morphology(input,output,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
       }
       break;
     case(filter2d::close):{//closing
@@ -359,10 +349,7 @@ int main(int argc,char **argv) {
           filter1d.morphology(input,tmpout,"dilate",dimZ_opt[0]);
         }
         else{
-          if(angle_opt.size())
-            filter2d.morphology(input,tmpout,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
-          else
-            filter2d.morphology(input,tmpout,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0]);
+	  filter2d.morphology(input,tmpout,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
         }
       }
       catch(std::string errorString){
@@ -376,10 +363,7 @@ int main(int argc,char **argv) {
         filter1d.morphology(tmpin,output,"erode",dimZ_opt[0]);
       }
       else{
-        if(angle_opt.size())
-          filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
-        else
-          filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0]);
+          filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
       }
       tmpin.close();
       if(remove(tmps.str().c_str( )) !=0){
@@ -396,10 +380,7 @@ int main(int argc,char **argv) {
         filter1d.morphology(input,tmpout,"erode",dimZ_opt[0]);
       }
       else{
-        if(angle_opt.size())
-          filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
-        else
-          filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0]);
+	filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
       }
       tmpout.close();
       ImgReaderGdal tmpin;
@@ -408,10 +389,7 @@ int main(int argc,char **argv) {
         filter1d.morphology(tmpin,output,"dilate",dimZ_opt[0]);
       }
       else{
-        if(angle_opt.size())
-          filter2d.morphology(tmpin,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
-        else
-          filter2d.morphology(tmpin,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0]);
+	filter2d.morphology(tmpin,output,"dilate",dimX_opt[0],dimY_opt[0],angle_opt,disc_opt[0]);
       }
       tmpin.close();
       if(remove(tmps.str().c_str( )) !=0){
diff --git a/src/apps/pklas2img.cc b/src/apps/pklas2img.cc
index fecb07a..79f73e1 100644
--- a/src/apps/pklas2img.cc
+++ b/src/apps/pklas2img.cc
@@ -286,6 +286,10 @@ int main(int argc,char **argv) {
       outputWriter.geo2image(thePoint.GetX(),thePoint.GetY(),dcol,drow);
       int icol=static_cast<int>(dcol);
       int irow=static_cast<int>(drow);
+      if(irow<0||irow>=nrow)
+	continue;
+      if(icol<0||icol>=ncol)
+	continue;
       assert(irow>=0);
       assert(irow<nrow);
       assert(icol>=0);
diff --git a/src/apps/pksetmask.cc b/src/apps/pksetmask.cc
index 9b85a7f..f4cab7f 100644
--- a/src/apps/pksetmask.cc
+++ b/src/apps/pksetmask.cc
@@ -34,7 +34,7 @@ int main(int argc, char *argv[])
   Optionpk<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> option_opt("co", "co", "options: NAME=VALUE [-co COMPRESS=LZW] [-co INTERLEAVE=BAND]");
-  Optionpk<unsigned short> invalid_opt("t", "invalid", "Mask value(s) where image is invalid. Use one value for each mask, or multiple values for a single mask.", 1);
+  Optionpk<int> invalid_opt("t", "invalid", "Mask value(s) where image is invalid. Use one value for each mask, or multiple values for a single mask.", 1);
   Optionpk<char> operator_opt("p", "operator", "Operator: < = > !. Use operator for each invalid option", '=');
   Optionpk<int> flag_opt("f", "flag", "Flag value to put in image if not valid", 0);
   Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)", "");

-- 
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