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 e93a8388d46cc553b26f827aded82d1f9bd89e52 Author: Pieter Kempeneers <kempe...@gmail.com> Date: Fri Jan 25 17:39:34 2013 +0100 see Changelog --- ChangeLog | 11 +- src/algorithms/Filter2d.cc | 76 ++++++++--- src/algorithms/Filter2d.h | 281 +++++++++++++++++++++++++++++++++++++-- src/apps/pkclassify_svm.cc | 1 + src/apps/pkcrop.cc | 4 +- src/apps/pkfilter.cc | 24 ++-- src/apps/pklas2img.cc | 26 ++-- src/fileclasses/FileReaderLas.cc | 39 ++++++ src/fileclasses/FileReaderLas.h | 12 ++ 9 files changed, 422 insertions(+), 52 deletions(-) diff --git a/ChangeLog b/ChangeLog index 93b1c13..f907a49 100755 --- a/ChangeLog +++ b/ChangeLog @@ -42,14 +42,20 @@ version 2.4 --enable-nlopt (when NLOPT is installed, needed for pksensormodel and pkgetchandelier) - OptFactory.h (added) factory class for nlopt::opt (selecting algorithm via string) + - Filter2d + support filtering of matrix (doit with Vector2d arguments) (to create DTM according to P. Bunting) + - FileReaderLas + support class filters - ImgReaderGdal.cc in addition to internal setNoData member variable, also support GDALSetNoData - ImgWriterGdal.cc in addition to internal setNoData member variable, also support GDALSetNoData - ImgWriterOGR.h renamed ascii2shape to ascii2ogr, support csv file in ascii2ogr - - ImgWriterOGR.cc + - ImgWriterOgr.cc renamed ascii2shape to ascii2ogr, support csv file in ascii2ogr + - pkcrop + changed default option for color table to empty vector - pkinfo support of nodata value when calculating histogram and image statistics options min and max are now set with -min (--min) and -max (--max) @@ -58,6 +64,7 @@ version 2.4 - pkfilter bug solved in update of progress bar support of standard deviation + default empty classes - pkgetmask options min and max are now set with -min (--min) and -max (--max) - pkstatogr @@ -81,3 +88,5 @@ version 2.4 tool to model pushbroom sensor (with optimization of boresight angles using NLOPT) - pkascii2ogr support csv input file + - pklas2img + support DTM according to Bunting (slightly adapted version of Chang2003, including a median filter after morphological operator) diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc index 12dd080..e46af52 100644 --- a/src/algorithms/Filter2d.cc +++ b/src/algorithms/Filter2d.cc @@ -503,53 +503,85 @@ void Filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, } switch(method){ case(MEDIAN): - outBuffer[x/down]=hist.median(windowBuffer); + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.median(windowBuffer); break; case(VAR):{ - outBuffer[x/down]=hist.var(windowBuffer); + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.var(windowBuffer); break; } case(STDEV):{ - outBuffer[x/down]=sqrt(hist.var(windowBuffer)); + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=sqrt(hist.var(windowBuffer)); break; } case(MEAN):{ - outBuffer[x/down]=hist.mean(windowBuffer); + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.mean(windowBuffer); break; } case(MIN):{ - outBuffer[x/down]=hist.min(windowBuffer); + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.min(windowBuffer); break; } case(ISMIN):{ - outBuffer[x/down]=(hist.min(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0; + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=(hist.min(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0; break; } case(MINMAX):{ double min=0; double max=0; - hist.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max); - if(min!=max) - outBuffer[x/down]=0; - else - outBuffer[x/down]=windowBuffer[dimX*dimY/2];//centre pixels + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else{ + hist.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max); + if(min!=max) + outBuffer[x/down]=0; + else + outBuffer[x/down]=windowBuffer[dimX*dimY/2];//centre pixels + } break; } case(MAX):{ - outBuffer[x/down]=hist.max(windowBuffer); + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.max(windowBuffer); break; } case(ISMAX):{ - outBuffer[x/down]=(hist.max(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0; + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=(hist.max(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0; break; } case(ORDER):{ - double lbound=0; - double ubound=dimX*dimY; - double theMin=hist.min(windowBuffer); - double theMax=hist.max(windowBuffer); - double scale=(ubound-lbound)/(theMax-theMin); - outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[dimX*dimY/2]-theMin)+lbound); + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else{ + double lbound=0; + double ubound=dimX*dimY; + double theMin=hist.min(windowBuffer); + double theMax=hist.max(windowBuffer); + double scale=(ubound-lbound)/(theMax-theMin); + outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[dimX*dimY/2]-theMin)+lbound); + } break; } case(SUM):{ @@ -582,7 +614,7 @@ void Filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, outBuffer[x/down]+=100.0*occurrence[*(vit++)]/windowBuffer.size(); } else - outBuffer[x/down]=0; + outBuffer[x/down]=m_noValue; break; } case(MAJORITY):{ @@ -598,7 +630,7 @@ void Filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, outBuffer[x/down]=inBuffer[dimY/2][x]; } else - outBuffer[x/down]=0; + outBuffer[x/down]=m_noValue; break; } case(THRESHOLD):{ @@ -611,7 +643,7 @@ void Filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, } } else - outBuffer[x/down]=0; + outBuffer[x/down]=m_noValue; break; } case(MIXED):{ diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h index 3a76df0..954f1fc 100644 --- a/src/algorithms/Filter2d.h +++ b/src/algorithms/Filter2d.h @@ -69,6 +69,7 @@ public: /* void homogeneousSpatial(const string& inputFilename, const string& outputFilename, int dim, bool disc=false, int noValue=0); */ void doit(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dim, short down=2, bool disc=false); void doit(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dimX, int dimY, short down=1, bool disc=false); + template<class T1, class T2> void doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, int method, int dimX, int dimY, short down=1, bool disc=false); void median(const string& inputFilename, const string& outputFilename, int dim, bool disc=false); void var(const string& inputFilename, const string& outputFilename, int dim, bool disc=false); void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dimX, int dimY, bool disc=false, double angle=-190); @@ -101,13 +102,59 @@ private: filter(inputVector,outputVector); } - template<class T1, class T2> void Filter2d::filter(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector) + template<class T1, class T2> void Filter2d::filter(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector) { - outputVector.resize(inputVector.size()); - int dimX=m_taps[0].size();//horizontal!!! - int dimY=m_taps.size();//vertical!!! + outputVector.resize(inputVector.size()); + int dimX=m_taps[0].size();//horizontal!!! + int dimY=m_taps.size();//vertical!!! + Vector2d<T1> inBuffer(dimY); + vector<T2> outBuffer(inputVector[0].size()); + //initialize last half of inBuffer + int indexI=0; + int indexJ=0; + for(int y=0;y<dimY;++y){ + if(y<dimY/2) + continue;//skip first half + inBuffer[y]=inputVector[indexJ++]; + } + for(int y=0;y<inputVector.size();++y){ + if(y){//inBuffer already initialized for y=0 + //erase first line from inBuffer + inBuffer.erase(inBuffer.begin()); + //read extra line and push back to inBuffer if not out of bounds + if(y+dimY/2<inputVector.size()) + inBuffer.push_back(inputVector[y+dimY/2]); + } + for(int x=0;x<inputVector[0].size();++x){ + outBuffer[x]=0; + for(int j=-dimY/2;j<(dimY+1)/2;++j){ + for(int i=-dimX/2;i<(dimX+1)/2;++i){ + indexI=x+i; + indexJ=dimY/2+j; + //check if out of bounds + if(x<dimX/2) + indexI=x+abs(i); + else if(x>=inputVector[0].size()-dimX/2) + indexI=x-abs(i); + if(y<dimY/2) + indexJ=dimY/2+abs(j); + else if(y>=inputVector.size()-dimY/2) + indexJ=dimY/2-abs(j); + outBuffer[x]+=(m_taps[dimY/2+j][dimX/2+i]*inBuffer[indexJ][indexI]); + } + } + } + //copy outBuffer to outputVector + outputVector[y]=outBuffer; + } + } + +template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, int method, int dimX, int dimY, short down, bool disc) +{ + Histogram hist; + outputVector.resize((inputVector.size()+down-1)/down); Vector2d<T1> inBuffer(dimY); - vector<T2> outBuffer(inputVector[0].size()); + vector<T2> outBuffer((inputVector[0].size()+down-1)/down); //initialize last half of inBuffer int indexI=0; int indexJ=0; @@ -116,6 +163,11 @@ private: continue;//skip first half inBuffer[y]=inputVector[indexJ++]; } + const char* pszMessage; + void* pProgressArg=NULL; + GDALProgressFunc pfnProgress=GDALTermProgress; + double progress=0; + pfnProgress(progress,pszMessage,pProgressArg); for(int y=0;y<inputVector.size();++y){ if(y){//inBuffer already initialized for y=0 //erase first line from inBuffer @@ -123,9 +175,22 @@ private: //read extra line and push back to inBuffer if not out of bounds if(y+dimY/2<inputVector.size()) inBuffer.push_back(inputVector[y+dimY/2]); + else{ + int over=y+dimY/2-inputVector.size(); + int index=(inBuffer.size()-1)-over; + assert(index>=0); + assert(index<inBuffer.size()); + inBuffer.push_back(inBuffer[index]); + } } + if((y+1+down/2)%down) + continue; for(int x=0;x<inputVector[0].size();++x){ - outBuffer[x]=0; + if((x+1+down/2)%down) + continue; + outBuffer[x/down]=0; + vector<double> windowBuffer; + map<int,int> occurrence; for(int j=-dimY/2;j<(dimY+1)/2;++j){ for(int i=-dimX/2;i<(dimX+1)/2;++i){ indexI=x+i; @@ -139,12 +204,210 @@ private: indexJ=dimY/2+abs(j); else if(y>=inputVector.size()-dimY/2) indexJ=dimY/2-abs(j); - outBuffer[x]+=(m_taps[dimY/2+j][dimX/2+i]*inBuffer[indexJ][indexI]); - } + bool masked=false; + for(int imask=0;imask<m_mask.size();++imask){ + if(inBuffer[indexJ][indexI]==m_mask[imask]){ + masked=true; + break; + } + } + if(!masked){ + vector<short>::const_iterator vit=m_class.begin(); + //todo: test if this works (only add occurrence if within defined classes)! + if(!m_class.size()) + ++occurrence[inBuffer[indexJ][indexI]]; + else{ + while(vit!=m_class.end()){ + if(inBuffer[indexJ][indexI]==*(vit++)) + ++occurrence[inBuffer[indexJ][indexI]]; + } + } + windowBuffer.push_back(inBuffer[indexJ][indexI]); + } + } + } + switch(method){ + case(MEDIAN): + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.median(windowBuffer); + break; + case(VAR):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.var(windowBuffer); + break; + } + case(STDEV):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=sqrt(hist.var(windowBuffer)); + break; + } + case(MEAN):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.mean(windowBuffer); + break; + } + case(MIN):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.min(windowBuffer); + break; + } + case(ISMIN):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=(hist.min(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0; + break; + } + case(MINMAX):{ + double min=0; + double max=0; + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else{ + hist.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max); + if(min!=max) + outBuffer[x/down]=0; + else + outBuffer[x/down]=windowBuffer[dimX*dimY/2];//centre pixels + } + break; + } + case(MAX):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=hist.max(windowBuffer); + break; + } + case(ISMAX):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else + outBuffer[x/down]=(hist.max(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0; + break; + } + case(ORDER):{ + if(windowBuffer.empty()) + outBuffer[x/down]=m_noValue; + else{ + double lbound=0; + double ubound=dimX*dimY; + double theMin=hist.min(windowBuffer); + double theMax=hist.max(windowBuffer); + double scale=(ubound-lbound)/(theMax-theMin); + outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[dimX*dimY/2]-theMin)+lbound); + } + break; + } + case(SUM):{ + outBuffer[x/down]=hist.sum(windowBuffer); + break; + } + case(HOMOG): + if(occurrence.size()==1)//all values in window must be the same + outBuffer[x/down]=inBuffer[dimY/2][x]; + else//favorize original value in case of ties + outBuffer[x/down]=m_noValue; + break; + case(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/2][x]) + outBuffer[x/down]=1; + else if(*wit==inBuffer[dimY/2][x]){//todo:wit mag niet central pixel zijn + outBuffer[x/down]=m_noValue; + break; + } + } + break; + } + case(DENSITY):{ + if(windowBuffer.size()){ + 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; + break; + } + case(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){ + if(mit->second>maxit->second) + maxit=mit; + } + if(occurrence[inBuffer[dimY/2][x]]<maxit->second)// + outBuffer[x/down]=maxit->first; + else//favorize original value in case of ties + outBuffer[x/down]=inBuffer[dimY/2][x]; + } + else + outBuffer[x/down]=m_noValue; + break; + } + case(THRESHOLD):{ + assert(m_class.size()==m_threshold.size()); + if(windowBuffer.size()){ + outBuffer[x/down]=inBuffer[dimY/2][x];//initialize with original value (in case thresholds not met) + for(int iclass=0;iclass<m_class.size();++iclass){ + if(100.0*(occurrence[m_class[iclass]])/windowBuffer.size()>m_threshold[iclass]) + outBuffer[x/down]=m_class[iclass]; + } + } + else + outBuffer[x/down]=m_noValue; + break; + } + case(MIXED):{ + enum Type { BF=11, CF=12, MF=13, NF=20, W=30 }; + double nBF=occurrence[BF]; + double nCF=occurrence[CF]; + double nMF=occurrence[MF]; + double nNF=occurrence[NF]; + double nW=occurrence[W]; + if(windowBuffer.size()){ + if((nBF+nCF+nMF)&&(nBF+nCF+nMF>=nNF+nW)){//forest + if(nBF/(nBF+nCF)>=0.75) + outBuffer[x/down]=BF; + else if(nCF/(nBF+nCF)>=0.75) + outBuffer[x/down]=CF; + else + outBuffer[x/down]=MF; + } + else{//non-forest + if(nW&&(nW>=nNF)) + outBuffer[x/down]=W; + else + outBuffer[x/down]=NF; + } + } + else + outBuffer[x/down]=inBuffer[indexJ][indexI]; + break; + } + default: + break; } } + progress=(1.0+y/down); + progress+=(outputVector.size()); + progress/=outputVector.size(); + pfnProgress(progress,pszMessage,pProgressArg); //copy outBuffer to outputVector - outputVector[y]=outBuffer; + outputVector[y/down]=outBuffer; } } diff --git a/src/apps/pkclassify_svm.cc b/src/apps/pkclassify_svm.cc index 3c86021..94139de 100644 --- a/src/apps/pkclassify_svm.cc +++ b/src/apps/pkclassify_svm.cc @@ -17,6 +17,7 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with pktools. If not, see <http://www.gnu.org/licenses/>. ***********************************************************************/ + #include <vector> #include <map> #include <algorithm> diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc index f17dc28..263f0ca 100644 --- a/src/apps/pkcrop.cc +++ b/src/apps/pkcrop.cc @@ -60,7 +60,7 @@ int main(int argc, char *argv[]) Optionpk<double> offset_opt("off", "offset", "output=scale*input+offset", 0); Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image",""); Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image"); - Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)", ""); + Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)"); Optionpk<string> option_opt("co", "co", "options: NAME=VALUE [-co COMPRESS=LZW] [-co INTERLEAVE=BAND]"); Optionpk<short> flag_opt("f", "flag", "Flag value to put in image if out of bounds.", 0); Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbour, bilinear: bi-linear interpolation).", "near"); @@ -327,7 +327,7 @@ int main(int argc, char *argv[]) } else if(imgReader.isGeoRef()) imgWriter.setProjection(imgReader.getProjection()); - if(colorTable_opt[0]!=""){ + if(colorTable_opt.size()){ if(verbose_opt[0]) cout << "set colortable " << colorTable_opt[0] << endl; assert(imgWriter.getDataType()==GDT_Byte); diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc index f77fb21..7beadb0 100644 --- a/src/apps/pkfilter.cc +++ b/src/apps/pkfilter.cc @@ -50,8 +50,8 @@ int main(int argc,char **argv) { Optionpk<bool> license_opt("lic","license","show license information",false); Optionpk<bool> help_opt("h","help","shows this help info",false); Optionpk<bool> todo_opt("\0","todo","",false); - Optionpk<std::string> input_opt("i","input","input image file",""); - Optionpk<std::string> output_opt("o", "output", "Output image file", ""); + Optionpk<std::string> input_opt("i","input","input image file"); + Optionpk<std::string> output_opt("o", "output", "Output image file"); Optionpk<bool> disc_opt("c", "circular", "circular disc kernel for dilation and erosion", false); Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation."); Optionpk<int> function_opt("f", "filter", "filter function (0: median, 1: variance, 2: min, 3: max, 4: sum, 5: mean, 6: minmax, 7: dilation, 8: erosion, 9: closing, 10: opening, 11: spatially homogeneous (central pixel must be identical to all other pixels within window), 12: SobelX edge detection in X, 13: SobelY edge detection in Y, 14: SobelXY, -14: SobelYX, 15: smooth, 16: density, 17: majority voting (only for classes), 18: forest aggregation (mixed), 19: smooth no data (mask) val [...] @@ -60,12 +60,12 @@ int main(int argc,char **argv) { 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<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)", 0); - Optionpk<short> mask_opt("\0", "mask", "mask value(s) "); + Optionpk<short> mask_opt("m", "mask", "mask value(s) "); Optionpk<std::string> tap_opt("tap", "tap", "text file containing taps used for spatial filtering (from ul to lr). Use dimX and dimY to specify tap dimensions in x and y. Leave empty for not using taps"); Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering"); - Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", ""); + Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image",""); Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image"); - Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)", ""); + Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)"); Optionpk<std::string> option_opt("co", "co", "options: NAME=VALUE [-co COMPRESS=LZW] [-co INTERLEAVE=BAND]"); Optionpk<short> down_opt("d", "down", "down sampling factor. Use value 1 for no downsampling)", 1); Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0); @@ -111,6 +111,7 @@ int main(int argc,char **argv) { ImgReaderGdal input; ImgWriterGdal output; + assert(input_opt.size()); input.open(input_opt[0]); // output.open(output_opt[0],input); GDALDataType theType=GDT_Unknown; @@ -140,6 +141,7 @@ int main(int argc,char **argv) { option_opt.push_back(theInterleave); } try{ + assert(output_opt.size()); output.open(output_opt[0],(input.nrOfCol()+down_opt[0]-1)/down_opt[0],(input.nrOfRow()+down_opt[0]-1)/down_opt[0],input.nrOfBand(),theType,imageType,option_opt); } catch(string errorstring){ @@ -152,7 +154,13 @@ int main(int argc,char **argv) { input.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2); output.setGeoTransform(ulx,uly,deltaX*down_opt[0],deltaY*down_opt[0],rot1,rot2); } - if(input.getColorTable()!=NULL) + if(colorTable_opt.size()){ + if(verbose_opt[0]) + cout << "set colortable " << colorTable_opt[0] << endl; + assert(output.getDataType()==GDT_Byte); + output.setColorTable(colorTable_opt[0]); + } + else if(input.getColorTable()!=NULL) output.setColorTable(input.getColorTable()); // if(input.isGeoRef()){ @@ -162,7 +170,7 @@ int main(int argc,char **argv) { Filter2d::Filter2d filter2d; Filter filter1d; - if(verbose_opt[0]) + if(class_opt.size()&&verbose_opt[0]) std::cout<< "class values: "; for(int iclass=0;iclass<class_opt.size();++iclass){ if(!dimZ_opt.size()) @@ -211,7 +219,7 @@ int main(int argc,char **argv) { filter1d.doit(input,output,down_opt[0]); } else{ - if(colorTable_opt[0]!="") + if(colorTable_opt.size()) output.setColorTable(colorTable_opt[0]); switch(function_opt[0]){ case(Filter2d::MEDIAN): diff --git a/src/apps/pklas2img.cc b/src/apps/pklas2img.cc index fe26312..4233052 100644 --- a/src/apps/pklas2img.cc +++ b/src/apps/pklas2img.cc @@ -50,9 +50,10 @@ int main(int argc,char **argv) { Optionpk<double> hThreshold_opt("ht", "maxHeight", "initial and maximum height threshold for progressive morphological filtering (e.g., -ht 0.2 -ht 2.5)", 0.2); Optionpk<short> maxIter_opt("\0", "maxIter", "Maximum number of iterations in post filter (default is 100)", 100.0); Optionpk<short> nbin_opt("nb", "nbin", "Number of percentile bins for calculating profile (=number of output bands) (default is 10)", 10.0); - Optionpk<unsigned short> returns_opt("r", "returns", "number(s) of returns to include (use -r -1 for last return only). Default is 0 (include all returns)", 0); - Optionpk<string> composite_opt("c", "composite", "composite for multiple points in cell (min, max, median, mean, sum, first, last, profile, number (point density)). Default is last (overwrite cells with latest point", "last"); - Optionpk<string> filter_opt("fir", "filter", "filter las points (single,multiple,all). Default is all", "all"); + Optionpk<unsigned short> returns_opt("r", "returns", "number(s) of returns to include"); + Optionpk<unsigned short> classes_opt("c", "classes", "classes to keep: 0 (created, never classified), 1 (unclassified), 2 (ground), 3 (low vegetation), 4 (medium vegetation), 5 (high vegetation), 6 (building), 7 (low point, noise), 8 (model key-point), 9 (water), 10 (reserved), 11 (reserved), 12 (overlap)"); + Optionpk<string> composite_opt("comp", "comp", "composite for multiple points in cell (min, max, median, mean, sum, first, last, profile, number (point density)). Default is last (overwrite cells with latest point", "last"); + Optionpk<string> filter_opt("fir", "filter", "filter las points (last,single,multiple,all). Default is all", "all"); Optionpk<string> postFilter_opt("pf", "pfilter", "post processing filter (etew_min,promorph (progressive morphological filter),bunting (adapted promorph),open,close,none) . Default is none", "none"); Optionpk<short> dimx_opt("\0", "dimX", "Dimension X of postFilter (default is 3)", 3); Optionpk<short> dimy_opt("\0", "dimY", "Dimension Y of postFilter (default is 3)", 3); @@ -94,6 +95,7 @@ int main(int argc,char **argv) { maxIter_opt.retrieveOption(argc,argv); nbin_opt.retrieveOption(argc,argv); returns_opt.retrieveOption(argc,argv); + classes_opt.retrieveOption(argc,argv); composite_opt.retrieveOption(argc,argv); filter_opt.retrieveOption(argc,argv); postFilter_opt.retrieveOption(argc,argv); @@ -256,8 +258,10 @@ int main(int argc,char **argv) { //set bounding filter // lasReader.addBoundsFilter(minULX,maxULY,maxLRX,minLRY); //set returns filter - if(returns_opt[0]) + if(returns_opt.size()) lasReader.addReturnsFilter(returns_opt); + if(classes_opt.size()) + lasReader.addClassFilter(classes_opt); lasReader.setFilters(); if(attribute_opt[0]!="z"){ @@ -295,6 +299,10 @@ int main(int argc,char **argv) { continue; if((filter_opt[0]=="multiple")&&(thePoint.GetNumberOfReturns()<2)) continue; + if(filter_opt[0]=="last"){ + if(thePoint.GetReturnNumber()!=thePoint.GetNumberOfReturns()) + continue; + } double dcol,drow; outputWriter.geo2image(thePoint.GetX(),thePoint.GetY(),dcol,drow); int icol=static_cast<int>(dcol); @@ -454,12 +462,10 @@ int main(int argc,char **argv) { try{ theFilter.morphology(outputData,currentOutput,Filter2d::ERODE,dimx,dimy,disc_opt[0],maxSlope_opt[0]); theFilter.morphology(currentOutput,outputData,Filter2d::DILATE,dimx,dimy,disc_opt[0],maxSlope_opt[0]); - // if(postFilter_opt[0]=="bunting"){//todo: implement doit in Filter2d on Vector2d - // theFilter.doit(outputData,currentOutput,Filter2d::MEDIAN,dimx,dimy,1,disc_opt[0]); - // filter2d.doit(input,output,Filter2d::MEDIAN,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]); - - // outputData=currentOutput; - // } + if(postFilter_opt[0]=="bunting"){//todo: implement doit in Filter2d on Vector2d + theFilter.doit(outputData,currentOutput,Filter2d::MEDIAN,dimx,dimy,1,disc_opt[0]); + outputData=currentOutput; + } } catch(std::string errorString){ cout << errorString << endl; diff --git a/src/fileclasses/FileReaderLas.cc b/src/fileclasses/FileReaderLas.cc index 25f1d42..aceaac1 100644 --- a/src/fileclasses/FileReaderLas.cc +++ b/src/fileclasses/FileReaderLas.cc @@ -22,6 +22,30 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>. #include <fstream> #include "FileReaderLas.h" //--------------------------------------------------------------------------- +LastReturnFilter::LastReturnFilter( ) : liblas::FilterI(eInclusion) {} + +bool LastReturnFilter::filter(const liblas::Point& p) +{ + + // If the GetReturnNumber equals the GetNumberOfReturns, + // we're a last return + + bool output = false; + if (p.GetReturnNumber() == p.GetNumberOfReturns()) + { + output = true; + } + + // If the type is switched to eExclusion, we'll throw out all last returns. + if (GetType() == eExclusion && output == true) + { + output = false; + } else { + output = true; + } + return output; +} + FileReaderLas::FileReaderLas(void) { m_reader=NULL; @@ -138,3 +162,18 @@ void FileReaderLas::addReturnsFilter(std::vector<unsigned short> const& returns) m_filters.push_back(liblas::FilterPtr(return_filter)); } +void FileReaderLas::addClassFilter(std::vector<unsigned short> const& classes){ + + std::vector<liblas::FilterPtr> filters; + std::vector<liblas::Classification> theClasses; + for(int iclass=0;iclass<classes.size();++iclass){ + liblas::Classification aClass(classes[iclass]); + theClasses.push_back(aClass); + } + liblas::FilterPtr class_filter = liblas::FilterPtr(new liblas::ClassificationFilter(theClasses)); + // eInclusion means to keep the classes that match. eExclusion would + // throw out those that matched + class_filter->SetType(liblas::FilterI::eInclusion); + m_filters.push_back(class_filter); +} + diff --git a/src/fileclasses/FileReaderLas.h b/src/fileclasses/FileReaderLas.h index 93023dd..8435d01 100644 --- a/src/fileclasses/FileReaderLas.h +++ b/src/fileclasses/FileReaderLas.h @@ -25,6 +25,17 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>. #include "liblas/liblas.hpp" //-------------------------------------------------------------------------- +class LastReturnFilter: public liblas::FilterI +{ +public: + LastReturnFilter(); + bool filter(const liblas::Point& point); + +private: + LastReturnFilter(LastReturnFilter const& other); + LastReturnFilter& operator=(LastReturnFilter const& rhs); +}; + class FileReaderLas { public: @@ -47,6 +58,7 @@ public: liblas::Point const& readPointAt(std::size_t n){m_reader->ReadPointAt(n);return m_reader->GetPoint();}; // void addBoundsFilter(double ulx, double uly, double lrx, double lry); void addReturnsFilter(std::vector<unsigned short> const& returns); + void addClassFilter(std::vector<unsigned short> const& classes); void setFilters(const std::vector<liblas::FilterPtr>& filters){m_filters=filters;setFilters();}; void setFilters(){m_reader->SetFilters(m_filters);}; protected: -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git _______________________________________________ Pkg-grass-devel mailing list Pkg-grass-devel@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel