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 6321cde98402f18cf6c50e0c2aa30e3bf706e41b
Author: Pieter Kempeneers <kempe...@gmail.com>
Date: Fri Oct 18 18:14:18 2013 +0200
support for even kernel size in Filter2d
---
src/algorithms/Filter2d.cc | 329 ++++++++++++++++++++-------------------------
src/algorithms/Filter2d.h | 154 ++++++++++++---------
src/apps/pkdiff.cc | 119 +++++++++++-----
src/apps/pkextract.cc | 52 ++++---
4 files changed, 351 insertions(+), 303 deletions(-)
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index 8c6e9ef..17df45f 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -84,23 +84,21 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input,
ImgWriterGdal& outpu
double progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
for(int iband=0;iband<input.nrOfBand();++iband){
- Vector2d<double> inBuffer(dimY);
+ Vector2d<double> inBuffer(dimY,input.nrOfCol());
vector<double> outBuffer(input.nrOfCol());
int indexI=0;
int indexJ=0;
//initialize last half of inBuffer
- for(int y=0;y<dimY;++y){
- inBuffer[y].resize(input.nrOfCol());
- if(y<dimY/2)
- continue;//skip first half
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
try{
- input.readData(inBuffer[y],GDT_Float64,indexJ,iband);
- ++indexJ;
+ input.readData(inBuffer[indexJ],GDT_Float64,abs(j),iband);
}
catch(string errorstring){
- cerr << errorstring << "in band " << iband << ", line " << indexJ <<
endl;
+ cerr << errorstring << "in line " << indexJ << endl;
}
+ ++indexJ;
}
+
for(int y=0;y<input.nrOfRow();++y){
if(y){//inBuffer already initialized for y=0
//erase first line from inBuffer
@@ -116,6 +114,13 @@ void filter2d::Filter2d::filter(const ImgReaderGdal&
input, ImgWriterGdal& outpu
cerr << errorstring << "in band " << iband << ", line " << y <<
endl;
}
}
+ else{
+ int over=y+dimY/2-input.nrOfRow();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
}
for(int x=0;x<input.nrOfCol();++x){
outBuffer[x]=0;
@@ -123,30 +128,30 @@ void filter2d::Filter2d::filter(const ImgReaderGdal&
input, ImgWriterGdal& outpu
bool masked=false;
if(noData){//only filter noData values
for(int imask=0;imask<m_mask.size();++imask){
- if(inBuffer[dimY/2][x]==m_mask[imask]){
+ if(inBuffer[(dimY-1)/2][x]==m_mask[imask]){
masked=true;
break;
}
}
if(!masked){
- outBuffer[x]=inBuffer[dimY/2][x];
+ outBuffer[x]=inBuffer[(dimY-1)/2][x];
continue;
}
}
assert(!noData||masked);
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
- for(int i=-dimX/2;i<(dimX+1)/2;++i){
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
indexI=x+i;
- indexJ=dimY/2+j;
+ indexJ=(dimY-1)/2+j;
//check if out of bounds
- if(x<dimX/2)
+ if(x<(dimX-1)/2)
indexI=x+abs(i);
- else if(x>=input.nrOfCol()-dimX/2)
+ else if(x>=input.nrOfCol()-(dimX-1)/2)
indexI=x-abs(i);
- if(y<dimY/2)
- indexJ=dimY/2+abs(j);
- else if(y>=input.nrOfRow()-dimY/2)
- indexJ=dimY/2-abs(j);
+ if(y<(dimY-1)/2)
+ indexJ=(dimY-1)/2+abs(j);
+ else if(y>=input.nrOfRow()-(dimY-1)/2)
+ indexJ=(dimY-1)/2-abs(j);
//do not take masked values into account
masked=false;
for(int imask=0;imask<m_mask.size();++imask){
@@ -156,8 +161,8 @@ void filter2d::Filter2d::filter(const ImgReaderGdal& input,
ImgWriterGdal& outpu
}
}
if(!masked){
-
outBuffer[x]+=(m_taps[dimY/2+j][dimX/2+i]*inBuffer[indexJ][indexI]);
- norm+=m_taps[dimY/2+j][dimX/2+i];
+
outBuffer[x]+=(m_taps[(dimY-1)/2+j][(dimX-1)/2+i]*inBuffer[indexJ][indexI]);
+ norm+=m_taps[(dimY-1)/2+j][(dimX-1)/2+i];
}
}
}
@@ -184,6 +189,12 @@ 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)
{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
bool usePriors=true;
if(prior.empty()){
cout << "no prior information" << endl;
@@ -195,6 +206,7 @@ void filter2d::Filter2d::majorVoting(const string&
inputFilename, const string&
cout << " " << static_cast<short>(prior[iclass]);
cout << endl;
}
+
ImgReaderGdal input;
ImgWriterGdal output;
input.open(inputFilename);
@@ -209,28 +221,25 @@ void filter2d::Filter2d::majorVoting(const string&
inputFilename, const string&
dimX=m_taps[0].size();
dimY=m_taps.size();
}
- Vector2d<double> inBuffer(dimY);
+
+ assert(dimX);
+ assert(dimY);
+
+ Vector2d<double> inBuffer(dimY,input.nrOfCol());
vector<double> outBuffer(input.nrOfCol());
- //initialize last half of inBuffer
int indexI=0;
int indexJ=0;
-// byte* tmpbuf=new byte[input.rowSize()];
- for(int y=0;y<dimY;++y){
- inBuffer[y].resize(input.nrOfCol());
- if(y<dimY/2)
- continue;//skip first half
- try{
- input.readData(inBuffer[y],GDT_Float64,indexJ++);
- }
- catch(string errorstring){
- cerr << errorstring << "in line " << indexJ << endl;
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ try{
+ input.readData(inBuffer[indexJ],GDT_Float64,abs(j));
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in line " << indexJ << endl;
+ }
+ ++indexJ;
}
- }
- const char* pszMessage;
- void* pProgressArg=NULL;
- GDALProgressFunc pfnProgress=GDALTermProgress;
- double progress=0;
- pfnProgress(progress,pszMessage,pProgressArg);
+
for(int y=0;y<input.nrOfRow();++y){
if(y){//inBuffer already initialized for y=0
//erase first line from inBuffer
@@ -246,23 +255,41 @@ void filter2d::Filter2d::majorVoting(const string&
inputFilename, const string&
cerr << errorstring << "in line" << y << endl;
}
}
+ else{
+ int over=y+dimY/2-input.nrOfRow();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
}
for(int x=0;x<input.nrOfCol();++x){
outBuffer[x]=0;
map<int,int> occurrence;
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
- for(int i=-dimX/2;i<(dimX+1)/2;++i){
+ 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){
indexI=x+i;
- indexJ=dimY/2+j;
//check if out of bounds
- if(x<dimX/2)
- indexI=x+abs(i);
- else if(x>=input.nrOfCol()-dimX/2)
- indexI=x-abs(i);
- if(y<dimY/2)
- indexJ=dimY/2+abs(j);
- else if(y>=input.nrOfRow()-dimY/2)
- indexJ=dimY/2-abs(j);
+ if(indexI<0)
+ indexI=-indexI;
+ else if(indexI>=input.nrOfCol())
+ indexI=input.nrOfCol()-i;
+ if(y+j<0)
+ indexJ=-j;
+ else if(y+j>=input.nrOfRow())
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ else
+ indexJ=(dimY-1)/2+j;
+
+ // if(x<dimX/2)
+ // indexI=x+abs(i);
+ // else if(x>=input.nrOfCol()-dimX/2)
+ // indexI=x-abs(i);
+ // if(y<dimY/2)
+ // indexJ=dimY/2+abs(j);
+ // else if(y>=input.nrOfRow()-dimY/2)
+ // indexJ=dimY/2-abs(j);
if(usePriors){
occurrence[inBuffer[indexJ][indexI]]+=prior[inBuffer[indexJ][indexI]-1];
}
@@ -275,10 +302,10 @@ void filter2d::Filter2d::majorVoting(const string&
inputFilename, const string&
if(mit->second>maxit->second)
maxit=mit;
}
- if(occurrence[inBuffer[dimY/2][x]]<maxit->second)//
+ if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
outBuffer[x]=maxit->first;
else//favorize original value in case of ties
- outBuffer[x]=inBuffer[dimY/2][x];
+ outBuffer[x]=inBuffer[(dimY-1)/2][x];
}
//write outBuffer to file
try{
@@ -294,95 +321,6 @@ void filter2d::Filter2d::majorVoting(const string&
inputFilename, const string&
output.close();
}
-// void Filter2d::homogeneousSpatial(const string& inputFilename, const
string& outputFilename, int dim, bool disc, int noValue)
-// {
-// ImgReaderGdal input;
-// ImgWriterGdal output;
-// input.open(inputFilename);
-// output.open(outputFilename,input);
-// int dimX=0;//horizontal!!!
-// int dimY=0;//vertical!!!
-// assert(dim);
-// dimX=dim;
-// dimY=dim;
-// Vector2d<double> inBuffer(dimY);
-// vector<double> outBuffer(input.nrOfCol());
-// //initialize last half of inBuffer
-// int indexI=0;
-// int indexJ=0;
-// for(int y=0;y<dimY;++y){
-// inBuffer[y].resize(input.nrOfCol());
-// if(y<dimY/2)
-// continue;//skip first half
-// try{
-// input.readData(inBuffer[y],GDT_Float64,indexJ++);
-// }
-// catch(string errorstring){
-// cerr << errorstring << "in line " << indexJ << endl;
-// }
-// }
-// const char* pszMessage;
-// void* pProgressArg=NULL;
-// GDALProgressFunc pfnProgress=GDALTermProgress;
-// double progress=0;
-// pfnProgress(progress,pszMessage,pProgressArg);
-// for(int y=0;y<input.nrOfRow();++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<input.nrOfRow()){
-// //allocate buffer
-// inBuffer.push_back(inBuffer.back());
-// try{
-// input.readData(inBuffer[inBuffer.size()-1],GDT_Float64,y+dimY/2);
-// }
-// catch(string errorstring){
-// cerr << errorstring << "in line " << y << endl;
-// }
-// }
-// }
-// for(int x=0;x<input.nrOfCol();++x){
-// outBuffer[x]=0;
-// map<int,int> occurrence;
-// for(int j=-dimY/2;j<(dimY+1)/2;++j){
-// for(int i=-dimX/2;i<(dimX+1)/2;++i){
-// if(disc&&(i*i+j*j>(dim/2)*(dim/2)))
-// continue;
-// indexI=x+i;
-// //check if out of bounds
-// if(indexI<0)
-// indexI=-indexI;
-// else if(indexI>=input.nrOfCol())
-// indexI=input.nrOfCol()-indexI;
-// if(y+j<0)
-// indexJ=-j;
-// else if(y+j>=input.nrOfRow())
-// indexJ=dimY/2-j;
-// else
-// indexJ=dimY/2+j;
-// ++occurrence[inBuffer[indexJ][indexI]];
-// }
-// }
-// if(occurrence.size()==1)//all values in window must be the same
-// outBuffer[x]=inBuffer[dimY/2][x];
-// else//favorize original value in case of ties
-// outBuffer[x]=noValue;
-// }
-// //write outBuffer to file
-// try{
-// output.writeData(outBuffer,GDT_Float64,y);
-// }
-// catch(string errorstring){
-// cerr << errorstring << "in line" << y << endl;
-// }
-// progress=(1.0+y)/input.nrOfRow();
-// pfnProgress(progress,pszMessage,pProgressArg);
-// }
-// input.close();
-// output.close();
-// }
-
void filter2d::Filter2d::median(const string& inputFilename, const string&
outputFilename,int dim, bool disc)
{
ImgReaderGdal input;
@@ -408,21 +346,23 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal&
output, const std::string& method, int dimX, int dimY, short down, bool disc)
{
- assert(dimX);
- assert(dimY);
- statfactory::StatFactory stat;
const char* pszMessage;
void* pProgressArg=NULL;
GDALProgressFunc pfnProgress=GDALTermProgress;
double progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
+
+ assert(dimX);
+ assert(dimY);
+
+ 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);
- //initialize last half of inBuffer
int indexI=0;
int indexJ=0;
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
try{
input.readData(inBuffer[indexJ],GDT_Float64,abs(j),iband);
}
@@ -462,8 +402,9 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
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){
+ 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){
double d2=i*i+j*j;//square distance
if(disc&&(d2>(dimX/2)*(dimY/2)))
continue;
@@ -476,9 +417,9 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
if(y+j<0)
indexJ=-j;
else if(y+j>=input.nrOfRow())
- indexJ=dimY/2-j;
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
else
- indexJ=dimY/2+j;
+ indexJ=(dimY-1)/2+j;
bool masked=false;
for(int imask=0;imask<m_mask.size();++imask){
if(inBuffer[indexJ][indexI]==m_mask[imask]){
@@ -539,7 +480,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
-
outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
+ outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[centre])?
1:0;
break;
}
case(filter2d::minmax):{
@@ -552,7 +493,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
if(min!=max)
outBuffer[x/down]=0;
else
- outBuffer[x/down]=windowBuffer[dimX*dimY/2];//centre pixels
+ outBuffer[x/down]=windowBuffer[centre];//centre pixels
}
break;
}
@@ -567,7 +508,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
-
outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
+ outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[centre])?
1:0;
break;
}
case(filter2d::order):{
@@ -579,7 +520,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
double theMin=stat.min(windowBuffer);
double theMax=stat.max(windowBuffer);
double scale=(ubound-lbound)/(theMax-theMin);
-
outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[dimX*dimY/2]-theMin)+lbound);
+
outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[centre]-theMin)+lbound);
}
break;
}
@@ -589,7 +530,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
}
case(filter2d::homog):
if(occurrence.size()==1)//all values in window must be the same
- outBuffer[x/down]=inBuffer[dimY/2][x];
+ outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
else//favorize original value in case of ties
outBuffer[x/down]=m_noValue;
break;
@@ -597,9 +538,9 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
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])
+ else if(*wit!=inBuffer[(dimY-1)/2][x])
outBuffer[x/down]=1;
- else if(*wit==inBuffer[dimY/2][x]){//todo:wit mag niet central
pixel zijn
+ else if(*wit==inBuffer[(dimY-1)/2][x]){//todo:wit mag niet central
pixel zijn
outBuffer[x/down]=m_noValue;
break;
}
@@ -623,10 +564,10 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
if(mit->second>maxit->second)
maxit=mit;
}
- if(occurrence[inBuffer[dimY/2][x]]<maxit->second)//
+ if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
outBuffer[x/down]=maxit->first;
else//favorize original value in case of ties
- outBuffer[x/down]=inBuffer[dimY/2][x];
+ outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
}
else
outBuffer[x/down]=m_noValue;
@@ -635,7 +576,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input,
ImgWriterGdal& output,
case(filter2d::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)
+ outBuffer[x/down]=inBuffer[(dimY-1)/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];
@@ -712,23 +653,25 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input,
ImgWriterGdal& output,
//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);
const char* pszMessage;
void* pProgressArg=NULL;
GDALProgressFunc pfnProgress=GDALTermProgress;
double progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
+
+ assert(dimX);
+ assert(dimY);
+
Vector2d<short> inBuffer(dimY,input.nrOfCol());
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;
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
try{
input.readData(inBuffer[indexJ],GDT_Int16,abs(j));
}
@@ -771,8 +714,9 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input,
ImgWriterGdal& output,
outBuffer[iclass][x/down]=0;
}
vector<double> windowBuffer;
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
- for(int i=-dimX/2;i<(dimX+1)/2;++i){
+ 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){
if(i!=0&&j!=0&&!eightConnectivity)
continue;
if(i==0&&j==0)
@@ -786,9 +730,9 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input,
ImgWriterGdal& output,
if(y+j<0)
indexJ=-j;
else if(y+j>=input.nrOfRow())
- indexJ=dimY/2-j;
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
else
- indexJ=dimY/2+j;
+ indexJ=(dimY-1)/2+j;
bool masked=false;
for(int imask=0;imask<m_mask.size();++imask){
if(inBuffer[indexJ][indexI]==m_mask[imask]){
@@ -797,7 +741,6 @@ void filter2d::Filter2d::mrf(const ImgReaderGdal& input,
ImgWriterGdal& output,
}
}
if(!masked){
- // if(inBuffer[dimY/2][x]//centre pixel
for(int iclass=0;iclass<m_class.size();++iclass){
if(inBuffer[indexJ][indexI]==m_class[iclass])
potential[iclass]+=1;
@@ -967,16 +910,23 @@ void filter2d::Filter2d::shift(const ImgReaderGdal&
input, ImgWriterGdal& output
void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal&
output, const std::string& method, int dimX, int dimY, const vector<double>
&angle, bool disc)
{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
assert(dimX);
assert(dimY);
+
statfactory::StatFactory stat;
for(int iband=0;iband<input.nrOfBand();++iband){
Vector2d<double> inBuffer(dimY,input.nrOfCol());
vector<double> outBuffer(input.nrOfCol());
- //initialize last half of inBuffer
int indexI=0;
int indexJ=0;
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
try{
input.readData(inBuffer[indexJ],GDT_Float64,abs(j),iband);
++indexJ;
@@ -985,11 +935,6 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal&
input, ImgWriterGdal& o
cerr << errorstring << "in line " << indexJ << endl;
}
}
- const char* pszMessage;
- void* pProgressArg=NULL;
- GDALProgressFunc pfnProgress=GDALTermProgress;
- double progress=0;
- pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
if(y){//inBuffer already initialized for y=0
//erase first line from inBuffer
@@ -1005,12 +950,20 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal&
input, ImgWriterGdal& o
cerr << errorstring << "in band " << iband << ", line " << y <<
endl;
}
}
+ else{
+ int over=y+dimY/2-input.nrOfRow();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
}
for(int x=0;x<input.nrOfCol();++x){
- double currentValue=inBuffer[dimY/2][x];
+ double currentValue=inBuffer[(dimY-1)/2][x];
outBuffer[x]=currentValue;
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]){
currentMasked=true;
@@ -1021,11 +974,11 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal&
input, ImgWriterGdal& o
outBuffer[x]=currentValue;
}
else{
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
- for(int i=-dimX/2;i<(dimX+1)/2;++i){
- if(disc&&(i*i+j*j>(dimX/2)*(dimY/2)))
- continue;
- bool masked=false;
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ double d2=i*i+j*j;//square distance
+ if(disc&&(d2>(dimX/2)*(dimY/2)))
+ continue;
if(angle.size()){
double theta;
//use polar coordinates in radians
@@ -1069,10 +1022,14 @@ void filter2d::Filter2d::morphology(const
ImgReaderGdal& input, ImgWriterGdal& o
indexI=input.nrOfCol()-i;
if(y+j<0)
indexJ=-j;
- else if(y+j>=input.nrOfRow())
- indexJ=dimY/2-j;//indexJ=inBuffer.size()-1-j;
- else
- indexJ=dimY/2+j;
+ else if(y+j>=input.nrOfRow())
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ else
+ indexJ=(dimY-1)/2+j;
+ //todo: introduce novalue as this: ?
+ // if(inBuffer[indexJ][indexI]==m_noValue)
+ // continue;
+ bool masked=false;
for(int imask=0;imask<m_mask.size();++imask){
if(inBuffer[indexJ][indexI]==m_mask[imask]){
masked=true;
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 4b7a8a0..ea5d690 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -182,35 +182,45 @@ private:
//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++];
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ inBuffer[indexJ]=inputVector[abs(j)];
+ ++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())
+ if(y+dimY/2<inputVector.size()){
+ //allocate buffer
inBuffer.push_back(inputVector[y+dimY/2]);
+ }
+ else{
+ int over=y+dimY/2-inputVector.nRows();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
}
- for(int x=0;x<inputVector[0].size();++x){
+ for(int x=0;x<inputVector.nCols();++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]);
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ indexI=x+i;
+ indexJ=(dimY-1)/2+j;
+ //check if out of bounds
+ if(x<(dimX-1)/2)
+ indexI=x+abs(i);
+ else if(x>=inputVector.nCols()-(dimX-1)/2)
+ indexI=x-abs(i);
+ if(y<(dimY-1)/2)
+ indexJ=(dimY-1)/2+abs(j);
+ else if(y>=inputVector.nRows()-(dimY-1)/2)
+ indexJ=(dimY-1)/2-abs(j);
+
outBuffer[x]+=(m_taps[(dimY-1)/2+j][(dimX-1)/2+i]*inBuffer[indexJ][indexI]);
}
}
}
@@ -221,23 +231,26 @@ private:
template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>&
inputVector, Vector2d<T2>& outputVector, const std::string& method, int dimX,
int dimY, short down, bool disc)
{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ assert(dimX);
+ assert(dimY);
+
statfactory::StatFactory stat;
outputVector.resize((inputVector.size()+down-1)/down);
Vector2d<T1> inBuffer(dimY);
std::vector<T2> outBuffer((inputVector[0].size()+down-1)/down);
- //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++];
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ inBuffer[indexJ]=inputVector[abs(j)];
+ ++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
@@ -261,19 +274,21 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
outBuffer[x/down]=0;
std::vector<double> windowBuffer;
std::map<int,int> occurrence;
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
- for(int i=-dimX/2;i<(dimX+1)/2;++i){
+ 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){
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);
+ if(indexI<0)
+ indexI=-indexI;
+ else if(indexI>=inputVector[0].size())
+ indexI=inputVector[0].size()-i;
+ if(y+j<0)
+ indexJ=-j;
+ else if(y+j>=inputVector.size())
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ 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]){
@@ -335,7 +350,7 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
-
outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
+ outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[centre])?
1:0;
break;
}
case(filter2d::minmax):{
@@ -348,7 +363,7 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
if(min!=max)
outBuffer[x/down]=0;
else
- outBuffer[x/down]=windowBuffer[dimX*dimY/2];//centre pixels
+ outBuffer[x/down]=windowBuffer[centre];//centre pixels
}
break;
}
@@ -363,7 +378,7 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
-
outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
+ outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[centre])?
1:0;
break;
}
case(filter2d::order):{
@@ -375,7 +390,7 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
double theMin=stat.min(windowBuffer);
double theMax=stat.max(windowBuffer);
double scale=(ubound-lbound)/(theMax-theMin);
-
outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[dimX*dimY/2]-theMin)+lbound);
+
outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[centre]-theMin)+lbound);
}
break;
}
@@ -385,7 +400,7 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
}
case(filter2d::homog):
if(occurrence.size()==1)//all values in window must be the same
- outBuffer[x/down]=inBuffer[dimY/2][x];
+ outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
else//favorize original value in case of ties
outBuffer[x/down]=m_noValue;
break;
@@ -393,9 +408,9 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
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/2][x])
+ else if(*wit!=inBuffer[(dimY-1)/2][x])
outBuffer[x/down]=1;
- else if(*wit==inBuffer[dimY/2][x]){//todo:wit mag niet central pixel
zijn
+ else if(*wit==inBuffer[(dimY-1)/2][x]){//todo:wit mag niet central
pixel zijn
outBuffer[x/down]=m_noValue;
break;
}
@@ -419,10 +434,10 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
if(mit->second>maxit->second)
maxit=mit;
}
- if(occurrence[inBuffer[dimY/2][x]]<maxit->second)//
+ if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
outBuffer[x/down]=maxit->first;
else//favorize original value in case of ties
- outBuffer[x/down]=inBuffer[dimY/2][x];
+ outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
}
else
outBuffer[x/down]=m_noValue;
@@ -431,7 +446,7 @@ template<class T1, class T2> void Filter2d::doit(const
Vector2d<T1>& inputVector
case(filter2d::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)
+ outBuffer[x/down]=inBuffer[(dimY-1)/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];
@@ -573,6 +588,12 @@ template<class T> void Filter2d::shift(const Vector2d<T>&
input, Vector2d<T>& ou
template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>&
input, Vector2d<T>& output, const std::string& method, int dimX, int dimY, bool
disc, double hThreshold)
{
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
unsigned long int nchange=0;
assert(dimX);
assert(dimY);
@@ -580,19 +601,14 @@ template<class T> unsigned long int
Filter2d::morphology(const Vector2d<T>& inpu
Vector2d<T> inBuffer(dimY,input.nCols());
output.clear();
output.resize(input.nRows(),input.nCols());
- //initialize last half of inBuffer
int indexI=0;
int indexJ=0;
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
+ //initialize last half of inBuffer
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
for(int i=0;i<input.nCols();++i)
inBuffer[indexJ][i]=input[abs(j)][i];
++indexJ;
}
- const char* pszMessage;
- void* pProgressArg=NULL;
- GDALProgressFunc pfnProgress=GDALTermProgress;
- double progress=0;
- pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nRows();++y){
if(y){//inBuffer already initialized for y=0
//erase first line from inBuffer
@@ -604,10 +620,17 @@ template<class T> unsigned long int
Filter2d::morphology(const Vector2d<T>& inpu
for(int i=0;i<input.nCols();++i)
inBuffer[inBuffer.size()-1][i]=input[y+dimY/2][i];
}
+ else{
+ int over=y+dimY/2-input.nRows();
+ int index=(inBuffer.size()-1)-over;
+ assert(index>=0);
+ assert(index<inBuffer.size());
+ inBuffer.push_back(inBuffer[index]);
+ }
}
for(int x=0;x<input.nCols();++x){
output[y][x]=0;
- double currentValue=inBuffer[dimY/2][x];
+ double currentValue=inBuffer[(dimY-1)/2][x];
std::vector<double> statBuffer;
bool currentMasked=false;
for(int imask=0;imask<m_mask.size();++imask){
@@ -621,9 +644,10 @@ template<class T> unsigned long int
Filter2d::morphology(const Vector2d<T>& inpu
output[y][x]=currentValue;
}
else{
- for(int j=-dimY/2;j<(dimY+1)/2;++j){
- for(int i=-dimX/2;i<(dimX+1)/2;++i){
- if(disc&&(i*i+j*j>(dimX/2)*(dimY/2)))
+ for(int j=-(dimY-1)/2;j<=dimY/2;++j){
+ for(int i=-(dimX-1)/2;i<=dimX/2;++i){
+ double d2=i*i+j*j;//square distance
+ if(disc&&(d2>(dimX/2)*(dimY/2)))
continue;
indexI=x+i;
//check if out of bounds
@@ -634,9 +658,9 @@ template<class T> unsigned long int
Filter2d::morphology(const Vector2d<T>& inpu
if(y+j<0)
indexJ=-j;
else if(y+j>=input.nRows())
- indexJ=dimY/2-j;//indexJ=inBuffer.size()-1-j;
- else
- indexJ=dimY/2+j;
+ indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
+ else
+ indexJ=(dimY-1)/2+j;
if(inBuffer[indexJ][indexI]==m_noValue)
continue;
bool masked=false;
diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc
index 7b1b1f8..7a123d8 100644
--- a/src/apps/pkdiff.cc
+++ b/src/apps/pkdiff.cc
@@ -42,11 +42,13 @@ int main(int argc, char *argv[])
Optionpk<short> lzw_opt("\0", "lzw", "compression (default value is 1)", 1);
Optionpk<string> labelref_opt("lr", "lref", "name of the reference label in
case reference is shape file(default is label)", "label");
Optionpk<string> labelclass_opt("lc", "lclass", "name of the classified
label in case output is shape file (default is class)", "class");
- Optionpk<short> class_opt("c", "class", "numeric classes used (must cover
range in input and reference raster image. Leave empty if range must be read
from first input image (default)", 0);
+ // Optionpk<short> class_opt("c", "class", "numeric classes used (must cover
range in input and reference raster image. Leave empty if range must be read
from first input image (default)", 0);
Optionpk<short> boundary_opt("\0", "boundary", "boundary for selecting the
sample (default: 1)", 1);
Optionpk<bool> disc_opt("\0", "circular", "use circular disc kernel
boundary)", false);
Optionpk<bool> homogeneous_opt("\0", "homogeneous", "only take homogeneous
regions into account", false);
Optionpk<string> option_opt("co", "co", "options: NAME=VALUE [-co
COMPRESS=LZW] [-co INTERLEAVE=BAND]");
+ Optionpk<string> classname_opt("\0", "class", "list of class names.");
+ Optionpk<short> classvalue_opt("\0", "reclass", "list of class values (use
same order as in classname opt.");
Optionpk<short> verbose_opt("v", "verbose", "verbose (default value is 0)",
0);
bool doProcess;//stop process when program was invoked with help option (-h
--help)
@@ -67,10 +69,12 @@ int main(int argc, char *argv[])
lzw_opt.retrieveOption(argc,argv);
labelref_opt.retrieveOption(argc,argv);
labelclass_opt.retrieveOption(argc,argv);
- class_opt.retrieveOption(argc,argv);
+ // class_opt.retrieveOption(argc,argv);
boundary_opt.retrieveOption(argc,argv);
disc_opt.retrieveOption(argc,argv);
homogeneous_opt.retrieveOption(argc,argv);
+ classname_opt.retrieveOption(argc,argv);
+ classvalue_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
catch(string predefinedString){
@@ -97,11 +101,15 @@ int main(int argc, char *argv[])
vector<short> referenceRange;
ConfusionMatrix cm;
int nclass=0;
+ map<string,short> classValueMap;
+ vector<std::string> nameVector(255);//the inverse of the classValueMap
vector<string> classNames;
if(confusion_opt[0]){
- if(class_opt.size()>1)
- inputRange=class_opt;
- else{
+ // if(class_opt.size()>1)
+ // inputRange=class_opt;
+ // if(classvalue_opt.size()>1)
+ // inputRange=classvalue_opt;
+ // else{
try{
if(verbose_opt[0])
cout << "opening input image file " << input_opt[0] << endl;
@@ -113,7 +121,7 @@ int main(int argc, char *argv[])
}
inputReader.getRange(inputRange,band_opt[0]);
inputReader.close();
- }
+ // }
for(int iflag=0;iflag<flag_opt.size();++iflag){
vector<short>::iterator fit;
@@ -126,6 +134,15 @@ int main(int argc, char *argv[])
cout << "nclass (inputRange.size()): " << nclass << endl;
cout << "input range: " << endl;
}
+ if(classname_opt.size()){
+ assert(classname_opt.size()==classvalue_opt.size());
+ for(int iclass=0;iclass<classname_opt.size();++iclass){
+ classValueMap[classname_opt[iclass]]=classvalue_opt[iclass];
+ assert(classvalue_opt[iclass]<nameVector.size());
+ nameVector[classvalue_opt[iclass]]=classname_opt[iclass];
+ }
+ }
+ // nclass=classValueMap.size();
for(int rc=0;rc<inputRange.size();++rc){
classNames.push_back(type2string(inputRange[rc]));
if(verbose_opt[0])
@@ -138,6 +155,7 @@ int main(int argc, char *argv[])
cout << iclass << " " << cm.getClass(iclass) << endl;
}
}
+
unsigned int ntotalValidation=0;
unsigned int nflagged=0;
Vector2d<int> resultClass(nclass,nclass);
@@ -243,10 +261,24 @@ int main(int argc, char *argv[])
//get x and y from readFeature
double x,y;
OGRGeometry *poGeometry;
+ OGRPolygon readPolygon;
+ OGRPoint centroidPoint;
+ OGRPoint *poPoint;
poGeometry = readFeature->GetGeometryRef();
- assert( poGeometry != NULL &&
wkbFlatten(poGeometry->getGeometryType()) == wkbPoint );
- OGRPoint *poPoint = (OGRPoint *) poGeometry;
- // writeFeature->SetGeometry(poPoint);
+ // assert( poGeometry != NULL &&
wkbFlatten(poGeometry->getGeometryType()) == wkbPoint );
+ if(poGeometry==NULL)
+ continue;
+ else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
+ readPolygon = *((OGRPolygon *) poGeometry);
+ readPolygon.Centroid(¢roidPoint);
+ poPoint=¢roidPoint;
+ }
+ else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
+ poPoint = (OGRPoint *) poGeometry;
+ else{
+ std::cerr << "Warning: skipping feature (not of type point or
polygon)" << std::endl;
+ continue;
+ }
x=poPoint->getX();
y=poPoint->getY();
short inputValue;
@@ -255,7 +287,14 @@ int main(int argc, char *argv[])
short maskValue;
short outputValue;
//read referenceValue from feature
- short
referenceValue=readFeature->GetFieldAsInteger(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
+ short referenceValue;
+ string referenceClassName;
+ if(classValueMap.size()){
+
referenceClassName=readFeature->GetFieldAsString(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
+ referenceValue=classValueMap[referenceClassName];
+ }
+ else
+
referenceValue=readFeature->GetFieldAsInteger(readFeature->GetFieldIndex(labelref_opt[0].c_str()));
if(verbose_opt[0])
cout << "reference value: " << referenceValue << endl;
bool pixelFlagged=false;
@@ -358,15 +397,22 @@ int main(int argc, char *argv[])
writeFeature->SetField(labelclass_opt[0].c_str(),static_cast<int>(inputValue));
if(confusion_opt[0]){
++ntotalValidation;
- int
rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),referenceValue));
- int
ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),inputValue));
- assert(rc<nclass);
- assert(ic<nclass);
- ++nvalidation[rc];
- ++resultClass[rc][ic];
- if(verbose_opt[0]>1)
- cout << "increment: " << rc << " " << referenceRange[rc] <<
" " << ic << " " << inputRange[ic] << endl;
- cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
+ if(classValueMap.size()){
+ assert(inputValue<nameVector.size());
+ string className=nameVector[inputValue];
+
cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1);
+ }
+ else{
+ int
rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),referenceValue));
+ int
ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),inputValue));
+ assert(rc<nclass);
+ assert(ic<nclass);
+ ++nvalidation[rc];
+ ++resultClass[rc][ic];
+ if(verbose_opt[0]>1)
+ cout << "increment: " << rc << " " << referenceRange[rc]
<< " " << ic << " " << inputRange[ic] << endl;
+ cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
+ }
}
if(inputValue==referenceValue){//correct
if(valueE_opt[0]!=flag_opt[0])
@@ -404,19 +450,26 @@ int main(int argc, char *argv[])
if(!windowJ&&!windowI){//centre pixel
if(confusion_opt[0]){
++ntotalValidation;
- int
rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),referenceValue));
- int
ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),inputValue));
- if(rc>=nclass)
- continue;
- if(ic>=nclass)
- continue;
- // assert(rc<nclass);
- // assert(ic<nclass);
- ++nvalidation[rc];
- ++resultClass[rc][ic];
- if(verbose_opt[0]>1)
- cout << "increment: " << rc << " " <<
referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
- cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
+ if(classValueMap.size()){
+ assert(inputValue<nameVector.size());
+ string className=nameVector[inputValue];
+
cm.incrementResult(type2string<short>(classValueMap[referenceClassName]),type2string<short>(classValueMap[className]),1);
+ }
+ else{
+ int
rc=distance(referenceRange.begin(),find(referenceRange.begin(),referenceRange.end(),referenceValue));
+ int
ic=distance(inputRange.begin(),find(inputRange.begin(),inputRange.end(),inputValue));
+ if(rc>=nclass)
+ continue;
+ if(ic>=nclass)
+ continue;
+ // assert(rc<nclass);
+ // assert(ic<nclass);
+ ++nvalidation[rc];
+ ++resultClass[rc][ic];
+ if(verbose_opt[0]>1)
+ cout << "increment: " << rc << " " <<
referenceRange[rc] << " " << ic << " " << inputRange[ic] << endl;
+ cm.incrementResult(cm.getClass(rc),cm.getClass(ic),1);
+ }
}
if(inputValue==referenceValue){//correct
if(valueE_opt[0]!=flag_opt[0])
@@ -453,7 +506,7 @@ int main(int argc, char *argv[])
maskReader.close();
}
}
- }
+ }//reference is shape file
else{
ImgWriterGdal imgWriter;
try{
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index 90b3028..5129206 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -34,7 +34,7 @@ along with pktools. If not, see
<http://www.gnu.org/licenses/>.
#endif
namespace rule{
- enum RULE_TYPE {point=0, mean=1, proportion=2, custom=3, minimum=4,
maximum=5, maxvote=6};
+ enum RULE_TYPE {point=0, mean=1, proportion=2, custom=3, minimum=4,
maximum=5, maxvote=6, centroid=7};
}
int main(int argc, char *argv[])
@@ -62,7 +62,7 @@ int main(int argc, char *argv[])
Optionpk<string> label_opt("cn", "cname", "name of the class label in the
output vector file", "label");
Optionpk<bool> polygon_opt("l", "line", "create OGRPolygon as geometry
instead of OGRPoint. Only if sample option is also of polygon type.", false);
Optionpk<int> band_opt("b", "band", "band index to crop. Use -1 to use all
bands)", -1);
- Optionpk<string> rule_opt("r", "rule", "rule how to report image information
per feature. point (value at each point or at centroid of the polygon if line
is set), mean (mean value written to centroid of polygon if line is set),
proportion, minimum (of polygon), maximum (of polygon), maxvote.", "point");
+ Optionpk<string> rule_opt("r", "rule", "rule how to report image information
per feature. point (value at each point or at centroid of the polygon if line
is set), centroid, mean (mean value written to centroid of polygon if line is
set), proportion, minimum (of polygon), maximum (of polygon), maxvote.",
"point");
Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
bool doProcess;//stop process when program was invoked with help option (-h
--help)
@@ -106,6 +106,7 @@ int main(int argc, char *argv[])
//initialize ruleMap
ruleMap["point"]=rule::point;
ruleMap["mean"]=rule::mean;
+ ruleMap["centroid"]=rule::centroid;
ruleMap["proportion"]=rule::proportion;
ruleMap["custom"]=rule::custom;
ruleMap["maxvote"]=rule::maxvote;
@@ -797,6 +798,7 @@ int main(int argc, char *argv[])
break;
case(rule::point):
case(rule::mean):
+ case(rule::centroid):
default:{
if(keepFeatures_opt[0]){
ogrWriter.createField("origId",OFTInteger);//the fieldId of the
original feature
@@ -1147,7 +1149,7 @@ int main(int argc, char *argv[])
double ulx,uly,lrx,lry;
double uli,ulj,lri,lrj;
- if(polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point){
+
if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
ulx=writeCentroidPoint.getX();
uly=writeCentroidPoint.getY();
lrx=ulx;
@@ -1194,7 +1196,7 @@ int main(int argc, char *argv[])
else
writePolygonFeature =
OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- else if(ruleMap[rule_opt[0]]==rule::mean){
+ else
if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
if(writeTest)
writeCentroidFeature =
OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
@@ -1205,6 +1207,7 @@ int main(int argc, char *argv[])
switch(ruleMap[rule_opt[0]]){
case(rule::point):
case(rule::mean):
+ case(rule::centroid):
default:
polyValues.resize(nband);
break;
@@ -1318,7 +1321,7 @@ int main(int argc, char *argv[])
OGRFeature *writePointFeature;
if(!polygon_opt[0]){
//create feature
- if(ruleMap[rule_opt[0]]!=rule::mean){//do not create in
case of mean value (only create point at centroid
+
if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do
not create in case of mean value (only create point at centroid
if(writeTest)
writePointFeature =
OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
@@ -1343,16 +1346,17 @@ int main(int argc, char *argv[])
imgReader.readData(value,GDT_Float64,i,j,theBand);
if(verbose_opt[0]>1)
std::cout << ": " << value << std::endl;
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean){
+
if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
int iclass=0;
switch(ruleMap[rule_opt[0]]){
case(rule::point)://in centroid if polygon_opt==true or
all values as points if polygon_opt!=true
+ case(rule::centroid):
default:
polyValues[iband]=value;
break;
case(rule::mean)://mean as polygon if polygon_opt==true
or as point in centroid if polygon_opt!=true
polyValues[iband]+=value;
- break;
+ break;
case(rule::proportion):
case(rule::custom):
case(rule::minimum):
@@ -1410,7 +1414,7 @@ int main(int argc, char *argv[])
if(!polygon_opt[0]){
// if(keepFeatures_opt[0])
//
writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(ruleMap[rule_opt[0]]!=rule::mean){//do not create in
case of mean value (only at centroid)
+
if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do
not create in case of mean value (only at centroid)
if(keepFeatures_opt[0])
writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
//write feature
@@ -1439,7 +1443,7 @@ int main(int argc, char *argv[])
}
}
}
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean){
+
if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
//do not create if no points found within polygon
if(!nPointPolygon)
continue;
@@ -1471,7 +1475,7 @@ int main(int argc, char *argv[])
std::cout << "write feature has " <<
writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
}
switch(ruleMap[rule_opt[0]]){
- case(rule::point)://value at each point (or at centroid of
polygon if line is not set
+ case(rule::point)://value at each point (or at centroid of
polygon if line is set)
default:{
if(verbose_opt[0])
std::cout << "number of points in polygon: " <<
nPointPolygon << std::endl;
@@ -1543,9 +1547,13 @@ int main(int argc, char *argv[])
}
break;
}
- case(rule::mean):{//mean value (written to centroid of polygon
if line is not set
+ case(rule::mean):
+ case(rule::centroid):{//mean value (written to centroid of
polygon if line is not set)
if(verbose_opt[0])
std::cout << "number of points in polygon: " <<
nPointPolygon << std::endl;
+ //test
+ if(ruleMap[rule_opt[0]]==rule::centroid)
+ assert(nPointPolygon<=1);
for(int index=0;index<polyValues.size();++index){
double theValue=polyValues[index];
// ostringstream fs;
@@ -1756,7 +1764,7 @@ int main(int argc, char *argv[])
double ulx,uly,lrx,lry;
double uli,ulj,lri,lrj;
- if(polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point){
+
if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
ulx=writeCentroidPoint.getX();
uly=writeCentroidPoint.getY();
lrx=ulx;
@@ -1803,7 +1811,7 @@ int main(int argc, char *argv[])
else
writePolygonFeature =
OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- else if(ruleMap[rule_opt[0]]==rule::mean){
+ else
if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
if(writeTest)
writeCentroidFeature =
OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
@@ -1814,6 +1822,7 @@ int main(int argc, char *argv[])
switch(ruleMap[rule_opt[0]]){
case(rule::point):
case(rule::mean):
+ case(rule::centroid):
default:
polyValues.resize(nband);
break;
@@ -1927,7 +1936,7 @@ int main(int argc, char *argv[])
OGRFeature *writePointFeature;
if(!polygon_opt[0]){
//create feature
- if(ruleMap[rule_opt[0]]!=rule::mean){//do not create in
case of mean value (only create point at centroid
+
if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do
not create in case of mean value (only create point at centroid)
if(writeTest)
writePointFeature =
OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
@@ -1952,10 +1961,11 @@ int main(int argc, char *argv[])
imgReader.readData(value,GDT_Float64,i,j,theBand);
if(verbose_opt[0]>1)
std::cout << ": " << value << std::endl;
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean){
+
if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
int iclass=0;
switch(ruleMap[rule_opt[0]]){
case(rule::point)://in centroid if polygon_opt==true or
all values as points if polygon_opt!=true
+ case(rule::centroid):
default:
polyValues[iband]=value;
break;
@@ -2019,7 +2029,7 @@ int main(int argc, char *argv[])
if(!polygon_opt[0]){
if(keepFeatures_opt[0])
writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(ruleMap[rule_opt[0]]!=rule::mean){//do not create in
case of mean value (only at centroid)
+
if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do
not create in case of mean value (only at centroid)
//write feature
if(verbose_opt[0]>1)
std::cout << "creating point feature" << std::endl;
@@ -2046,7 +2056,7 @@ int main(int argc, char *argv[])
}
}
}
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean){
+
if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
//add ring to polygon
if(polygon_opt[0]){
writePolygon.addRing(&writeRing);
@@ -2075,7 +2085,7 @@ int main(int argc, char *argv[])
std::cout << "write feature has " <<
writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
}
switch(ruleMap[rule_opt[0]]){
- case(0)://value at each point (or at centroid of polygon if line
is not set
+ case(rule::point)://value at each point (or at centroid of
polygon if line is set)
default:{
if(verbose_opt[0])
std::cout << "number of points in polygon: " <<
nPointPolygon << std::endl;
@@ -2143,9 +2153,13 @@ int main(int argc, char *argv[])
}//for index
break;
}//case 0 and default
- case(rule::mean):{//mean value (written to centroid of polygon
if line is not set
+ case(rule::mean):
+ case(rule::centroid):{//mean value (written to centroid of
polygon if line is not set
if(verbose_opt[0])
std::cout << "number of points in polygon: " <<
nPointPolygon << std::endl;
+ //test
+ if(ruleMap[rule_opt[0]]==rule::centroid)
+ assert(nPointPolygon<=1);
for(int index=0;index<polyValues.size();++index){
double theValue=polyValues[index];
ostringstream fs;
--
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
