[pktools] 163/375: redesign of geotransform in ImgReaderGdal and ImgWriterGdal

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 17fe2c8ad3c0dbf769ef9a84ec02a1266873cf9f
Author: Pieter Kempeneers <kempenep at gmail.com>
Date:   Tue Dec 24 00:15:58 2013 +0100

    redesign of geotransform in ImgReaderGdal and ImgWriterGdal
---
 src/algorithms/Filter.cc          |   2 +
 src/algorithms/Filter.h           |  22 ++--
 src/algorithms/Filter2d.cc        | 251 ++++++++++++++++++++------------------
 src/algorithms/Filter2d.h         |  69 +++++------
 src/apps/pkascii2img.cc           |  10 +-
 src/apps/pkcrop.cc                |  21 +++-
 src/apps/pkdsm2shadow.cc          |   8 +-
 src/apps/pkdumpimg.cc             |   5 +-
 src/apps/pkenhance.cc             |   2 +
 src/apps/pkfilter.cc              | 125 ++++++++++++++++---
 src/apps/pkfilterascii.cc         |   2 +
 src/apps/pkgetmask.cc             |   6 +-
 src/apps/pkinfo.cc                |   4 +-
 src/apps/pklas2img.cc             |  19 ++-
 src/apps/pkmosaic.cc              |  10 +-
 src/imageclasses/ImgReaderGdal.cc |  96 ++++++++-------
 src/imageclasses/ImgReaderGdal.h  |  28 ++---
 src/imageclasses/ImgWriterGdal.cc | 206 ++++++++++++++++---------------
 src/imageclasses/ImgWriterGdal.h  | 103 ++++++++--------
 19 files changed, 582 insertions(+), 407 deletions(-)

diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 40357ac..6a6ffbb 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -22,6 +22,8 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include <math.h>
 #include <iostream>
 
+using namespace std;
+
 filter::Filter::Filter(void)
 {
 }
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index ecca488..4af1420 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -480,15 +480,15 @@ template<class T> void Filter::morphology(const std::vector<T>& input, std::vect
       output[(i-offset+down-1)/down]=stat.min(statBuffer);
       break;
     default:
-      string errorString="method not supported";
+      std::string errorString="method not supported";
       throw(errorString);
       break;
     }
     if(verbose){
-      cout << "buffer: ";
+      std::cout << "buffer: ";
       for(int ibuf=0;ibuf<statBuffer.size();++ibuf)
-        cout << statBuffer[ibuf] << " ";
-      cout << "->" << output[(i-offset+down-1)/down] << endl;
+        std::cout << statBuffer[ibuf] << " ";
+      std::cout << "->" << output[(i-offset+down-1)/down] << std::endl;
     }
   }
   //main
@@ -520,15 +520,15 @@ template<class T> void Filter::morphology(const std::vector<T>& input, std::vect
       output[(i-offset+down-1)/down]=stat.min(statBuffer);
       break;
     default:
-      string errorString="method not supported";
+      std::string errorString="method not supported";
       throw(errorString);
       break;
     }
     if(verbose){
-      cout << "buffer: ";
+      std::cout << "buffer: ";
       for(int ibuf=0;ibuf<statBuffer.size();++ibuf)
-        cout << statBuffer[ibuf] << " ";
-      cout << "->" << output[(i-offset+down-1)/down] << endl;
+        std::cout << statBuffer[ibuf] << " ";
+      std::cout << "->" << output[(i-offset+down-1)/down] << std::endl;
     }
     statBuffer.clear();
   }
@@ -579,10 +579,10 @@ template<class T> void Filter::morphology(const std::vector<T>& input, std::vect
       break;
     }
     if(verbose){
-      cout << "buffer: ";
+      std::cout << "buffer: ";
       for(int ibuf=0;ibuf<statBuffer.size();++ibuf)
-        cout << statBuffer[ibuf] << " ";
-      cout << "->" << output[(i-offset+down-1)/down] << endl;
+        std::cout << statBuffer[ibuf] << " ";
+      std::cout << "->" << output[(i-offset+down-1)/down] << std::endl;
     }
   }
 }
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index 26faaad..7171c72 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -26,15 +26,21 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 // #include "imageclasses/ImgUtils.h"
 
 filter2d::Filter2d::Filter2d(void)
-  : m_noValue(0)
 {
 }
 
 filter2d::Filter2d::Filter2d(const Vector2d<double> &taps)
-  : m_taps(taps), m_noValue(0)
+  : m_taps(taps)
 {
 }
 
