[pktools] 81/375: filter and pkclassify_svm at work
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:01 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 d37f521248a9073a2243f94d5e7a72750e34af89
Author: user <user at osgeolive.(none)>
Date: Tue Mar 26 17:59:08 2013 +0100
filter and pkclassify_svm at work
---
src/algorithms/Filter2d.cc | 28 ++++++++++----
src/algorithms/Filter2d.h | 1 +
src/apps/pkclassify_svm.cc | 94 +++++++++++++++++++++++++---------------------
src/apps/pkfilter.cc | 29 +++++++++++++-
4 files changed, 101 insertions(+), 51 deletions(-)
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index 562bf88..3ed1033 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -692,11 +692,20 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
}
}
-void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, double beta, bool eightConnectivity, short down, bool verbose)
+void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, double beta, bool eightConnectivity, short down, bool verbose){
+ assert(m_class.size()>1);
+ Vector2d<double> fullBeta(m_class.size(),m_class.size());
+ for(int iclass1=0;iclass1<m_class.size();++iclass1)
+ for(int iclass2=0;iclass2<m_class.size();++iclass2)
+ fullBeta[iclass1][iclass2]=beta;
+ mrf(input,output,dimX,dimY,fullBeta,eightConnectivity,down,verbose);
+}
+
+//beta[classTo][classFrom]
+void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, Vector2d<double> beta, bool eightConnectivity, short down, bool verbose)
{
assert(dimX);
assert(dimY);
- assert(m_class.size()>1);
const char* pszMessage;
void* pProgressArg=NULL;
GDALProgressFunc pfnProgress=GDALTermProgress;
@@ -706,6 +715,8 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
Vector2d<double> outBuffer(m_class.size(),(input.nrOfCol()+down-1)/down);
assert(input.nrOfBand()==1);
assert(output.nrOfBand()==m_class.size());
+ assert(m_class.size()>1);
+ assert(beta.size()==m_class.size());
//initialize last half of inBuffer
int indexI=0;
int indexJ=0;
@@ -785,11 +796,14 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input, ImgWriterGdal& output,
}
}
}
- for(int iclass=0;iclass<m_class.size();++iclass){
- if(eightConnectivity)//todo: normalize by 1/z ?
- outBuffer[iclass][x/down]=exp(-beta*(8-potential[iclass]));
- else
- outBuffer[iclass][x/down]=exp(-beta*(4-potential[iclass]));
+ for(int iclass1=0;iclass1<m_class.size();++iclass1){
+ assert(beta[iclass1].size()==m_class.size());
+ for(int iclass2=0;iclass2<m_class.size();++iclass2){
+ if(eightConnectivity)
+ outBuffer[iclass1][x/down]=exp(-beta[iclass1][iclass2]*(dimX*dimY-1-potential[iclass1]))/(exp(-beta[iclass1][iclass2]*(dimX*dimY-1-potential[iclass1]))+exp(-beta[iclass1][iclass2]*(potential[iclass1])));
+ else
+ outBuffer[iclass1][x/down]=exp(-beta[iclass1][iclass2]*(dimX+dimY-1-potential[iclass1]))/(exp(-beta[iclass1][iclass2]*(dimX+dimY-1-potential[iclass1]))+exp(-beta[iclass1][iclass2]*(potential[iclass1])));
+ }
}
}
progress=(1.0+y/down)/output.nrOfRow();
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index cfd05e7..0f9d592 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -77,6 +77,7 @@ public:
void doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=2, bool disc=false);
void doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, short down=1, bool disc=false);
void mrf(const ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, double beta, bool eightConnectivity=true, short down=1, bool verbose=false);
+ void mrf(const ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, Vector2d<double> beta, bool eightConnectivity=true, short down=1, bool verbose=false);
template<class T1, class T2> void doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, const std::string& method, int dimX, int dimY, short down=1, bool disc=false);
void median(const string& inputFilename, const string& outputFilename, int dim, bool disc=false);
void var(const string& inputFilename, const string& outputFilename, int dim, bool disc=false);
diff --git a/src/apps/pkclassify_svm.cc b/src/apps/pkclassify_svm.cc
index 893f9b6..bd265dc 100644
--- a/src/apps/pkclassify_svm.cc
+++ b/src/apps/pkclassify_svm.cc
@@ -563,12 +563,11 @@ int main(int argc, char *argv[])
if(classname_opt.empty()){
std::cerr << "Warning: no class name and value pair provided for all " << nclass << " classes, using string2type<int> instead!" << std::endl;
for(int iclass=0;iclass<nclass;++iclass){
- if(verbose_opt[0]>0)
+ if(verbose_opt[0])
std::cout << iclass << " " << cm.getClass(iclass) << " -> " << string2type<short>(cm.getClass(iclass)) << std::endl;
classValueMap[cm.getClass(iclass)]=string2type<short>(cm.getClass(iclass));
}
}
-
ImgReaderGdal testImage;
try{
if(verbose_opt[0]>=1)
@@ -672,8 +671,7 @@ int main(int argc, char *argv[])
if(priorimg_opt.size())
linePrior.resize(nclass,ncol);//prior prob for each class
Vector2d<float> hpixel(ncol);
- vector<float> prOut(nclass);//posterior prob for each (internal) class
- Vector2d<float> probOut(nclass,ncol);//posterior prob for each output class
+ Vector2d<float> probOut(nclass,ncol);//posterior prob for each (internal) class
vector<float> entropy(ncol);
Vector2d<char> classBag;//classified line for writing to image file
if(classBag_opt.size())
@@ -717,7 +715,7 @@ int main(int argc, char *argv[])
exit(3);
}
assert(nband==hpixel[0].size());
- if(verbose_opt[0]==2)
+ if(verbose_opt[0]>1)
std::cout << "used bands: " << nband << std::endl;
//read mask
if(!lineMask.empty()){
@@ -736,11 +734,14 @@ int main(int argc, char *argv[])
//read prior
if(priorimg_opt.size()){
try{
- for(int iclass=0;iclass<nclass;++iclass)
+ for(short iclass=0;iclass<nclass;++iclass){
+ if(verbose_opt.size()>1)
+ std::cout << "Reading " << priorimg_opt[0] << " band " << iclass << " line " << iline << std::endl;
priorReader.readData(linePrior[iclass],GDT_Float32,iline,iclass);
+ }
}
catch(string theError){
- cerr << "Error reading " << priorimg_opt[0] << ": " << theError << std::endl;
+ std::cerr << "Error reading " << priorimg_opt[0] << ": " << theError << std::endl;
exit(3);
}
catch(...){
@@ -817,7 +818,9 @@ int main(int argc, char *argv[])
continue;//next column
}
for(short iclass=0;iclass<nclass;++iclass)
- prOut[iclass]=0;
+ probOut[iclass][icol]=0;
+ if(verbose_opt[0]>1)
+ std::cout << "begin classification " << std::endl;
//----------------------------------- classification -------------------
for(int ibag=0;ibag<nbag;++ibag){
//calculate image features
@@ -871,25 +874,32 @@ int main(int argc, char *argv[])
maxP=0;
classBag[ibag][icol]=0;
}
+ double normPrior=0;
+ if(priorimg_opt.size()){
+ for(short iclass=0;iclass<nclass;++iclass)
+ normPrior+=linePrior[iclass][icol];
+ }
for(short iclass=0;iclass<nclass;++iclass){
if(priorimg_opt.size())
- priors[iclass]=linePrior[classValueMap[cm.getClass(iclass)]][icol];//todo: check if correct for all cases... (automatic classValueMap and manual input for names and values)
- // priors[iclass]=linePrior[iclass][icol];
+ priors[iclass]=linePrior[iclass][icol]/normPrior;//todo: check if correct for all cases... (automatic classValueMap and manual input for names and values)
switch(comb_opt[0]){
default:
case(0)://sum rule
- // prOut[iclass][icol]+=prValues[iclass]+static_cast<float>(1.0-nbag)/nbag*priors[iclass];//add probabilities for each bag
- prOut[iclass]+=result[iclass]+static_cast<float>(1.0-nbag)/nbag*priors[iclass];//add probabilities for each bag
+ // probOut[iclass][icol]+=prValues[iclass]+static_cast<float>(1.0-nbag)/nbag*priors[iclass];//add probabilities for each bag
+ probOut[iclass][icol]+=result[iclass]*priors[iclass];//add probabilities for each bag
+ // probOut[iclass][icol]+=result[iclass]+static_cast<float>(1.0-nbag)/nbag*priors[iclass];//add probabilities for each bag
break;
case(1)://product rule
- // prOut[iclass][icol]*=pow(priors[iclass],static_cast<float>(1.0-nbag)/nbag)*prValues[iclass];//add probabilities for each bag
- prOut[iclass]*=pow(priors[iclass],static_cast<float>(1.0-nbag)/nbag)*result[iclass];//add probabilities for each bag
+ // probOut[iclass][icol]*=pow(priors[iclass],static_cast<float>(1.0-nbag)/nbag)*prValues[iclass];//add probabilities for each bag
+ probOut[iclass][icol]*=pow(priors[iclass],static_cast<float>(1.0-nbag)/nbag)*result[iclass];//add probabilities for each bag
break;
case(2)://max rule
- // if(prValues[iclass]>prOut[iclass][icol])
- // prOut[iclass][icol]=prValues[iclass];
- if(result[iclass]>prOut[iclass])
- prOut[iclass]=result[iclass];
+ // if(prValues[iclass]>probOut[iclass][icol])
+ // probOut[iclass][icol]=prValues[iclass];
+ if(priors[iclass]*result[iclass]>probOut[iclass][icol])
+ probOut[iclass][icol]=priors[iclass]*result[iclass];
+ // if(result[iclass]>probOut[iclass][icol])
+ // probOut[iclass][icol]=result[iclass];
break;
}
if(classBag_opt.size()){
@@ -912,28 +922,28 @@ int main(int argc, char *argv[])
float maxBag2=0;//second max probability
float normBag=0;
for(short iclass=0;iclass<nclass;++iclass){
- if(prOut[iclass]>maxBag1){
- maxBag1=prOut[iclass];
+ if(probOut[iclass][icol]>maxBag1){
+ maxBag1=probOut[iclass][icol];
// classOut[icol]=classValueMap[type2string<short>(iclass)];
classOut[icol]=classValueMap[nameVector[iclass]];
// classOut[icol]=classValueMap[cm.getClass(iclass)];
}
- else if(prOut[iclass]>maxBag2)
- maxBag2=prOut[iclass];
- normBag+=prOut[iclass];
+ else if(probOut[iclass][icol]>maxBag2)
+ maxBag2=probOut[iclass][icol];
+ normBag+=probOut[iclass][icol];
}
- //normalize prOut and convert to percentage
+ //normalize probOut and convert to percentage
entropy[icol]=0;
for(short iclass=0;iclass<nclass;++iclass){
- float prv=prOut[iclass];
+ float prv=probOut[iclass][icol];
prv/=normBag;
entropy[icol]-=prv*log(prv)/log(2);
prv*=100.0;
- prOut[iclass]=static_cast<short>(prv+0.5);
- assert(classValueMap[nameVector[iclass]]<probOut.size());
- assert(classValueMap[nameVector[iclass]]>=0);
- probOut[classValueMap[nameVector[iclass]]][icol]=static_cast<short>(prv+0.5);
+ probOut[iclass][icol]=static_cast<short>(prv+0.5);
+ // assert(classValueMap[nameVector[iclass]]<probOut.size());
+ // assert(classValueMap[nameVector[iclass]]>=0);
+ // probOut[classValueMap[nameVector[iclass]]][icol]=static_cast<short>(prv+0.5);
}
entropy[icol]/=log(nclass)/log(2);
entropy[icol]=static_cast<short>(100*entropy[icol]+0.5);
@@ -1059,9 +1069,9 @@ int main(int argc, char *argv[])
imgReaderOgr.readData(validationPixel,OFTReal,fields,poFeature);
assert(validationPixel.size()==nband);
- vector<float> prOut(nclass);//posterior prob for each reclass
+ vector<float> probOut(nclass);//posterior prob for each reclass
for(short iclass=0;iclass<nclass;++iclass)
- prOut[iclass]=0;
+ probOut[iclass]=0;
for(int ibag=0;ibag<nbag;++ibag){
for(int iband=0;iband<nband;++iband){
validationFeature.push_back((validationPixel[iband]-offset[ibag][iband])/scale[ibag][iband]);
@@ -1111,14 +1121,14 @@ int main(int argc, char *argv[])
switch(comb_opt[0]){
default:
case(0)://sum rule
- prOut[iclass]+=result[iclass]+static_cast<float>(1.0-nbag)/nbag*priors[iclass];//add probabilities for each bag
+ probOut[iclass]+=result[iclass]+static_cast<float>(1.0-nbag)/nbag*priors[iclass];//add probabilities for each bag
break;
case(1)://product rule
- prOut[iclass]*=pow(priors[iclass],static_cast<float>(1.0-nbag)/nbag)*result[iclass];//add probabilities for each bag
+ probOut[iclass]*=pow(priors[iclass],static_cast<float>(1.0-nbag)/nbag)*result[iclass];//add probabilities for each bag
break;
case(2)://max rule
- if(result[iclass]>prOut[iclass])
- prOut[iclass]=result[iclass];
+ if(result[iclass]>probOut[iclass])
+ probOut[iclass]=result[iclass];
break;
}
}
@@ -1131,15 +1141,15 @@ int main(int argc, char *argv[])
string classOut="Unclassified";
for(short iclass=0;iclass<nclass;++iclass){
if(verbose_opt[0]>1)
- std::cout << prOut[iclass] << " ";
- if(prOut[iclass]>maxBag){
- maxBag=prOut[iclass];
+ std::cout << probOut[iclass] << " ";
+ if(probOut[iclass]>maxBag){
+ maxBag=probOut[iclass];
classOut=nameVector[iclass];
// classOut=cm.getClass(iclass);
//test
// assert(iclass==cm.getClassIndex(cm.getClass(iclass)));
}
- // normBag+=prOut[iclass];
+ // normBag+=probOut[iclass];
}
//look for class name
if(verbose_opt[0]>1){
@@ -1148,12 +1158,12 @@ int main(int argc, char *argv[])
else
std::cout << "->" << classOut << std::endl;
}
- //normalize prOut and convert to percentage
+ //normalize probOut and convert to percentage
// for(int iclass=0;iclass<nreclass;++iclass){
- // float prv=prOut[iclass];
+ // float prv=probOut[iclass];
// prv/=normBag;
// prv*=100.0;
- // prOut[iclass]=static_cast<short>(prv+0.5);
+ // probOut[iclass]=static_cast<short>(prv+0.5);
// }
if(output_opt.size()){
if(classValueMap.size())
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 12a8f3a..213fbe1 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -28,6 +28,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "base/Vector2d.h"
#include "algorithms/Filter2d.h"
#include "algorithms/Filter.h"
+#include "fileclasses/FileReaderAscii.h"
#include "imageclasses/ImgReaderGdal.h"
#include "imageclasses/ImgWriterGdal.h"
@@ -59,7 +60,7 @@ int main(int argc,char **argv) {
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<double> beta_opt("beta", "beta", "beta for Markov Random Field", 1.0);
+ Optionpk<string> beta_opt("beta", "beta", "ASCII file with beta for each class transition in Markov Random Field");
Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
bool doProcess;//stop process when program was invoked with help option (-h --help)
@@ -430,7 +431,31 @@ int main(int argc,char **argv) {
case(filter2d::mrf):{//Markov Random Field
if(verbose_opt[0])
std::cout << "Markov Random Field filtering" << std::endl;
- filter2d.mrf(input, output, dimX_opt[0], dimY_opt[0], beta_opt[0], true, down_opt[0], verbose_opt[0]);
+ if(beta_opt.size()){
+ //in file: classFrom classTo
+ //in variable: beta[classTo][classFrom]
+ FileReaderAscii betaReader(beta_opt[0]);
+ Vector2d<double> beta(class_opt.size(),class_opt.size());
+ vector<int> cols(class_opt.size());
+ for(int iclass=0;iclass<class_opt.size();++iclass)
+ cols[iclass]=iclass;
+ betaReader.readData(beta,cols);
+ if(verbose_opt[0]){
+ std::cout << "using values for beta:" << std::endl;
+ for(int iclass1=0;iclass1<class_opt.size();++iclass1)
+ std::cout << " " << iclass1 << " (" << class_opt[iclass1] << ")";
+ std::cout << std::endl;
+ for(int iclass1=0;iclass1<class_opt.size();++iclass1){
+ std::cout << iclass1 << " (" << class_opt[iclass1] << ")";
+ for(int iclass2=0;iclass2<class_opt.size();++iclass2)
+ std::cout << " " << beta[iclass2][iclass1] << " (" << class_opt[iclass2] << ")";
+ std::cout << std::endl;
+ }
+ }
+ filter2d.mrf(input, output, dimX_opt[0], dimY_opt[0], beta, true, down_opt[0], verbose_opt[0]);
+ }
+ else
+ filter2d.mrf(input, output, dimX_opt[0], dimY_opt[0], 1, true, down_opt[0], verbose_opt[0]);
break;
}
case(filter2d::sobelx):{//Sobel edge detection in X
--
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