+int filter2d::Filter2d::pushNoDataValue(double noDataValue)
+{
+  if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
+    m_noDataValues.push_back(noDataValue);
+  return(m_noDataValues.size());
+}
+
 void filter2d::Filter2d::setTaps(const Vector2d<double> &taps)
 {
   m_taps=taps;
@@ -85,7 +91,7 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input, ImgWriterGdal& outpu
   pfnProgress(progress,pszMessage,pProgressArg);
   for(int iband=0;iband<input.nrOfBand();++iband){
     Vector2d<double> inBuffer(dimY,input.nrOfCol());
-    vector<double> outBuffer(input.nrOfCol());
+    std::vector<double> outBuffer(input.nrOfCol());
     int indexI=0;
     int indexJ=0;
     //initialize last half of inBuffer
@@ -93,8 +99,8 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input, ImgWriterGdal& outpu
       try{
         input.readData(inBuffer[indexJ],GDT_Float64,abs(j),iband);
       }
-      catch(string errorstring){
-	cerr << errorstring << "in line " << indexJ << endl;
+      catch(std::string errorstring){
+	std::cerr << errorstring << "in line " << indexJ << std::endl;
       }
       ++indexJ;
     }
@@ -110,8 +116,8 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input, ImgWriterGdal& outpu
 	  try{
             input.readData(inBuffer[inBuffer.size()-1],GDT_Float64,y+dimY/2,iband);
 	  }
-	  catch(string errorstring){
-	    cerr << errorstring << "in band " << iband << ", line " << y << endl;
+	  catch(std::string errorstring){
+	    std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
 	  }
 	}
         else{
@@ -127,8 +133,8 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input, ImgWriterGdal& outpu
         double norm=0;
         bool masked=false;
         if(noData){//only filter noData values
-          for(int imask=0;imask<m_mask.size();++imask){
-            if(inBuffer[(dimY-1)/2][x]==m_mask[imask]){
+          for(int imask=0;imask<m_noDataValues.size();++imask){
+            if(inBuffer[(dimY-1)/2][x]==m_noDataValues[imask]){
               masked=true;
               break;
             }
@@ -154,8 +160,8 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input, ImgWriterGdal& outpu
 	      indexJ=(dimY-1)/2-abs(j);
             //do not take masked values into account
             masked=false;
-	    for(int imask=0;imask<m_mask.size();++imask){
-	      if(inBuffer[indexJ][indexI]==m_mask[imask]){
+	    for(int imask=0;imask<m_noDataValues.size();++imask){
+	      if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
 		masked=true;
 		break;
 	      }
@@ -175,8 +181,8 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input, ImgWriterGdal& outpu
       try{
         output.writeData(outBuffer,GDT_Float64,y,iband);
       }
-      catch(string errorstring){
-	    cerr << errorstring << "in band " << iband << ", line " << y << endl;
+      catch(std::string errorstring){
+	    std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
       }
       progress=(1.0+y);
       progress+=(output.nrOfRow()*iband);
@@ -187,7 +193,7 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input, ImgWriterGdal& outpu
 }
 
 
-void filter2d::Filter2d::majorVoting(const string& inputFilename, const string& outputFilename,int dim,const vector<int> &prior)
+void filter2d::Filter2d::majorVoting(const std::string& inputFilename, const std::string& outputFilename,int dim,const std::vector<int> &prior)
 {
   const char* pszMessage;
   void* pProgressArg=NULL;
@@ -197,14 +203,14 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
 
   bool usePriors=true;  
   if(prior.empty()){
-    cout << "no prior information" << endl;
+    std::cout << "no prior information" << std::endl;
     usePriors=false;
   }
   else{
-    cout << "using priors ";    
+    std::cout << "using priors ";    
     for(int iclass=0;iclass<prior.size();++iclass)
-      cout << " " << static_cast<short>(prior[iclass]);
-    cout << endl;    
+      std::cout << " " << static_cast<short>(prior[iclass]);
+    std::cout << std::endl;    
   }  
 
   ImgReaderGdal input;
@@ -226,7 +232,7 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
   assert(dimY);
 
   Vector2d<double> inBuffer(dimY,input.nrOfCol());
-  vector<double> outBuffer(input.nrOfCol());
+  std::vector<double> outBuffer(input.nrOfCol());
   int indexI=0;
   int indexJ=0;
   //initialize last half of inBuffer
@@ -234,8 +240,8 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
       try{
         input.readData(inBuffer[indexJ],GDT_Float64,abs(j));
       }
-      catch(string errorstring){
-	cerr << errorstring << "in line " << indexJ << endl;
+      catch(std::string errorstring){
+	std::cerr << errorstring << "in line " << indexJ << std::endl;
       }
       ++indexJ;
     }
@@ -251,8 +257,8 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
 	try{
           input.readData(inBuffer[inBuffer.size()-1],GDT_Float64,y+dimY/2);
 	}
-	catch(string errorstring){
-	  cerr << errorstring << "in line" << y << endl;
+	catch(std::string errorstring){
+	  std::cerr << errorstring << "in line" << y << std::endl;
 	}
       }
       else{
@@ -265,7 +271,7 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
     }
     for(int x=0;x<input.nrOfCol();++x){
       outBuffer[x]=0;
-      map<int,int> occurrence;
+      std::map<int,int> occurrence;
       int centre=dimX*(dimY-1)/2+(dimX-1)/2;
       for(int j=-(dimY-1)/2;j<=dimY/2;++j){
         for(int i=-(dimX-1)/2;i<=dimX/2;++i){
@@ -297,8 +303,8 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
 	    ++occurrence[inBuffer[indexJ][indexI]];
 	}
       }
-      map<int,int>::const_iterator maxit=occurrence.begin();
-      for(map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
+      std::map<int,int>::const_iterator maxit=occurrence.begin();
+      for(std::map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
 	if(mit->second>maxit->second)
 	  maxit=mit;
       }
@@ -311,8 +317,8 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
     try{
       output.writeData(outBuffer,GDT_Float64,y);
     }
-    catch(string errorstring){
-      cerr << errorstring << "in line" << y << endl;
+    catch(std::string errorstring){
+      std::cerr << errorstring << "in line" << y << std::endl;
     }
     progress=(1.0+y)/output.nrOfRow();
     pfnProgress(progress,pszMessage,pProgressArg);
@@ -321,7 +327,7 @@ void filter2d::Filter2d::majorVoting(const string& inputFilename, const string&
   output.close();
 }
 
-void filter2d::Filter2d::median(const string& inputFilename, const string& outputFilename,int dim, bool disc)
+void filter2d::Filter2d::median(const std::string& inputFilename, const std::string& outputFilename,int dim, bool disc)
 {
   ImgReaderGdal input;
   ImgWriterGdal output;
@@ -330,7 +336,7 @@ void filter2d::Filter2d::median(const string& inputFilename, const string& outpu
   doit(input,output,"median",dim,disc);
 }
 
-void filter2d::Filter2d::var(const string& inputFilename, const string& outputFilename,int dim, bool disc)
+void filter2d::Filter2d::var(const std::string& inputFilename, const std::string& outputFilename,int dim, bool disc)
 {
   ImgReaderGdal input;
   ImgWriterGdal output;
@@ -358,7 +364,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
   statfactory::StatFactory stat;
   for(int iband=0;iband<input.nrOfBand();++iband){
     Vector2d<double> inBuffer(dimY,input.nrOfCol());
-    vector<double> outBuffer((input.nrOfCol()+down-1)/down);
+    std::vector<double> outBuffer((input.nrOfCol()+down-1)/down);
     int indexI=0;
     int indexJ=0;
     //initialize last half of inBuffer
@@ -366,8 +372,8 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
       try{
         input.readData(inBuffer[indexJ],GDT_Float64,abs(j),iband);
       }
-      catch(string errorstring){
-	cerr << errorstring << "in line " << indexJ << endl;
+      catch(std::string errorstring){
+	std::cerr << errorstring << "in line " << indexJ << std::endl;
       }
       ++indexJ;
     }
@@ -382,8 +388,8 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
 	  try{
             input.readData(inBuffer[inBuffer.size()-1],GDT_Float64,y+dimY/2,iband);
 	  }
-	  catch(string errorstring){
-	    cerr << errorstring << "in band " << iband << ", line " << y << endl;
+	  catch(std::string errorstring){
+	    std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
 	  }
 	}
         else{
@@ -400,8 +406,8 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
         if((x+1+down/2)%down)
           continue;
 	outBuffer[x/down]=0;
-	vector<double> windowBuffer;
-	map<int,int> occurrence;
+	std::vector<double> windowBuffer;
+	std::map<int,int> occurrence;
         int centre=dimX*(dimY-1)/2+(dimX-1)/2;
 	for(int j=-(dimY-1)/2;j<=dimY/2;++j){
 	  for(int i=-(dimX-1)/2;i<=dimX/2;++i){
@@ -421,14 +427,14 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
 	    else
 	      indexJ=(dimY-1)/2+j;
 	    bool masked=false;
-	    for(int imask=0;imask<m_mask.size();++imask){
-	      if(inBuffer[indexJ][indexI]==m_mask[imask]){
+	    for(int imask=0;imask<m_noDataValues.size();++imask){
+	      if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
 		masked=true;
 		break;
 	      }
 	    }
 	    if(!masked){
-              vector<short>::const_iterator vit=m_class.begin();
+              std::vector<short>::const_iterator vit=m_class.begin();
               if(!m_class.size())
                 ++occurrence[inBuffer[indexJ][indexI]];
               else{
@@ -444,50 +450,50 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
         switch(getFilterType(method)){
         case(filter2d::median):
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
             outBuffer[x/down]=stat.median(windowBuffer);
           break;
         case(filter2d::var):{
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
             outBuffer[x/down]=stat.var(windowBuffer);
           break;
         }
         case(filter2d::stdev):{
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
             outBuffer[x/down]=sqrt(stat.var(windowBuffer));
           break;
         }
         case(filter2d::mean):{
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
             outBuffer[x/down]=stat.mean(windowBuffer);
           break;
         }
         case(filter2d::min):{
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
            outBuffer[x/down]=stat.min(windowBuffer);
           break;
         }
         case(filter2d::ismin):{
            if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
             outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[centre])? 1:0;
           break;
         }
-        case(filter2d::minmax):{
+        case(filter2d::minmax):{//is the same as homog?
           double min=0;
           double max=0;
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else{
             stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
             if(min!=max)
@@ -499,21 +505,21 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
         }
         case(filter2d::max):{
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
             outBuffer[x/down]=stat.max(windowBuffer);
           break;
         }
         case(filter2d::ismax):{
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else
             outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[centre])? 1:0;
           break;
         }
         case(filter2d::order):{
           if(windowBuffer.empty())
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           else{
             double lbound=0;
             double ubound=dimX*dimY;
@@ -529,38 +535,49 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
           break;
         }
         case(filter2d::homog):
-	  if(occurrence.size()==1)//all values in window must be the same
+	  if(occurrence.size()==1)//all values in window are the same
 	    outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
-	  else//favorize original value in case of ties
-	    outBuffer[x/down]=m_noValue;
+	  else
+	    outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           break;
         case(filter2d::heterog):{
-          for(vector<double>::const_iterator wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
-            if(wit==windowBuffer.begin()+windowBuffer.size()/2)
-              continue;
-            else if(*wit!=inBuffer[(dimY-1)/2][x])
-              outBuffer[x/down]=1;
-            else if(*wit==inBuffer[(dimY-1)/2][x]){//todo:wit mag niet central pixel zijn
-              outBuffer[x/down]=m_noValue;
-              break;
-            }
-          }
-          break;
+	  if(occurrence.size()==windowBuffer.size())
+	    outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
+	  else	    
+	    outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
+	  // if(occurrence.size()==1)//all values in window are the same
+	  //   outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
+	  // else
+	  //   outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
+          // break;
+          // for(std::vector<double>::const_iterator wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
+          //   if(wit==windowBuffer.begin()+windowBuffer.size()/2)
+          //     continue;
+          //   else if(*wit!=inBuffer[(dimY-1)/2][x]){
+          //     outBuffer[x/down]=1;
+	  //     break;
+	  //   }
+          //   else if(*wit==inBuffer[(dimY-1)/2][x]){//todo:wit mag niet central pixel zijn
+          //     outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
+          //     break;
+          //   }
+          // }
+          // break;
         }
         case(filter2d::density):{
 	  if(windowBuffer.size()){
-	    vector<short>::const_iterator vit=m_class.begin();
+	    std::vector<short>::const_iterator vit=m_class.begin();
 	    while(vit!=m_class.end())
 	      outBuffer[x/down]+=100.0*occurrence[*(vit++)]/windowBuffer.size();
 	  }
 	  else
-	    outBuffer[x/down]=m_noValue;
+	    outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           break;
 	}
         case(filter2d::majority):{
 	  if(occurrence.size()){
-            map<int,int>::const_iterator maxit=occurrence.begin();
-            for(map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
+            std::map<int,int>::const_iterator maxit=occurrence.begin();
+            for(std::map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
               if(mit->second>maxit->second)
                 maxit=mit;
             }
@@ -570,7 +587,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
               outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
 	  }
 	  else
-	    outBuffer[x/down]=m_noValue;
+	    outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           break;
         }
         case(filter2d::threshold):{
@@ -583,7 +600,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
             }
           }
           else
-	    outBuffer[x/down]=m_noValue;
+	    outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           break;
         }
         case(filter2d::scramble):{//could be done more efficiently window by window with random shuffling entire buffer and assigning entire buffer at once to output image...
@@ -592,7 +609,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
             outBuffer[x/down]=windowBuffer[randomIndex];
           }
           else
-	    outBuffer[x/down]=m_noValue;
+	    outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
           break;
         }
         case(filter2d::mixed):{
@@ -634,8 +651,8 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
       try{
         output.writeData(outBuffer,GDT_Float64,y/down,iband);
       }
-      catch(string errorstring){
-	cerr << errorstring << "in band " << iband << ", line " << y << endl;
+      catch(std::string errorstring){
+	std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
       }
     }
   }
@@ -675,8 +692,8 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
     try{
       input.readData(inBuffer[indexJ],GDT_Int16,abs(j));
     }
-    catch(string errorstring){
-      cerr << errorstring << "in line " << indexJ << endl;
+    catch(std::string errorstring){
+      std::cerr << errorstring << "in line " << indexJ << std::endl;
     }
     ++indexJ;
   }
@@ -691,8 +708,8 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
         try{
           input.readData(inBuffer[inBuffer.size()-1],GDT_Int16,y+dimY/2);
         }
-        catch(string errorstring){
-          cerr << errorstring << "in line " << y << endl;
+        catch(std::string errorstring){
+          std::cerr << errorstring << "in line " << y << std::endl;
         }
       }
       else{
@@ -708,12 +725,12 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
     for(int x=0;x<input.nrOfCol();++x){
       if((x+1+down/2)%down)
         continue;
-      vector<short> potential(m_class.size());
+      std::vector<short> potential(m_class.size());
       for(int iclass=0;iclass<m_class.size();++iclass){
         potential[iclass]=0;
         outBuffer[iclass][x/down]=0;
       }
-      vector<double> windowBuffer;
+      std::vector<double> windowBuffer;
       int centre=dimX*(dimY-1)/2+(dimX-1)/2;
       for(int j=-(dimY-1)/2;j<=dimY/2;++j){
         for(int i=-(dimX-1)/2;i<=dimX/2;++i){
@@ -734,8 +751,8 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
           else
             indexJ=(dimY-1)/2+j;
           bool masked=false;
-          for(int imask=0;imask<m_mask.size();++imask){
-            if(inBuffer[indexJ][indexI]==m_mask[imask]){
+          for(int imask=0;imask<m_noDataValues.size();++imask){
+            if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
               masked=true;
               break;
             }
@@ -774,14 +791,14 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
       try{
         output.writeData(outBuffer[iclass],GDT_Float64,y/down,iclass);
       }
-      catch(string errorstring){
-        cerr << errorstring << "in class " << iclass << ", line " << y << endl;
+      catch(std::string errorstring){
+        std::cerr << errorstring << "in class " << iclass << ", line " << y << std::endl;
       }
     }
   }
 }
 
-void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output, int offsetX, int offsetY, double randomSigma, RESAMPLE resample, bool verbose)
+void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output, double offsetX, double offsetY, double randomSigma, RESAMPLE resample, bool verbose)
 {
   assert(input.nrOfCol()==output.nrOfCol());
   assert(input.nrOfRow()==output.nrOfRow());
@@ -802,12 +819,12 @@ void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output
 }
 
 //todo: re-implement without dependency of CImg and reg libraries
-// void filter2d::Filter2d::dwt_texture(const string& inputFilename, const string& outputFilename,int dim, int scale, int down, int iband, bool verbose)
+// void filter2d::Filter2d::dwt_texture(const std::string& inputFilename, const std::string& outputFilename,int dim, int scale, int down, int iband, bool verbose)
 // {
 //   ImgReaderGdal input;
 //   ImgWriterGdal output;
 //   if(verbose)
-//     cout << "opening file " << inputFilename << endl;
+//     std::cout << "opening file " << inputFilename << std::endl;
 //   input.open(inputFilename);
 //   double magicX=1,magicY=1;
 //   output.open(outputFilename,(input.nrOfCol()+down-1)/down,(input.nrOfRow()+down-1)/down,scale*3,GDT_Float32,input.getImageType());
@@ -816,7 +833,7 @@ void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output
 //     output.copyGeoTransform(input);
 //   }
 //   if(verbose)
-//     cout << "Dimension texture (row x col x band) = " << (input.nrOfCol()+down-1)/down << " x " << (input.nrOfRow()+down-1)/down << " x " << scale*3 << endl;
+//     std::cout << "Dimension texture (row x col x band) = " << (input.nrOfCol()+down-1)/down << " x " << (input.nrOfRow()+down-1)/down << " x " << scale*3 << std::endl;
 //   assert(dim%2);
 //   int dimX=dim;
 //   int dimY=dim;
@@ -828,12 +845,12 @@ void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output
 //   for(int j=-dimY/2;j<(dimY+1)/2;++j){
 //     try{
 //       if(verbose)
-// 	cout << "reading input line " << abs(j) << endl;
+// 	cout << "reading input line " << abs(j) << std::endl;
 //       input.readData(inBuffer[indexJ],GDT_Float32,abs(j),iband);
 //       ++indexJ;
 //     }
-//     catch(string errorstring){
-//       cerr << errorstring << "in band " << iband << ", line " << indexJ << endl;
+//     catch(std::string errorstring){
+//       std::cerr << errorstring << "in band " << iband << ", line " << indexJ << std::endl;
 //     }
 //   }
 //   const char* pszMessage;
@@ -843,7 +860,7 @@ void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output
 //   pfnProgress(progress,pszMessage,pProgressArg);
 //   for(int y=0;y<input.nrOfRow();y+=down){
 //     if(verbose)
-//       cout << "calculating line " << y/down << endl;
+//       std::cout << "calculating line " << y/down << std::endl;
 //     if(y){//inBuffer already initialized for y=0
 //       //erase first line from inBuffer
 //       inBuffer.erase(inBuffer.begin());
@@ -853,11 +870,11 @@ void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output
 // 	inBuffer.push_back(inBuffer.back());
 // 	try{
 // 	  if(verbose)
-// 	    cout << "reading input line " << y+dimY/2 << endl;
+// 	    std::cout << "reading input line " << y+dimY/2 << std::endl;
 //           input.readData(inBuffer[inBuffer.size()-1],GDT_Float32,y+dimY/2,iband);
 // 	}
-// 	catch(string errorstring){
-// 	  cerr << errorstring << "in band " << iband << ", line " << y << endl;
+// 	catch(std::string errorstring){
+// 	  std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
 // 	}
 //       }
 //     }
@@ -894,12 +911,12 @@ void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output
 //     //write outBuffer to file
 //     try{
 //       if(verbose)
-//         cout << "writing line " << y/down << endl;
+//         std::cout << "writing line " << y/down << std::endl;
 //       for(int v=0;v<scale*3;++v)
 //         output.writeData(outBuffer[v],GDT_Float32,y/down,v);
 //     }
-//     catch(string errorstring){
-//       cerr << errorstring << "in band " << iband << ", line " << y << endl;
+//     catch(std::string errorstring){
+//       std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
 //     }
 //     progress=(1.0+y)/output.nrOfRow();
 //     pfnProgress(progress,pszMessage,pProgressArg);
@@ -908,7 +925,7 @@ void filter2d::Filter2d::shift(const ImgReaderGdal& input, ImgWriterGdal& output
 //   output.close();
 // }
 
-void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, const vector<double> &angle, bool disc)
+void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, const std::vector<double> &angle, bool disc)
 {
   const char* pszMessage;
   void* pProgressArg=NULL;
@@ -922,7 +939,7 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
   statfactory::StatFactory stat;
   for(int iband=0;iband<input.nrOfBand();++iband){
     Vector2d<double> inBuffer(dimY,input.nrOfCol());
-    vector<double> outBuffer(input.nrOfCol());
+    std::vector<double> outBuffer(input.nrOfCol());
     int indexI=0;
     int indexJ=0;
     //initialize last half of inBuffer
@@ -931,8 +948,8 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
 	input.readData(inBuffer[indexJ],GDT_Float64,abs(j),iband);
 	++indexJ;
       }
-      catch(string errorstring){
-	cerr << errorstring << "in line " << indexJ << endl;
+      catch(std::string errorstring){
+	std::cerr << errorstring << "in line " << indexJ << std::endl;
       }
     }
     for(int y=0;y<input.nrOfRow();++y){
@@ -946,8 +963,8 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
 	  try{
             input.readData(inBuffer[inBuffer.size()-1],GDT_Float64,y+dimY/2,iband);
 	  }
-	  catch(string errorstring){
-	    cerr << errorstring << "in band " << iband << ", line " << y << endl;
+	  catch(std::string errorstring){
+	    std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
 	  }
 	}
         else{
@@ -961,11 +978,11 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
       for(int x=0;x<input.nrOfCol();++x){
         double currentValue=inBuffer[(dimY-1)/2][x];
 	outBuffer[x]=currentValue;
-	vector<double> statBuffer;
+	std::vector<double> statBuffer;
 	bool currentMasked=false;
         int centre=dimX*(dimY-1)/2+(dimX-1)/2;
-	for(int imask=0;imask<m_mask.size();++imask){
-	  if(currentValue==m_mask[imask]){
+	for(int imask=0;imask<m_noDataValues.size();++imask){
+	  if(currentValue==m_noDataValues[imask]){
 	    currentMasked=true;
 	    break;
 	  }
@@ -1027,11 +1044,11 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
               else
                 indexJ=(dimY-1)/2+j;
               //todo: introduce novalue as this: ?
-              // if(inBuffer[indexJ][indexI]==m_noValue)
+              // if(inBuffer[indexJ][indexI]==(m_noDataValues.size())? m_noDataValues[0] : 0)
               //   continue;
               bool masked=false;
-	      for(int imask=0;imask<m_mask.size();++imask){
-		if(inBuffer[indexJ][indexI]==m_mask[imask]){
+	      for(int imask=0;imask<m_noDataValues.size();++imask){
+		if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
 		  masked=true;
 		  break;
 		}
@@ -1060,8 +1077,8 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
               outBuffer[x]=stat.min(statBuffer);
               break;
             default:
-              ostringstream ess;
-              ess << "Error:  morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << endl;
+              std::ostringstream ess;
+              ess << "Error:  morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << std::endl;
               throw(ess.str());
               break;
             }
@@ -1074,8 +1091,8 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
       try{
         output.writeData(outBuffer,GDT_Float64,y,iband);
       }
-      catch(string errorstring){
-	cerr << errorstring << "in band " << iband << ", line " << y << endl;
+      catch(std::string errorstring){
+	std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
       }
       progress=(1.0+y);
       progress+=(output.nrOfRow()*iband);
@@ -1125,7 +1142,7 @@ void filter2d::Filter2d::dwtCut(const ImgReaderGdal& input, ImgWriterGdal& outpu
 
 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;
+  std::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());
@@ -1134,10 +1151,10 @@ void filter2d::Filter2d::linearFeature(const ImgReaderGdal& input, ImgWriterGdal
     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)
+void filter2d::Filter2d::linearFeature(const Vector2d<float>& input, std::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;
+  int nband=0;//linear feature
   if(l1)
     ++nband;
   if(a1)
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 066c546..873290e 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -80,9 +80,9 @@ public:
   };
 
   void setTaps(const Vector2d<double> &taps);
-  void setNoValue(double noValue=0){m_noValue=noValue;};
+  /* void setNoValue(double noValue=0){m_noValue=noValue;}; */
   void pushClass(short theClass=1){m_class.push_back(theClass);};
-  void pushMask(short theMask=0){m_mask.push_back(theMask);};
+  int pushNoDataValue(double noDataValue=0);//{m_mask.push_back(theMask);};
   void pushThreshold(double theThreshold){m_threshold.push_back(theThreshold);};
   void setThresholds(const std::vector<double>& theThresholds){m_threshold=theThresholds;};
   void setClasses(const std::vector<short>& theClasses){m_class=theClasses;};
@@ -114,9 +114,9 @@ public:
   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);
   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 shift(const ImgReaderGdal& input, ImgWriterGdal& output, double offsetX=0, double offsetY=0, double randomSigma=0, RESAMPLE resample=BILINEAR, bool verbose=false);
+  template<class T> void shift(const Vector2d<T>& input, Vector2d<T>& output, double offsetX=0, double offsetY=0, double randomSigma=0, RESAMPLE resample=0, bool verbose=false);
+  void linearFeature(const Vector2d<float>& input, std::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:
@@ -162,9 +162,10 @@ private:
   }
 
   Vector2d<double> m_taps;
-  double m_noValue;
+  /* double m_noValue; */
   std::vector<short> m_class;
-  std::vector<short> m_mask;
+  /* std::vector<short> m_mask; */
+  std::vector<double> m_noDataValues;
   std::vector<double> m_threshold;
 };
 
@@ -303,8 +304,8 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
           else
             indexJ=(dimY-1)/2+j;
           bool masked=false;
-          for(int imask=0;imask<m_mask.size();++imask){
-            if(inBuffer[indexJ][indexI]==m_mask[imask]){
+          for(int imask=0;imask<m_noDataValues.size();++imask){
+            if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
               masked=true;
               break;
             }
@@ -327,41 +328,41 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
       switch(getFilterType(method)){
       case(filter2d::median):
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=stat.median(windowBuffer);
         break;
       case(filter2d::var):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=stat.var(windowBuffer);
         break;
       }
       case(filter2d::stdev):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=sqrt(stat.var(windowBuffer));
         break;
       }
       case(filter2d::mean):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=stat.mean(windowBuffer);
         break;
       }
       case(filter2d::min):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=stat.min(windowBuffer);
         break;
       }
       case(filter2d::ismin):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[centre])? 1:0;
         break;
@@ -370,7 +371,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
         double min=0;
         double max=0;
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else{
           stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
           if(min!=max)
@@ -382,21 +383,21 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
       }
       case(filter2d::max):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=stat.max(windowBuffer);
         break;
       }
       case(filter2d::ismax):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else
           outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[centre])? 1:0;
         break;
       }
       case(filter2d::order):{
         if(windowBuffer.empty())
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         else{
           double lbound=0;
           double ubound=dimX*dimY;
@@ -415,7 +416,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
         if(occurrence.size()==1)//all values in window must be the same
           outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
         else//favorize original value in case of ties
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         break;
       case(filter2d::heterog):{
         for(std::vector<double>::const_iterator wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
@@ -424,7 +425,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
           else if(*wit!=inBuffer[(dimY-1)/2][x])
             outBuffer[x/down]=1;
           else if(*wit==inBuffer[(dimY-1)/2][x]){//todo:wit mag niet central pixel zijn
-            outBuffer[x/down]=m_noValue;
+            outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
             break;
           }
         }
@@ -437,7 +438,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
             outBuffer[x/down]+=100.0*occurrence[*(vit++)]/windowBuffer.size();
         }
         else
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         break;
       }
       case(filter2d::majority):{
@@ -453,7 +454,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
             outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
         }
         else
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         break;
       }
       case(filter2d::threshold):{
@@ -466,7 +467,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
           }
         }
         else
-          outBuffer[x/down]=m_noValue;
+          outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
         break;
       }
       case(filter2d::mixed):{
@@ -516,7 +517,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
 //   }
 // };
 
-template<class T> void Filter2d::shift(const Vector2d<T>& input, Vector2d<T>& output, int offsetX, int offsetY, double randomSigma, RESAMPLE resample, bool verbose){
+template<class T> void Filter2d::shift(const Vector2d<T>& input, Vector2d<T>& output, double offsetX, double offsetY, double randomSigma, RESAMPLE resample, bool verbose){
   output.resize(input.nRows(),input.nCols());
   const gsl_rng_type *rangenType;
   gsl_rng *rangen;
@@ -646,8 +647,8 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
       double currentValue=inBuffer[(dimY-1)/2][x];
       std::vector<double> statBuffer;
       bool currentMasked=false;
-      for(int imask=0;imask<m_mask.size();++imask){
-        if(currentValue==m_mask[imask]){
+      for(int imask=0;imask<m_noDataValues.size();++imask){
+        if(currentValue==m_noDataValues[imask]){
           currentMasked=true;
           break;
         }
@@ -674,11 +675,11 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
                 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
               else
                 indexJ=(dimY-1)/2+j;
-            if(inBuffer[indexJ][indexI]==m_noValue)
+            if(inBuffer[indexJ][indexI]==(m_noDataValues.size())? m_noDataValues[0] : 0)
               continue;
             bool masked=false;
-            for(int imask=0;imask<m_mask.size();++imask){
-              if(inBuffer[indexJ][indexI]==m_mask[imask]){
+            for(int imask=0;imask<m_noDataValues.size();++imask){
+              if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
                 masked=true;
                 break;
               }
@@ -714,19 +715,19 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
             break;
           default:
             std::ostringstream ess;
-            ess << "Error:  morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << endl;
+            ess << "Error:  morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << std::endl;
             throw(ess.str());
             break;
           }
           if(output[y][x]&&m_class.size())
             output[y][x]=m_class[0];
           // else{
-          //   assert(m_mask.size());
-          //   output[x]=m_mask[0];
+          //   assert(m_noDataValues.size());
+          //   output[x]=m_noDataValues[0];
           // }
         }
         else
-          output[y][x]=m_noValue;
+          output[y][x]=(m_noDataValues.size())? m_noDataValues[0] : 0;
       }
     }
     progress=(1.0+y);
diff --git a/src/apps/pkascii2img.cc b/src/apps/pkascii2img.cc
index 643c0d5..5d6d0d3 100644
--- a/src/apps/pkascii2img.cc
+++ b/src/apps/pkascii2img.cc
@@ -23,6 +23,7 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include <assert.h>
 #include "imageclasses/ImgWriterGdal.h"
 
+using namespace std;
 
 int main(int argc, char *argv[])
 {
@@ -138,7 +139,14 @@ int main(int argc, char *argv[])
     assert(dy_opt.size());
     if(verbose_opt[0])
       cout << output_opt[0] << " is georeferenced." << endl;
-    imgWriter.setGeoTransform(ulx_opt[0],uly_opt[0],dx_opt[0],dy_opt[0],0,0);
+    double gt[6];
+    gt[0]=ulx_opt[0];
+    gt[1]=dx_opt[0];
+    gt[2]=0;
+    gt[3]=uly_opt[0];
+    gt[4]=0;
+    gt[5]=dy_opt[0];
+    imgWriter.setGeoTransform(gt);
     imgWriter.setProjectionProj4(projection_opt[0]);
   }
   else{
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index e1cd306..32e986d 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -28,6 +28,8 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "base/Optionpk.h"
 #include "algorithms/Egcs.h"
 
+using namespace std;
+
 int main(int argc, char *argv[])
 {
   Optionpk<string>  input_opt("i", "input", "Input image file(s). If input contains multiple images, a multi-band output is created");
@@ -356,7 +358,14 @@ int main(int argc, char *argv[])
       }
       if(description_opt.size())
 	imgWriter.setImageDescription(description_opt[0]);
-      imgWriter.setGeoTransform(cropulx,cropuly,dx,dy,0,0);
+      double gt[6];
+      gt[0]=cropulx;
+      gt[1]=dx;
+      gt[2]=0;
+      gt[3]=cropuly;
+      gt[4]=0;
+      gt[5]=dy;
+      imgWriter.setGeoTransform(gt);
       if(projection_opt.size()){
 	if(verbose_opt[0])
 	  cout << "projection: " << projection_opt[0] << endl;
@@ -405,8 +414,12 @@ int main(int argc, char *argv[])
       }
       double theMin=0;
       double theMax=0;
-      if(autoscale_opt.size())
+      if(autoscale_opt.size()){
 	imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax);
+	if(verbose_opt[0])
+	  cout << "minmax: " << theMin << ", " << theMax << endl;
+      }	
+      
       double readRow=0;
       double readCol=0;
       double lowerCol=0;
@@ -417,7 +430,11 @@ int main(int argc, char *argv[])
 	//convert irow to geo
 	imgWriter.image2geo(0,irow,x,y);
 	//lookup corresponding row for irow in this file
+	  //test
+	  // cout << "x: " << x << ", y: " << y << endl;
 	imgReader.geo2image(x,y,readCol,readRow);
+	//test
+	// cout << "readRow: " << readRow << ", readCol: " << readCol << flush << endl;
 	// double lowerCol=0;
 	// double upperCol=0;
 	vector<double> writeBuffer;
diff --git a/src/apps/pkdsm2shadow.cc b/src/apps/pkdsm2shadow.cc
index 58c25ec..dc050d6 100644
--- a/src/apps/pkdsm2shadow.cc
+++ b/src/apps/pkdsm2shadow.cc
@@ -30,6 +30,8 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "imageclasses/ImgReaderGdal.h"
 #include "imageclasses/ImgWriterGdal.h"
 
+using namespace std;
+
 /*------------------
   Main procedure
   ----------------*/
@@ -108,9 +110,9 @@ int main(int argc,char **argv) {
   }
   if(input.isGeoRef()){
     output.setProjection(input.getProjection());
-    double ulx,uly,deltaX,deltaY,rot1,rot2;
-    input.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
-    output.setGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
+    double gt[6];
+    input.getGeoTransform(gt);
+    output.setGeoTransform(gt);
   }
   if(input.getColorTable()!=NULL)
     output.setColorTable(input.getColorTable());
diff --git a/src/apps/pkdumpimg.cc b/src/apps/pkdumpimg.cc
index 1a74b11..2ce93f4 100644
--- a/src/apps/pkdumpimg.cc
+++ b/src/apps/pkdumpimg.cc
@@ -201,8 +201,9 @@ int main(int argc, char *argv[])
     cout << "selected lower right column in input image: " << lri << endl;
     cout << "selected lower right row in input image: " << lrj << endl;
   }
-
-  imgWriter.setGeoTransform(cropulx,cropuly,dx,dy,0,0);
+  double gt[6];
+  imgReader.getGeoTransform(gt);
+  imgWriter.setGeoTransform(gt);
   // imgWriter.setProjection(imgReader.getProjection());
 
   double readRow=0;
diff --git a/src/apps/pkenhance.cc b/src/apps/pkenhance.cc
index 09fd5da..1fb9539 100644
--- a/src/apps/pkenhance.cc
+++ b/src/apps/pkenhance.cc
@@ -24,6 +24,8 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "imageclasses/ImgReaderGdal.h"
 #include "imageclasses/ImgWriterGdal.h"
 
+using namespace std;
+
 int main(int argc,char **argv) {
   Optionpk<std::string> input_opt("i","input","input image file");
   Optionpk<string> reference_opt("r", "reference", "Reference image file");
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 2ae4261..97824ee 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -32,7 +32,7 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "imageclasses/ImgReaderGdal.h"
 #include "imageclasses/ImgWriterGdal.h"
 
-
+using namespace std;
 /*------------------
   Main procedure
   ----------------*/
@@ -40,17 +40,17 @@ 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 (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<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> 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);
+  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);
   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<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 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<short> nodata_opt("nodata", "nodata", "nodata value(s) for smoothnodata filter");
   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");
   Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
@@ -66,7 +66,7 @@ int main(int argc,char **argv) {
   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> l2_opt("l2","l2", "obtain shortest object length for linear feature",false,2);
   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);
@@ -76,7 +76,7 @@ int main(int argc,char **argv) {
     doProcess=input_opt.retrieveOption(argc,argv);
     output_opt.retrieveOption(argc,argv);
     disc_opt.retrieveOption(argc,argv);
-    angle_opt.retrieveOption(argc,argv);
+    // angle_opt.retrieveOption(argc,argv);
     method_opt.retrieveOption(argc,argv);
     resample_opt.retrieveOption(argc,argv);
     dimX_opt.retrieveOption(argc,argv);
@@ -87,7 +87,7 @@ int main(int argc,char **argv) {
     family_opt.retrieveOption(argc,argv);
     class_opt.retrieveOption(argc,argv);
     threshold_opt.retrieveOption(argc,argv);
-    mask_opt.retrieveOption(argc,argv);
+    nodata_opt.retrieveOption(argc,argv);
     tap_opt.retrieveOption(argc,argv);
     tapz_opt.retrieveOption(argc,argv);
     fwhm_opt.retrieveOption(argc,argv);
@@ -116,11 +116,24 @@ int main(int argc,char **argv) {
     exit(0);//help was invoked, stop processing
   }
 
+  //not implemented yet, must debug first...
+  vector<double> angle_opt;
+
   ImgReaderGdal input;
   ImgWriterGdal output;
-  assert(input_opt.size());
+  if(input_opt.empty()){
+    cerr << "Error: no input file selected, use option -i" << endl;
+    exit(1);
+  }
+  if(output_opt.empty()){
+    cerr << "Error: no output file selected, use option -o" << endl;
+    exit(1);
+  }
+  if(method_opt.empty()){
+    cerr << "Error: no filter selected, use option -f" << endl;
+    exit(1);
+  }
   input.open(input_opt[0]);
-  // output.open(output_opt[0],input);
   GDALDataType theType=GDT_Unknown;
   if(verbose_opt[0])
     cout << "possible output data types: ";
@@ -148,7 +161,6 @@ int main(int argc,char **argv) {
     option_opt.push_back(theInterleave);
   }
   try{
-    assert(output_opt.size());
     if(filter2d::Filter2d::getFilterType(method_opt[0])==filter2d::mrf){
       assert(class_opt.size()>1);
       if(verbose_opt[0])
@@ -183,9 +195,9 @@ int main(int argc,char **argv) {
   }
   if(input.isGeoRef()){
     output.setProjection(input.getProjection());
-    double ulx,uly,deltaX,deltaY,rot1,rot2;
-    input.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
-    output.setGeoTransform(ulx,uly,deltaX*down_opt[0],deltaY*down_opt[0],rot1,rot2);
+    double gt[6];
+    input.getGeoTransform(gt);
+    output.setGeoTransform(gt);
   }
   if(colorTable_opt.size()){
     if(colorTable_opt[0]!="none"){
@@ -214,13 +226,13 @@ int main(int argc,char **argv) {
     if(verbose_opt[0])
       std::cout<< std::endl;
   }
-  if(mask_opt.size()){
+  if(nodata_opt.size()){
     if(verbose_opt[0])
       std::cout<< "mask values: ";
-    for(int imask=0;imask<mask_opt.size();++imask){
+    for(int imask=0;imask<nodata_opt.size();++imask){
       if(verbose_opt[0])
-        std::cout<< mask_opt[imask] << " ";
-      filter2d.pushMask(mask_opt[imask]);
+        std::cout<< nodata_opt[imask] << " ";
+      filter2d.pushNoDataValue(nodata_opt[imask]);
     }
     if(verbose_opt[0])
       std::cout<< std::endl;
@@ -350,6 +362,10 @@ int main(int argc,char **argv) {
   else{
     switch(filter2d::Filter2d::getFilterType(method_opt[0])){
     case(filter2d::dilate):
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for morphological operator" << std::endl;
+	exit(1);
+      }
       if(dimZ_opt.size()){
         if(verbose_opt[0])
           std::cout<< "1-D filtering: dilate" << std::endl;
@@ -360,6 +376,10 @@ int main(int argc,char **argv) {
       }
       break;
     case(filter2d::erode):
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for morphological operator" << std::endl;
+	exit(1);
+      }
       if(dimZ_opt.size()>0){
         if(verbose_opt[0])
           std::cout<< "1-D filtering: dilate" << std::endl;
@@ -370,6 +390,10 @@ int main(int argc,char **argv) {
       }
       break;
     case(filter2d::close):{//closing
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for morphological operator" << std::endl;
+	exit(1);
+      }
       ostringstream tmps;
       tmps << "/tmp/dilation_" << getpid() << ".tif";
       ImgWriterGdal tmpout;
@@ -402,6 +426,10 @@ int main(int argc,char **argv) {
       break;
     }
     case(filter2d::open):{//opening
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for morphological operator" << std::endl;
+	exit(1);
+      }
       ostringstream tmps;
       tmps << "/tmp/erosion_" << getpid() << ".tif";
       ImgWriterGdal tmpout;
@@ -436,7 +464,11 @@ int main(int argc,char **argv) {
       filter2d.doit(input,output,"heterog",dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
       break;
     }
-    case(filter2d::shift):{//spatially heterogeneous
+    case(filter2d::shift):{//shift
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for shift operator" << std::endl;
+	exit(1);
+      }
       assert(input.nrOfBand());
       assert(input.nrOfCol());
       assert(input.nrOfRow());
@@ -449,6 +481,10 @@ int main(int argc,char **argv) {
       break;
     }
     case(filter2d::linearfeature):{
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for linear feature" << std::endl;
+	exit(1);
+      }
       assert(input.nrOfBand());
       assert(input.nrOfCol());
       assert(input.nrOfRow());
@@ -458,6 +494,7 @@ int main(int argc,char **argv) {
       if(verbose_opt[0])
 	std::cout << "using angle " << theAngle << std::endl;
       try{
+	//using an angle step of 5 degrees and no maximum distance
         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){
@@ -496,6 +533,10 @@ int main(int argc,char **argv) {
       break;
     }
     case(filter2d::sobelx):{//Sobel edge detection in X
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
+	exit(1);
+      }
       Vector2d<double> theTaps(3,3);
       theTaps[0][0]=-1.0;
       theTaps[0][1]=0.0;
@@ -511,6 +552,10 @@ int main(int argc,char **argv) {
       break;
     }
     case(filter2d::sobely):{//Sobel edge detection in Y
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
+	exit(1);
+      }
       Vector2d<double> theTaps(3,3);
       theTaps[0][0]=1.0;
       theTaps[0][1]=2.0;
@@ -526,6 +571,10 @@ int main(int argc,char **argv) {
       break;
     }
     case(filter2d::sobelxy):{//Sobel edge detection in XY
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
+	exit(1);
+      }
       Vector2d<double> theTaps(3,3);
       theTaps[0][0]=0.0;
       theTaps[0][1]=1.0;
@@ -541,6 +590,10 @@ int main(int argc,char **argv) {
       break;
     }
     case(filter2d::sobelyx):{//Sobel edge detection in XY
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for sobel edge detection" << std::endl;
+	exit(1);
+      }
       Vector2d<double> theTaps(3,3);
       theTaps[0][0]=2.0;
       theTaps[0][1]=1.0;
@@ -556,29 +609,61 @@ int main(int argc,char **argv) {
       break;
     }
     case(filter2d::smooth):{//Smoothing filter
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for this filter" << std::endl;
+	exit(1);
+      }
       filter2d.smooth(input,output,dimX_opt[0],dimY_opt[0]);
       break;
     }
     case(filter2d::smoothnodata):{//Smoothing filter
+      if(down_opt[0]!=1){
+	std::cerr << "Error: down option not supported for this filter" << std::endl;
+	exit(1);
+      }
       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]);
       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]);
       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]);
       break;
     case(filter2d::threshold):
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
index 6ef0a82..6a6980c 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -33,6 +33,8 @@ extern "C" {
 #include <gsl/gsl_sort.h>
 }
 
+using namespace std;
+
 /*------------------
   Main procedure
   ----------------*/
diff --git a/src/apps/pkgetmask.cc b/src/apps/pkgetmask.cc
index 91f4e4e..041dd58 100644
--- a/src/apps/pkgetmask.cc
+++ b/src/apps/pkgetmask.cc
@@ -139,9 +139,9 @@ int main(int argc,char **argv) {
     imgWriter.setColorTable(imgReader.getColorTable());
   if(imgReader.isGeoRef()){
     imgWriter.setProjection(imgReader.getProjection());
-    double ulx,uly,lrx,lry;
-    imgReader.getBoundingBox(ulx,uly,lrx,lry);
-    imgWriter.setGeoTransform(ulx,uly,imgReader.getDeltaX(),imgReader.getDeltaY(),0,0);
+    double gt[6];
+    imgReader.getGeoTransform(gt);
+    imgWriter.setGeoTransform(gt);//ulx,uly,imgReader.getDeltaX(),imgReader.getDeltaY(),0,0);
   }
   vector<char> writeBuffer(imgWriter.nrOfCol());
   for(int irow=0;irow<imgReader.nrOfRow();++irow){
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index 00823fc..2aa9b6f 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -329,9 +329,7 @@ int main(int argc, char *argv[])
         std::cout << " -a_srs none" << " ";
     }
     if(geo_opt[0]&&!read_opt[0]){
-      double ulx,uly,deltaX,deltaY,rot1,rot2;
-      imgReader.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
-      std::cout << " --geo " << std::setprecision(12) << ulx << " " << uly << " " << deltaX << " " << deltaY << " " << rot1 << " " << rot2 << " ";
+      std::cout << " --geo " << std::setprecision(12) << imgReader.getGeoTransform();
     }
     if(interleave_opt[0]){
       std::cout << " --interleave " << imgReader.getInterleave() << " ";
diff --git a/src/apps/pklas2img.cc b/src/apps/pklas2img.cc
index 6c2b7d2..f7af41a 100644
--- a/src/apps/pklas2img.cc
+++ b/src/apps/pklas2img.cc
@@ -26,6 +26,8 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "algorithms/StatFactory.h"
 #include "algorithms/Filter2d.h"
 
+using namespace std;
+
 int main(int argc,char **argv) {
   Optionpk<string> input_opt("i", "input", "Input las file");
   Optionpk<short> nodata_opt("nodata", "nodata", "nodata value to put in image if not valid", 0);
@@ -203,7 +205,14 @@ int main(int argc,char **argv) {
   outputWriter.open(output_opt[0],ncol,nrow,nband,theType,oformat_opt[0],option_opt);
   outputWriter.GDALSetNoDataValue(nodata_opt[0]);
   //set projection
-  outputWriter.setGeoTransform(minULX,maxULY,dx_opt[0],dy_opt[0],0,0);
+  double gt[6];
+  gt[0]=minULX;
+  gt[1]=dx_opt[0];
+  gt[2]=0;
+  gt[3]=maxULY;
+  gt[4]=0;
+  gt[5]=dy_opt[0];
+  outputWriter.setGeoTransform(gt);
   if(projection_opt.size()){
     string projectionString=outputWriter.setProjectionProj4(projection_opt[0]);
     if(verbose_opt[0])
@@ -411,7 +420,7 @@ int main(int argc,char **argv) {
     int dimx=dimx_opt[0];
     int dimy=dimy_opt[0];
     filter2d::Filter2d morphFilter;
-    morphFilter.setNoValue(0);
+    // morphFilter.setNoValue(0);
     Vector2d<float> currentOutput=outputData;
     int iteration=1;
     while(nchange&&iteration<maxIter_opt[0]){
@@ -436,7 +445,7 @@ int main(int argc,char **argv) {
     int dimx=dimx_opt[0];
     int dimy=dimy_opt[0];
     filter2d::Filter2d theFilter;
-    theFilter.setNoValue(0);
+    // theFilter.setNoValue(0);
     Vector2d<float> currentOutput=outputData;
     double hThreshold=hThreshold_opt[0];
     int iteration=1;
@@ -470,7 +479,7 @@ int main(int argc,char **argv) {
     if(composite_opt[0]!="min")
       std::cout << "Warning: composite option is not set to min!" << std::endl;
     filter2d::Filter2d morphFilter;
-    morphFilter.setNoValue(0);
+    // morphFilter.setNoValue(0);
     Vector2d<float> filterInput=outputData;
     try{
       morphFilter.morphology(outputData,filterInput,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
@@ -485,7 +494,7 @@ int main(int argc,char **argv) {
     if(composite_opt[0]!="max")
       std::cout << "Warning: composite option is not set to max!" << std::endl;
     filter2d::Filter2d morphFilter;
-    morphFilter.setNoValue(0);
+    // morphFilter.setNoValue(0);
     Vector2d<float> filterInput=outputData;
     try{
       morphFilter.morphology(outputData,filterInput,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
diff --git a/src/apps/pkmosaic.cc b/src/apps/pkmosaic.cc
index a5a904b..4be22b3 100644
--- a/src/apps/pkmosaic.cc
+++ b/src/apps/pkmosaic.cc
@@ -395,7 +395,15 @@ int main(int argc, char *argv[])
   }
   if(description_opt.size())
     imgWriter.setImageDescription(description_opt[0]);
-  imgWriter.setGeoTransform(minULX,maxULY,dx,dy,0,0);
+  double gt[6];
+  gt[0]=minULX;
+  gt[1]=dx;
+  gt[2]=0;
+  gt[3]=maxULY;
+  gt[4]=0;
+  gt[5]=dy;
+  imgWriter.setGeoTransform(gt);
+  // imgWriter.setGeoTransform(minULX,maxULY,dx,dy,0,0);
   if(projection_opt.size()){
     if(verbose_opt[0])
       cout << "projection: " << projection_opt[0] << endl;
diff --git a/src/imageclasses/ImgReaderGdal.cc b/src/imageclasses/ImgReaderGdal.cc
index 963012d..f7d970a 100644
--- a/src/imageclasses/ImgReaderGdal.cc
+++ b/src/imageclasses/ImgReaderGdal.cc
@@ -59,22 +59,24 @@ void ImgReaderGdal::setCodec()//double magicX, double magicY)
   m_isGeoRef=( static_cast<std::string>(m_gds->GetProjectionRef())  != "" );
   // m_magic_x=magicX;
   // m_magic_y=magicY;
-  if(m_isGeoRef){
-    double adfGeoTransform[6];
-    if( m_gds->GetGeoTransform( adfGeoTransform ) == CE_None )
-    {
-      m_ulx=adfGeoTransform[0];
-      m_uly=adfGeoTransform[3];
-      m_delta_x=adfGeoTransform[1];
-      m_delta_y=-adfGeoTransform[5];
-    }
-  }
-  else{
-    m_ulx=0;
-    m_uly=nrOfRow();
-    m_delta_x=1;
-    m_delta_y=1;
-  }
+  double adfGeoTransform[6];
+  m_gds->GetGeoTransform( adfGeoTransform );
+  // if( m_gds->GetGeoTransform( adfGeoTransform ) == CE_None ){
+  m_gt[0]=adfGeoTransform[0];
+  m_gt[1]=adfGeoTransform[1];
+  m_gt[2]=adfGeoTransform[2];
+  m_gt[3]=adfGeoTransform[3];
+  m_gt[4]=adfGeoTransform[4];
+  m_gt[5]=adfGeoTransform[5];
+  // }
+  // else{
+  //   m_gt[0]=0;
+  //   m_gt[1]=1;
+  //   m_gt[2]=0;
+  //   m_gt[3]=0;
+  //   m_gt[4]=0;
+  //   m_gt[5]=1;
+  // }
 }
 
 std::string ImgReaderGdal::getProjection(void) const 
@@ -120,35 +122,44 @@ std::string ImgReaderGdal::getDriverDescription() const
   return m_gds->GetDriver()->GetDescription();
 }
 
-void ImgReaderGdal::getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const
-{
-  double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
-  m_gds->GetGeoTransform(adfGeoTransform);
-  ulx=adfGeoTransform[0];
-  deltaX=adfGeoTransform[1];
-  rot1=adfGeoTransform[2];
-  uly=adfGeoTransform[3];
-  rot2=adfGeoTransform[4];
-  deltaY=-adfGeoTransform[5];//convention of GDAL!
+void ImgReaderGdal::getGeoTransform(double* gt) const{
+  m_gds->GetGeoTransform(gt);
 }
 
+// void ImgReaderGdal::getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const
+// {
+//   double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
+//   m_gds->GetGeoTransform(adfGeoTransform);
+//   ulx=adfGeoTransform[0];
+//   deltaX=adfGeoTransform[1];
+//   rot1=adfGeoTransform[2];
+//   uly=adfGeoTransform[3];
+//   rot2=adfGeoTransform[4];
+//   deltaY=-adfGeoTransform[5];//convention of GDAL!
+// }
+
 std::string ImgReaderGdal::getGeoTransform() const
 {
-  if(!isGeoRef())
-    return("");
-  else{
-    double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
-    m_gds->GetGeoTransform(adfGeoTransform);
-    double ulx=adfGeoTransform[0];
-    double deltaX=adfGeoTransform[1];
-    double rot1=adfGeoTransform[2];
-    double uly=adfGeoTransform[3];
-    double rot2=adfGeoTransform[4];
-    double deltaY=-adfGeoTransform[5];//convention of GDAL!
-    std::ostringstream s;
-    s << "[" << ulx << "," << deltaX << "," << rot1 << "," << uly << "," << rot2 << "," << -deltaY << "]";
-    return(s.str());
-  }
+  double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+  m_gds->GetGeoTransform(gt);
+  std::ostringstream s;
+  s << "[" << gt[0] << "," << gt[1] << "," << gt[2] << "," << gt[3] << "," << gt[4] << "," << -gt[5] << "]";
+  return(s.str());
+  // if(!isGeoRef())
+  //   return("");
+  // else{
+  //   double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
+  //   m_gds->GetGeoTransform(adfGeoTransform);
+  //   double ulx=adfGeoTransform[0];
+  //   double deltaX=adfGeoTransform[1];
+  //   double rot1=adfGeoTransform[2];
+  //   double uly=adfGeoTransform[3];
+  //   double rot2=adfGeoTransform[4];
+  //   double deltaY=-adfGeoTransform[5];//convention of GDAL!
+  //   std::ostringstream s;
+  //   s << "[" << ulx << "," << deltaX << "," << rot1 << "," << uly << "," << rot2 << "," << -deltaY << "]";
+  //   return(s.str());
+  // }
 }
 
 char** ImgReaderGdal::getMetadata()
@@ -286,9 +297,10 @@ bool ImgReaderGdal::geo2image(double x, double y, double& i, double& j) const
   //adfGeotransform[3]: ULY (upper left Y coordinate)
   //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
   //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+
   double denom=(gt[1]-gt[2]*gt[4]/gt[5]);
   double eps=0.00001;
-  if(denom>eps){
+  if(fabs(denom)>eps){
     i=(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom;
     j=(y-gt[3]-gt[4]*(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom)/gt[5];
   }
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index b46a659..7e51074 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -27,8 +27,6 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "gdal_priv.h"
 #include "base/Vector2d.h"
 
-// using namespace std;
-
 enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };
 
 //--------------------------------------------------------------------------
@@ -47,16 +45,17 @@ public:
   std::string getProjection(void) const;
   std::string getProjectionRef(void) const;
   std::string getGeoTransform() const;
-  void getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const;
+  void getGeoTransform(double* gt) const;
+  /* void getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const; */
   std::string getDescription() const;
   std::string getMetadataItem() const;
   std::string getImageDescription() const;
   bool getBoundingBox (double& ulx, double& uly, double& lrx, double& lry) const;
   bool getCentrePos(double& x, double& y) const;
-  double getUlx() const {return (m_isGeoRef)? m_ulx : 0;};
-  double getUly() const {return (m_isGeoRef)? m_uly : nrOfRow()-1;};
-  double getLrx() const {return (m_isGeoRef)? m_ulx+nrOfCol()*m_delta_x : nrOfCol()-1;};
-  double getLry() const {return (m_isGeoRef)? m_uly-nrOfRow()*m_delta_y : 0;};
+  double getUlx() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(ulx);};
+  double getUly() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(uly);};
+  double getLrx() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(lrx);};
+  double getLry() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(lry);};
   // bool getMagicPixel(double& magicX, double& magicY) const {magicX=m_magic_x;magicY=m_magic_y;};
   int getNoDataValues(std::vector<double>& noDataValues) const;
   bool isNoData(double value) const{return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();};
@@ -66,8 +65,8 @@ public:
   bool covers(double ulx, double  uly, double lrx, double lry) const;
   bool geo2image(double x, double y, double& i, double& j) const;
   bool image2geo(double i, double j, double& x, double& y) const;
-  double getDeltaX(void) const {return m_delta_x;};
-  double getDeltaY(void) const {return m_delta_y;};
+  double getDeltaX(void) const {double gt[6];getGeoTransform(gt);return gt[1];};
+  double getDeltaY(void) const {double gt[6];getGeoTransform(gt);return -gt[5];};
   template<typename T> void readData(T& value, const GDALDataType& dataType, int col, int row, int band=0) const;
   template<typename T> void readData(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int row, int band=0) const;
   template<typename T> void readData(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, double row, int band=0, RESAMPLE resample=0) const;
@@ -103,12 +102,11 @@ protected:
   int m_ncol;
   int m_nrow;
   int m_nband;
-  double m_ulx;
-  double m_uly;
-  double m_delta_x;
-  double m_delta_y;
-  // double m_magic_x;
-  // double m_magic_y;
+  double m_gt[6];
+  /* double m_ulx; */
+  /* double m_uly; */
+  /* double m_delta_x; */
+  /* double m_delta_y; */
   bool m_isGeoRef;
   std::vector<double> m_noDataValues;
 };
diff --git a/src/imageclasses/ImgWriterGdal.cc b/src/imageclasses/ImgWriterGdal.cc
index 0b8da57..5569bb3 100644
--- a/src/imageclasses/ImgWriterGdal.cc
+++ b/src/imageclasses/ImgWriterGdal.cc
@@ -44,7 +44,7 @@ ImgWriterGdal::~ImgWriterGdal(void)
 }
 
 //---------------------------------------------------------------------------
-void ImgWriterGdal::open(const string& filename, const ImgReaderGdal& imgSrc, const vector<string>& options)
+void ImgWriterGdal::open(const std::string& filename, const ImgReaderGdal& imgSrc, const std::vector<std::string>& options)
 {
   m_isGeoRef=imgSrc.isGeoRef();
   m_filename=filename;
@@ -59,7 +59,7 @@ void ImgWriterGdal::open(const string& filename, const ImgReaderGdal& imgSrc, co
   setCodec(imgSrc);
 }
 
-// void ImgWriterGdal::open(const string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const string& imageType, const string& interleave, const string& compression, int magicX, int magicY)
+// void ImgWriterGdal::open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType, const std::string& interleave, const std::string& compression, int magicX, int magicY)
 // {
 //   m_isGeoRef=false;
 //   m_filename = filename;
@@ -74,7 +74,7 @@ void ImgWriterGdal::open(const string& filename, const ImgReaderGdal& imgSrc, co
 //   setCodec(imageType);
 // }
 
-void ImgWriterGdal::open(const string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const string& imageType, const vector<string>& options)
+void ImgWriterGdal::open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType, const std::vector<std::string>& options)
 {
   m_isGeoRef=false;
   m_filename = filename;
@@ -95,7 +95,7 @@ void ImgWriterGdal::close(void)
 {
   GDALClose(m_gds);
   char **papszOptions=NULL;
-  for(vector<string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
+  for(std::vector<std::string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
     papszOptions=CSLAddString(papszOptions,optionIt->c_str());
   CSLDestroy(papszOptions);
 }
@@ -106,7 +106,7 @@ void ImgWriterGdal::setCodec(const ImgReaderGdal& imgSrc){
   GDALDriver *poDriver;
   poDriver = GetGDALDriverManager()->GetDriverByName(imgSrc.getDriverDescription().c_str());
   if( poDriver == NULL ){
-    string errorString="FileOpenError";
+    std::string errorString="FileOpenError";
     throw(errorString);
   }
   char **papszMetadata;
@@ -114,21 +114,21 @@ void ImgWriterGdal::setCodec(const ImgReaderGdal& imgSrc){
   //todo: try and catch if CREATE is not supported (as in PNG)
   assert( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ));
   char **papszOptions=NULL;
-  for(vector<string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
+  for(std::vector<std::string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
     papszOptions=CSLAddString(papszOptions,optionIt->c_str());
   // char **papszOptions=NULL;
-  // ostringstream compressList;
+  // std::ostringstream compressList;
   // compressList << "COMPRESS=" << m_compression;
   // papszOptions = CSLAddString(papszOptions,(compressList.str()).c_str());
-  // ostringstream interleaveList;
+  // std::ostringstream interleaveList;
   // interleaveList << "INTERLEAVE=" << m_interleave;
   // papszOptions = CSLAddString(papszOptions,(interleaveList.str()).c_str());
   m_gds=poDriver->Create(m_filename.c_str(),m_ncol,m_nrow,m_nband,m_type,papszOptions);
   if(imgSrc.isGeoRef()){
     setProjection(imgSrc.getProjection());
-    double ulx,uly,deltaX,deltaY,rot1,rot2;
-    imgSrc.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
-    setGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
+    double gt[6];
+    imgSrc.getGeoTransform(gt);
+    setGeoTransform(gt);
   }
   m_gds->SetMetadata(imgSrc.getMetadata() ); 
   m_gds->SetMetadataItem( "TIFFTAG_DOCUMENTNAME", m_filename.c_str());
@@ -141,8 +141,8 @@ void ImgWriterGdal::setCodec(const ImgReaderGdal& imgSrc){
 
   time_t tim=time(NULL);
   tm *now=localtime(&tim);
-  ostringstream datestream;
-  //date string must be 20 characters long...
+  std::ostringstream datestream;
+  //date std::string must be 20 characters long...
   datestream << now->tm_year+1900;
   if(now->tm_mon+1<10)
     datestream << ":0" << now->tm_mon+1;
@@ -165,9 +165,9 @@ void ImgWriterGdal::setCodec(const ImgReaderGdal& imgSrc){
   else
     datestream << ":" << now->tm_sec;
   m_gds->SetMetadataItem( "TIFFTAG_DATETIME", datestream.str().c_str());
-//   list<string> lmeta;
+//   list<std::string> lmeta;
 //   imgReader.getMetadata(lmeta);
-//   list<string>::const_iterator lit=lmeta.begin();
+//   list<std::string>::const_iterator lit=lmeta.begin();
 //   while(lit!=lmeta.end()){
 //     cout << *lit << endl;
 //     ++lit;
@@ -178,13 +178,13 @@ void ImgWriterGdal::setCodec(const ImgReaderGdal& imgSrc){
     setColorTable(imgSrc.getColorTable());
 }
 
-void ImgWriterGdal::setCodec(const string& imageType)
+void ImgWriterGdal::setCodec(const std::string& imageType)
 {
   GDALAllRegister();
   GDALDriver *poDriver;
   poDriver = GetGDALDriverManager()->GetDriverByName(imageType.c_str());
   if( poDriver == NULL ){
-    ostringstream s;
+    std::ostringstream s;
     s << "FileOpenError (" << imageType << ")";
     throw(s.str());
   }
@@ -193,12 +193,12 @@ void ImgWriterGdal::setCodec(const string& imageType)
   //todo: try and catch if CREATE is not supported (as in PNG)
   assert( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ));
   char **papszOptions=NULL;
-  for(vector<string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
+  for(std::vector<std::string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
     papszOptions=CSLAddString(papszOptions,optionIt->c_str());
-  // ostringstream compressList;
+  // std::ostringstream compressList;
   // compressList << "COMPRESS=" << m_compression;
   // papszOptions = CSLAddString(papszOptions,(compressList.str()).c_str());
-  // ostringstream interleaveList;
+  // std::ostringstream interleaveList;
   // interleaveList << "INTERLEAVE=" << m_interleave;
   // papszOptions = CSLAddString(papszOptions,(interleaveList.str()).c_str());
   m_gds=poDriver->Create(m_filename.c_str(),m_ncol,m_nrow,m_nband,m_type,papszOptions);
@@ -215,8 +215,8 @@ void ImgWriterGdal::setCodec(const string& imageType)
 
   time_t tim=time(NULL);
   tm *now=localtime(&tim);
-  ostringstream datestream;
-  //date string must be 20 characters long...
+  std::ostringstream datestream;
+  //date std::string must be 20 characters long...
   datestream << now->tm_year+1900;
   if(now->tm_mon+1<10)
     datestream << ":0" << now->tm_mon+1;
@@ -248,10 +248,10 @@ void ImgWriterGdal::setMetadata(char** metadata)
   m_gds->SetMetadata(metadata); 
 }
 
-string ImgWriterGdal::getProjection(void) const 
+std::string ImgWriterGdal::getProjection(void) const 
 {
   assert(m_gds);
-  string theProjection=m_gds->GetProjectionRef();
+  std::string theProjection=m_gds->GetProjectionRef();
   //due to error in Gdal? AUTHORITY fields do not seem to work!
   // size_t startpos,endpos;
   // while((startpos=theProjection.find(",AUTHORITY"))!=string::npos){
@@ -262,33 +262,47 @@ string ImgWriterGdal::getProjection(void) const
 }
 
 //---------------------------------------------------------------------------
-void ImgWriterGdal::setGeoTransform(double ulx, double uly, double deltaX, double deltaY, double rot1, double rot2)
-{
+void ImgWriterGdal::setGeoTransform(double* gt){
   m_isGeoRef=true;
-  m_ulx=ulx;
-  m_uly=uly;
-  m_delta_x=deltaX;
-  m_delta_y=deltaY;
-  double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
-  adfGeoTransform[0]=ulx;
-  adfGeoTransform[1]=deltaX;
-  adfGeoTransform[2]=rot1;
-  adfGeoTransform[3]=uly;
-  adfGeoTransform[4]=rot2;
-  adfGeoTransform[5]=-deltaY;//convention of GDAL!
+  m_gt[0]=gt[0];
+  m_gt[1]=gt[1];
+  m_gt[2]=gt[2];
+  m_gt[3]=gt[3];
+  m_gt[4]=gt[4];
+  m_gt[5]=gt[5];
   if(m_gds)
-    m_gds->SetGeoTransform(adfGeoTransform);
+    m_gds->SetGeoTransform(m_gt);
 }
 
+// void ImgWriterGdal::setGeoTransform(double ulx, double uly, double deltaX, double deltaY, double rot1, double rot2)
+// {
+//   m_isGeoRef=true;
+//   m_ulx=ulx;
+//   m_uly=uly;
+//   m_delta_x=deltaX;
+//   m_delta_y=deltaY;
+//   double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
+//   adfGeoTransform[0]=ulx;
+//   adfGeoTransform[1]=deltaX;
+//   adfGeoTransform[2]=rot1;
+//   adfGeoTransform[3]=uly;
+//   adfGeoTransform[4]=rot2;
+//   adfGeoTransform[5]=-deltaY;//convention of GDAL!
+//   if(m_gds)
+//     m_gds->SetGeoTransform(adfGeoTransform);
+// }
+
 void ImgWriterGdal::copyGeoTransform(const ImgReaderGdal& imgSrc)
 {
   setProjection(imgSrc.getProjection());
-  double ulx,uly,deltaX,deltaY,rot1,rot2;
-  imgSrc.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
-  setGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
+  double gt[6];
+  imgSrc.getGeoTransform(gt);
+  setGeoTransform(gt);
+  // imgSrc.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
+  // setGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
 }
 
-string ImgWriterGdal::setProjectionProj4(const string& projection)
+std::string ImgWriterGdal::setProjectionProj4(const std::string& projection)
 {
   if(!m_isGeoRef)
     m_isGeoRef=true;
@@ -311,11 +325,11 @@ string ImgWriterGdal::setProjectionProj4(const string& projection)
     //     OSRExportToWkt( hSRS, &pszResult );  
     // else  
     // {  
-    //     ostringstream s;
+    //     std::ostringstream s;
     //     s << "Error in set projection " << projection;
     //     throw(s.str());
     // }  
-    // string theProjection=pszResult;
+    // std::string theProjection=pszResult;
     // assert(m_gds);
     // m_gds->SetProjection(theProjection.c_str());
     // OSRDestroySpatialReference( hSRS );  
@@ -323,7 +337,7 @@ string ImgWriterGdal::setProjectionProj4(const string& projection)
     // return theProjection;  
 }
 
-void ImgWriterGdal::setProjection(const string& projection)
+void ImgWriterGdal::setProjection(const std::string& projection)
 {
   if(!m_isGeoRef)
     m_isGeoRef=true;
@@ -335,9 +349,9 @@ void ImgWriterGdal::setProjection(const string& projection)
 }
 
 //default projection: ETSR-LAEA
-string ImgWriterGdal::setProjection(void)
+std::string ImgWriterGdal::setProjection(void)
 {
-  string theProjection;
+  std::string theProjection;
   OGRSpatialReference oSRS;
   char *pszSRS_WKT = NULL;
   //// oSRS.importFromEPSG(3035);
@@ -424,7 +438,7 @@ bool ImgWriterGdal::geo2image(double x, double y, double& i, double& j) const
   //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
   double denom=(gt[1]-gt[2]*gt[4]/gt[5]);
   double eps=0.00001;
-  if(denom>eps){
+  if(fabs(denom)>eps){
     i=(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom;
     j=(y-gt[3]-gt[4]*(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom)/gt[5];
   }
@@ -485,59 +499,59 @@ bool ImgWriterGdal::covers(double ulx, double  uly, double lrx, double lry) cons
   return((ulx < theLRX)&&(lrx > theULX)&&(lry < theULY)&&(uly > theLRY));
 }
 
-string ImgWriterGdal::getGeoTransform() const
+std::string ImgWriterGdal::getGeoTransform() const
 {
-  double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
-  double ulx;
-  double deltaX;
-  double rot1;
-  double uly;
-  double rot2;
-  double deltaY;
-  if(m_gds){
-    m_gds->GetGeoTransform(adfGeoTransform);
-    ulx=adfGeoTransform[0];
-    deltaX=adfGeoTransform[1];
-    rot1=adfGeoTransform[2];
-    uly=adfGeoTransform[3];
-    rot2=adfGeoTransform[4];
-    deltaY=-adfGeoTransform[5];//convention of GDAL!
-  }
+  double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+  if(m_gds)
+    m_gds->GetGeoTransform(gt);
   else{//virtual writer
-    ulx=m_ulx;
-    uly=m_uly;
-    deltaX=m_delta_x;
-    deltaY=m_delta_y;
-    rot1=0;
-    rot2=0;
+    gt[0]=m_gt[0];
+    gt[1]=m_gt[1];
+    gt[2]=m_gt[2];
+    gt[3]=m_gt[3];
+    gt[4]=m_gt[4];
+    gt[5]=m_gt[5];
   }
-  ostringstream s;
-  s << "[" << ulx << "," << deltaX << "," << rot1 << "," << uly << "," << rot2 << "," << -deltaY << "]";
+  std::ostringstream s;
+  s << "[" << gt[0] << "," << gt[1] << "," << gt[2] << "," << gt[3] << "," << gt[4] << "," << gt[5] << "]";
   return(s.str());
 }
 
-void ImgWriterGdal::getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const
-{
-  if(m_gds){
-    double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
-    m_gds->GetGeoTransform(adfGeoTransform);
-    ulx=adfGeoTransform[0];
-    deltaX=adfGeoTransform[1];
-    rot1=adfGeoTransform[2];
-    uly=adfGeoTransform[3];
-    rot2=adfGeoTransform[4];
-    deltaY=-adfGeoTransform[5];//convention of GDAL!
-  }
+void ImgWriterGdal::getGeoTransform(double* gt) const{
+  if(m_gds)
+    m_gds->GetGeoTransform(gt);
   else{//virtual writer
-    ulx=m_ulx;
-    uly=m_uly;
-    deltaX=m_delta_x;
-    deltaY=m_delta_y;
-    rot1=0;
-    rot2=0;
+    gt[0]=m_gt[0];
+    gt[1]=m_gt[1];
+    gt[2]=m_gt[2];
+    gt[3]=m_gt[3];
+    gt[4]=m_gt[4];
+    gt[5]=m_gt[5];
   }
 }
 
+// void ImgWriterGdal::getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const
+// {
+//   if(m_gds){
+//     double adfGeoTransform[6];// { 444720, 30, 0, 3751320, 0, -30 };
+//     m_gds->GetGeoTransform(adfGeoTransform);
+//     ulx=adfGeoTransform[0];
+//     deltaX=adfGeoTransform[1];
+//     rot1=adfGeoTransform[2];
+//     uly=adfGeoTransform[3];
+//     rot2=adfGeoTransform[4];
+//     deltaY=-adfGeoTransform[5];//convention of GDAL!
+//   }
+//   else{//virtual writer
+//     ulx=m_ulx;
+//     uly=m_uly;
+//     deltaX=m_delta_x;
+//     deltaY=m_delta_y;
+//     rot1=0;
+//     rot2=0;
+//   }
+// }
+
 GDALDataType ImgWriterGdal::getDataType(int band) const
 {
   assert(band<m_nband+1);
@@ -551,17 +565,17 @@ GDALRasterBand* ImgWriterGdal::getRasterBand(int band)
 }
 
 //filename is ascii file containing 5 columns: index R G B ALFA (0:transparent, 255:solid)
-void ImgWriterGdal::setColorTable(const string& filename, int band)
+void ImgWriterGdal::setColorTable(const std::string& filename, int band)
 {
   //todo: fool proof table in file (no checking currently done...)
-  ifstream ftable(filename.c_str(),ios::in);
-  string line;
+  std::ifstream ftable(filename.c_str(),std::ios::in);
+  std::string line;
 //   poCT=new GDALColorTable();
   GDALColorTable colorTable;
   short nline=0;
   while(getline(ftable,line)){
     ++nline;
-    istringstream ist(line);
+    std::istringstream ist(line);
     GDALColorEntry sEntry;
     short id;
     ist >> id >> sEntry.c1 >> sEntry.c2 >> sEntry.c3 >> sEntry.c4;
@@ -583,7 +597,7 @@ bool ImgWriterGdal::writeData(void* pdata, const GDALDataType& dataType, int ban
   //fetch raster band
   GDALRasterBand  *poBand;
   if(band>=nrOfBand()+1){
-    ostringstream s;
+    std::ostringstream s;
     s << "band (" << band << ") exceeds nrOfBand (" << nrOfBand() << ")";
     throw(s.str());
   }
diff --git a/src/imageclasses/ImgWriterGdal.h b/src/imageclasses/ImgWriterGdal.h
index d8c7d90..10fc7a5 100644
--- a/src/imageclasses/ImgWriterGdal.h
+++ b/src/imageclasses/ImgWriterGdal.h
@@ -3,7 +3,7 @@ ImgWriterGdal.h: class to write raster files using GDAL API library
 Copyright (C) 2008-2012 Pieter Kempeneers
 
 This file is part of pktools
-
+n
 pktools is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
@@ -28,32 +28,32 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "ImgReaderGdal.h"
 
 
-using namespace std;
-
 //--------------------------------------------------------------------------
 class ImgWriterGdal
 {
 public:
   ImgWriterGdal(void);
   ~ImgWriterGdal(void);
-  void open(const string& filename);
-  void open(const string& filename, const ImgReaderGdal& imgSrc, const vector<string>& options=vector<string>());
-  // void open(const string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const string& imageType="GTiff", const string& interleave="BAND", const string& compression="LZW", int magicX=1, int magicY=1);
-  void open(const string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const string& imageType, const vector<string>& options=vector<string>());
+  void open(const std::string& filename);
+  void open(const std::string& filename, const ImgReaderGdal& imgSrc, const std::vector<std::string>& options=std::vector<std::string>());
+  // void open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType="GTiff", const std::string& interleave="BAND", const std::string& compression="LZW", int magicX=1, int magicY=1);
+  void open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType, const std::vector<std::string>& options=std::vector<std::string>());
   void close(void);
   int nrOfCol(void) const { return m_ncol;};
   int nrOfRow(void) const { return m_nrow;};
   int nrOfBand(void) const { return m_nband;};
-  void setGeoTransform(double ulx, double uly, double deltaX, double deltaY, double rot1=0, double rot2=0);
   void copyGeoTransform(const ImgReaderGdal& imgSrc);
-  string setProjection(void);//set (and return) default projection ETSR-LAEA
-  void setProjection(const string& projection);
-  string setProjectionProj4(const string& projection);
-  void setImageDescription(const string& imageDescription){m_gds->SetMetadataItem( "TIFFTAG_IMAGEDESCRIPTION",imageDescription.c_str());};
+  std::string setProjection(void);//set (and return) default projection ETSR-LAEA
+  void setProjection(const std::string& projection);
+  std::string setProjectionProj4(const std::string& projection);
+  void setImageDescription(const std::string& imageDescription){m_gds->SetMetadataItem( "TIFFTAG_IMAGEDESCRIPTION",imageDescription.c_str());};
   CPLErr GDALSetNoDataValue(double noDataValue, int band=0) {getRasterBand(band)->SetNoDataValue(noDataValue);};
-  string getProjection(void) const;
-  string getGeoTransform() const;
-  void getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const;
+  std::string getProjection(void) const;
+  std::string getGeoTransform() const;
+  void getGeoTransform(double* gt) const;
+  void setGeoTransform(double* gt);
+  /* void setGeoTransform(double ulx, double uly, double deltaX, double deltaY, double rot1=0, double rot2=0); */
+  /* void getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const; */
   bool getBoundingBox(double& ulx, double& uly, double& lrx, double& lry) const;
   bool getCentrePos(double& x, double& y) const;
   bool covers(double x, double y) const;
@@ -62,41 +62,40 @@ public:
   bool image2geo(double i, double j, double& x, double& y) const;
   bool isGeoRef() const {return m_isGeoRef;};
   // void getMagicPixel(double& x, double& y) const {x=m_magic_x;y=m_magic_y;};
-  double getDeltaX(void) const {return m_delta_x;};
-  double getDeltaY(void) const{return m_delta_y;};
+  double getDeltaX(void) const {double gt[6];getGeoTransform(gt);return gt[1];};
+  double getDeltaY(void) const {double gt[6];getGeoTransform(gt);return -gt[5];};
   template<typename T> bool writeData(T& value, const GDALDataType& dataType, int col, int row, int band=0) const;
-  template<typename T> bool writeData(vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int row, int band=0) const;
-  template<typename T> bool writeData(vector<T>& buffer, const GDALDataType& dataType, int row, int band=0) const;
+  template<typename T> bool writeData(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int row, int band=0) const;
+  template<typename T> bool writeData(std::vector<T>& buffer, const GDALDataType& dataType, int row, int band=0) const;
   bool writeData(void* pdata, const GDALDataType& dataType, int band=0) const;
   template<typename T> bool writeDataBlock(Vector2d<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band=0) const;
-  // string getInterleave(){return m_interleave;};
-  // string getCompression(){return m_compression;};
+  // std::string getInterleave(){return m_interleave;};
+  // std::string getCompression(){return m_compression;};
   GDALDataType getDataType(int band=0) const;
   GDALRasterBand* getRasterBand(int band);
-  void setColorTable(const string& filename, int band=0);
+  void setColorTable(const std::string& filename, int band=0);
   void setColorTable(GDALColorTable* colorTable, int band=0);
   void setMetadata(char** metadata);
 
 protected:
-  void setCodec(const string& imageType);
+  void setCodec(const std::string& imageType);
   void setCodec(const ImgReaderGdal& ImgSrc);
 
-  string m_filename;
+  std::string m_filename;
   GDALDataset *m_gds;
   int m_ncol;
   int m_nrow;
   int m_nband;
   GDALDataType m_type;
-  double m_ulx;
-  double m_uly;
-  double m_delta_x;
-  double m_delta_y;
-  // double m_magic_x;
-  // double m_magic_y;
+  double m_gt[6];
+  /* double m_ulx; */
+  /* double m_uly; */
+  /* double m_delta_x; */
+  /* double m_delta_y; */
   bool m_isGeoRef;
-  // string m_interleave;
-  // string m_compression;
-  vector<string> m_options;
+  // std::string m_interleave;
+  // std::string m_compression;
+  std::vector<std::string> m_options;
 };
 
 template<typename T> bool ImgWriterGdal::writeData(T& value, const GDALDataType& dataType, int col, int row, int band) const
@@ -104,28 +103,28 @@ template<typename T> bool ImgWriterGdal::writeData(T& value, const GDALDataType&
   //fetch raster band
   GDALRasterBand  *poBand;
   if(band>=nrOfBand()+1){
-    ostringstream s;
+    std::ostringstream s;
     s << "band (" << band << ") exceeds nrOfBand (" << nrOfBand() << ")";
     throw(s.str());
   }
   poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
   if(col>=nrOfCol()){
-    ostringstream s;
+    std::ostringstream s;
     s << "col (" << col << ") exceeds nrOfCol (" << nrOfCol() << ")";
     throw(s.str());
   }
   if(col<0){
-    ostringstream s;
+    std::ostringstream s;
     s << "col (" << col << ") is negative";
     throw(s.str());
   }
   if(row>=nrOfRow()){
-    ostringstream s;
+    std::ostringstream s;
     s << "row (" << row << ") exceeds nrOfRow (" << nrOfRow() << ")";
     throw(s.str());
   }
   if(row<0){
-    ostringstream s;
+    std::ostringstream s;
     s << "row (" << row << ") is negative";
     throw(s.str());
   }
@@ -133,48 +132,48 @@ template<typename T> bool ImgWriterGdal::writeData(T& value, const GDALDataType&
   return true;
 }
 
-template<typename T> bool ImgWriterGdal::writeData(vector<T>& buffer, const GDALDataType& dataType, int minCol, int maxCol, int row, int band) const
+template<typename T> bool ImgWriterGdal::writeData(std::vector<T>& buffer, const GDALDataType& dataType, int minCol, int maxCol, int row, int band) const
 {
   //fetch raster band
   GDALRasterBand  *poBand;
   if(band>=nrOfBand()+1){
-    ostringstream s;
+    std::ostringstream s;
     s << "band (" << band << ") exceeds nrOfBand (" << nrOfBand() << ")";
     throw(s.str());
   }
   poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
   if(buffer.size()!=maxCol-minCol+1){
-    string errorstring="invalid buffer size";
+    std::string errorstring="invalid buffer size";
     throw(errorstring);
   }
   if(minCol>=nrOfCol()){
-    ostringstream s;
+    std::ostringstream s;
     s << "minCol (" << minCol << ") exceeds nrOfCol (" << nrOfCol() << ")";
     throw(s.str());
   }
   if(minCol<0){
-    ostringstream s;
+    std::ostringstream s;
     s << "mincol (" << minCol << ") is negative";
     throw(s.str());
   }
   if(maxCol>=nrOfCol()){
-    ostringstream s;
+    std::ostringstream s;
     s << "maxCol (" << maxCol << ") exceeds nrOfCol (" << nrOfCol() << ")";
     throw(s.str());
   }
   if(maxCol<minCol){
-    ostringstream s;
+    std::ostringstream s;
     s << "maxCol (" << maxCol << ") is less than minCol (" << minCol << ")";
     throw(s.str());
   }
 
   if(row>=nrOfRow()){
-    ostringstream s;
+    std::ostringstream s;
     s << "row (" << row << ") exceeds nrOfRow (" << nrOfRow() << ")";
     throw(s.str());
   }
   if(row<0){
-    ostringstream s;
+    std::ostringstream s;
     s << "row (" << row << ") is negative";
     throw(s.str());
   }
@@ -182,22 +181,22 @@ template<typename T> bool ImgWriterGdal::writeData(vector<T>& buffer, const GDAL
   return true;
 }
 
-template<typename T> bool ImgWriterGdal::writeData(vector<T>& buffer, const GDALDataType& dataType, int row, int band) const
+template<typename T> bool ImgWriterGdal::writeData(std::vector<T>& buffer, const GDALDataType& dataType, int row, int band) const
 {
   //fetch raster band
   GDALRasterBand  *poBand;
   if(band>=nrOfBand()+1){
-    ostringstream s;
+    std::ostringstream s;
     s << "band (" << band << ") exceeds nrOfBand (" << nrOfBand() << ")";
     throw(s.str());
   }
   poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
   if(buffer.size()!=nrOfCol()){
-    string errorstring="invalid buffer size";
+    std::string errorstring="invalid buffer size";
     throw(errorstring);
   }
   if(row>=nrOfRow()){
-    ostringstream s;
+    std::ostringstream s;
     s << "row (" << row << ") exceeds nrOfRow (" << nrOfRow() << ")";
     throw(s.str());
   }
@@ -210,7 +209,7 @@ template<typename T> bool ImgWriterGdal::writeDataBlock(Vector2d<T>& buffer, con
   //fetch raster band
   GDALRasterBand  *poBand;
   if(band>=nrOfBand()+1){
-    ostringstream s;
+    std::ostringstream s;
     s << "band (" << band << ") exceeds nrOfBand (" << nrOfBand() << ")";
     throw(s.str());
   }

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