pktools  2.6.7
Processing Kernel for geospatial data
Filter2d.cc
1 /**********************************************************************
2 Filter2d.cc: class for filtering images
3 Copyright (C) 2008-2012 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #include <sstream>
21 #include <iomanip>
22 #include <iostream>
23 #include <cmath>
24 #include "Filter2d.h"
25 #include "StatFactory.h"
26 // #include "imageclasses/ImgUtils.h"
27 
28 filter2d::Filter2d::Filter2d(void)
29 {
30 }
31 
32 filter2d::Filter2d::Filter2d(const Vector2d<double> &taps)
33  : m_taps(taps)
34 {
35 }
36 
37 int filter2d::Filter2d::pushNoDataValue(double noDataValue)
38 {
39  if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
40  m_noDataValues.push_back(noDataValue);
41  return(m_noDataValues.size());
42 }
43 
44 void filter2d::Filter2d::setTaps(const Vector2d<double> &taps)
45 {
46  m_taps=taps;
47 }
48 
49 void filter2d::Filter2d::smoothNoData(ImgReaderGdal& input, ImgWriterGdal& output, int dim)
50 {
51  smoothNoData(input, output,dim,dim);
52 }
53 
54 void filter2d::Filter2d::smooth(ImgReaderGdal& input, ImgWriterGdal& output, int dim)
55 {
56  smooth(input, output,dim,dim);
57 }
58 
59 void filter2d::Filter2d::smoothNoData(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY)
60 {
61  m_taps.resize(dimY);
62  for(int j=0;j<dimY;++j){
63  m_taps[j].resize(dimX);
64  for(int i=0;i<dimX;++i)
65  m_taps[j][i]=1.0;
66  }
67  filter(input,output,false,true,true);
68 }
69 
70 void filter2d::Filter2d::smooth(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY)
71 {
72  m_taps.resize(dimY);
73  for(int j=0;j<dimY;++j){
74  m_taps[j].resize(dimX);
75  for(int i=0;i<dimX;++i)
76  m_taps[j][i]=1.0;
77  }
78  filter(input,output,false,true,false);
79 }
80 
81 
82 void filter2d::Filter2d::filter(ImgReaderGdal& input, ImgWriterGdal& output, bool absolute, bool normalize, bool noData)
83 {
84  int dimX=m_taps[0].size();//horizontal!!!
85  int dimY=m_taps.size();//vertical!!!
86  // byte* tmpbuf=new byte[input.rowSize()];
87  const char* pszMessage;
88  void* pProgressArg=NULL;
89  GDALProgressFunc pfnProgress=GDALTermProgress;
90  double progress=0;
91  pfnProgress(progress,pszMessage,pProgressArg);
92  for(int iband=0;iband<input.nrOfBand();++iband){
93  Vector2d<double> inBuffer(dimY,input.nrOfCol());
94  std::vector<double> outBuffer(input.nrOfCol());
95  int indexI=0;
96  int indexJ=0;
97  //initialize last half of inBuffer
98  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
99  try{
100  input.readData(inBuffer[indexJ],abs(j),iband);
101  }
102  catch(std::string errorstring){
103  std::cerr << errorstring << "in line " << indexJ << std::endl;
104  }
105  ++indexJ;
106  }
107 
108  for(int y=0;y<input.nrOfRow();++y){
109  if(y){//inBuffer already initialized for y=0
110  //erase first line from inBuffer
111  if(dimY>1)
112  inBuffer.erase(inBuffer.begin());
113  //read extra line and push back to inBuffer if not out of bounds
114  if(y+dimY/2<input.nrOfRow()){
115  //allocate buffer
116  if(dimY>1)
117  inBuffer.push_back(inBuffer.back());
118  try{
119  input.readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
120  }
121  catch(std::string errorstring){
122  std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
123  }
124  }
125  else{
126  int over=y+dimY/2-input.nrOfRow();
127  int index=(inBuffer.size()-1)-over;
128  assert(index>=0);
129  assert(index<inBuffer.size());
130  inBuffer.push_back(inBuffer[index]);
131  }
132  }
133  for(int x=0;x<input.nrOfCol();++x){
134  outBuffer[x]=0;
135  double norm=0;
136  bool masked=false;
137  if(noData){//only filter noData values
138  for(int imask=0;imask<m_noDataValues.size();++imask){
139  if(inBuffer[(dimY-1)/2][x]==m_noDataValues[imask]){
140  masked=true;
141  break;
142  }
143  }
144  if(!masked){
145  outBuffer[x]=inBuffer[(dimY-1)/2][x];
146  continue;
147  }
148  }
149  assert(!noData||masked);
150  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
151  for(int i=-(dimX-1)/2;i<=dimX/2;++i){
152  indexI=x+i;
153  indexJ=(dimY-1)/2+j;
154  //check if out of bounds
155  if(x<(dimX-1)/2)
156  indexI=x+abs(i);
157  else if(x>=input.nrOfCol()-(dimX-1)/2)
158  indexI=x-abs(i);
159  if(y<(dimY-1)/2)
160  indexJ=(dimY-1)/2+abs(j);
161  else if(y>=input.nrOfRow()-(dimY-1)/2)
162  indexJ=(dimY-1)/2-abs(j);
163  //do not take masked values into account
164  masked=false;
165  for(int imask=0;imask<m_noDataValues.size();++imask){
166  if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
167  masked=true;
168  break;
169  }
170  }
171  if(!masked){
172  outBuffer[x]+=(m_taps[(dimY-1)/2+j][(dimX-1)/2+i]*inBuffer[indexJ][indexI]);
173  norm+=m_taps[(dimY-1)/2+j][(dimX-1)/2+i];
174  }
175  }
176  }
177  if(absolute)
178  outBuffer[x]=(normalize&&norm)? abs(outBuffer[x])/norm : abs(outBuffer[x]);
179  else if(normalize&&norm!=0)
180  outBuffer[x]=outBuffer[x]/norm;
181  }
182  //write outBuffer to file
183  try{
184  output.writeData(outBuffer,y,iband);
185  }
186  catch(std::string errorstring){
187  std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
188  }
189  progress=(1.0+y);
190  progress+=(output.nrOfRow()*iband);
191  progress/=output.nrOfBand()*output.nrOfRow();
192  pfnProgress(progress,pszMessage,pProgressArg);
193  }
194  }
195 }
196 
197 
198 void filter2d::Filter2d::majorVoting(const std::string& inputFilename, const std::string& outputFilename,int dim,const std::vector<int> &prior)
199 {
200  const char* pszMessage;
201  void* pProgressArg=NULL;
202  GDALProgressFunc pfnProgress=GDALTermProgress;
203  double progress=0;
204  pfnProgress(progress,pszMessage,pProgressArg);
205 
206  bool usePriors=true;
207  if(prior.empty()){
208  std::cout << "no prior information" << std::endl;
209  usePriors=false;
210  }
211  else{
212  std::cout << "using priors ";
213  for(int iclass=0;iclass<prior.size();++iclass)
214  std::cout << " " << static_cast<short>(prior[iclass]);
215  std::cout << std::endl;
216  }
217 
218  ImgReaderGdal input;
219  ImgWriterGdal output;
220  input.open(inputFilename);
221  output.open(outputFilename,input);
222  int dimX=0;//horizontal!!!
223  int dimY=0;//vertical!!!
224  if(dim){
225  dimX=dim;
226  dimY=dim;
227  }
228  else{
229  dimX=m_taps[0].size();
230  dimY=m_taps.size();
231  }
232 
233  assert(dimX);
234  assert(dimY);
235 
236  Vector2d<double> inBuffer(dimY,input.nrOfCol());
237  std::vector<double> outBuffer(input.nrOfCol());
238  int indexI=0;
239  int indexJ=0;
240  //initialize last half of inBuffer
241  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
242  try{
243  input.readData(inBuffer[indexJ],abs(j));
244  }
245  catch(std::string errorstring){
246  std::cerr << errorstring << "in line " << indexJ << std::endl;
247  }
248  ++indexJ;
249  }
250 
251  for(int y=0;y<input.nrOfRow();++y){
252  if(y){//inBuffer already initialized for y=0
253  //erase first line from inBuffer
254  if(dimY>1)
255  inBuffer.erase(inBuffer.begin());
256  //read extra line and push back to inBuffer if not out of bounds
257  if(y+dimY/2<input.nrOfRow()){
258  //allocate buffer
259  if(dimY>1)
260  inBuffer.push_back(inBuffer.back());
261  try{
262  input.readData(inBuffer[inBuffer.size()-1],y+dimY/2);
263  }
264  catch(std::string errorstring){
265  std::cerr << errorstring << "in line" << y << std::endl;
266  }
267  }
268  else{
269  int over=y+dimY/2-input.nrOfRow();
270  int index=(inBuffer.size()-1)-over;
271  assert(index>=0);
272  assert(index<inBuffer.size());
273  inBuffer.push_back(inBuffer[index]);
274  }
275  }
276  for(int x=0;x<input.nrOfCol();++x){
277  outBuffer[x]=0;
278  std::map<int,int> occurrence;
279  int centre=dimX*(dimY-1)/2+(dimX-1)/2;
280  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
281  for(int i=-(dimX-1)/2;i<=dimX/2;++i){
282  indexI=x+i;
283  //check if out of bounds
284  if(indexI<0)
285  indexI=-indexI;
286  else if(indexI>=input.nrOfCol())
287  indexI=input.nrOfCol()-i;
288  if(y+j<0)
289  indexJ=-j;
290  else if(y+j>=input.nrOfRow())
291  indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
292  else
293  indexJ=(dimY-1)/2+j;
294 
295  // if(x<dimX/2)
296  // indexI=x+abs(i);
297  // else if(x>=input.nrOfCol()-dimX/2)
298  // indexI=x-abs(i);
299  // if(y<dimY/2)
300  // indexJ=dimY/2+abs(j);
301  // else if(y>=input.nrOfRow()-dimY/2)
302  // indexJ=dimY/2-abs(j);
303  if(usePriors){
304  occurrence[inBuffer[indexJ][indexI]]+=prior[inBuffer[indexJ][indexI]-1];
305  }
306  else
307  ++occurrence[inBuffer[indexJ][indexI]];
308  }
309  }
310  std::map<int,int>::const_iterator maxit=occurrence.begin();
311  for(std::map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
312  if(mit->second>maxit->second)
313  maxit=mit;
314  }
315  if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
316  outBuffer[x]=maxit->first;
317  else//favorize original value in case of ties
318  outBuffer[x]=inBuffer[(dimY-1)/2][x];
319  }
320  //write outBuffer to file
321  try{
322  output.writeData(outBuffer,y);
323  }
324  catch(std::string errorstring){
325  std::cerr << errorstring << "in line" << y << std::endl;
326  }
327  progress=(1.0+y)/output.nrOfRow();
328  pfnProgress(progress,pszMessage,pProgressArg);
329  }
330  input.close();
331  output.close();
332 }
333 
334 void filter2d::Filter2d::median(const std::string& inputFilename, const std::string& outputFilename,int dim, bool disc)
335 {
336  ImgReaderGdal input;
337  ImgWriterGdal output;
338  input.open(inputFilename);
339  output.open(outputFilename,input);
340  doit(input,output,"median",dim,disc);
341 }
342 
343 void filter2d::Filter2d::var(const std::string& inputFilename, const std::string& outputFilename,int dim, bool disc)
344 {
345  ImgReaderGdal input;
346  ImgWriterGdal output;
347  input.open(inputFilename);
348  output.open(outputFilename,input);
349  doit(input,output,"var",dim,disc);
350 }
351 
352 void filter2d::Filter2d::doit(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, bool disc)
353 {
354  doit(input,output,method,dim,dim,down,disc);
355 }
356 
357 void filter2d::Filter2d::doit(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, short down, bool disc)
358 {
359  const char* pszMessage;
360  void* pProgressArg=NULL;
361  GDALProgressFunc pfnProgress=GDALTermProgress;
362  double progress=0;
363  pfnProgress(progress,pszMessage,pProgressArg);
364 
365  assert(dimX);
366  assert(dimY);
367 
369  for(int iband=0;iband<input.nrOfBand();++iband){
370  Vector2d<double> inBuffer(dimY,input.nrOfCol());
371  std::vector<double> outBuffer((input.nrOfCol()+down-1)/down);
372  int indexI=0;
373  int indexJ=0;
374  //initialize last half of inBuffer
375  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
376  try{
377  input.readData(inBuffer[indexJ],abs(j),iband);
378  }
379  catch(std::string errorstring){
380  std::cerr << errorstring << "in line " << indexJ << std::endl;
381  }
382  ++indexJ;
383  }
384  for(int y=0;y<input.nrOfRow();++y){
385  if(y){//inBuffer already initialized for y=0
386  //erase first line from inBuffer
387  if(dimY>1)
388  inBuffer.erase(inBuffer.begin());
389  //read extra line and push back to inBuffer if not out of bounds
390  if(y+dimY/2<input.nrOfRow()){
391  //allocate buffer
392  if(dimY>1)
393  inBuffer.push_back(inBuffer.back());
394  try{
395  input.readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
396  }
397  catch(std::string errorstring){
398  std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
399  }
400  }
401  else{
402  int over=y+dimY/2-input.nrOfRow();
403  int index=(inBuffer.size()-1)-over;
404  assert(index>=0);
405  assert(index<inBuffer.size());
406  inBuffer.push_back(inBuffer[index]);
407  }
408  }
409  if((y+1+down/2)%down)
410  continue;
411  for(int x=0;x<input.nrOfCol();++x){
412  if((x+1+down/2)%down)
413  continue;
414  outBuffer[x/down]=0;
415  std::vector<double> windowBuffer;
416  std::map<long int,int> occurrence;
417  int centre=dimX*(dimY-1)/2+(dimX-1)/2;
418  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
419  for(int i=-(dimX-1)/2;i<=dimX/2;++i){
420  double d2=i*i+j*j;//square distance
421  if(disc&&(d2>(dimX/2)*(dimY/2)))
422  continue;
423  indexI=x+i;
424  //check if out of bounds
425  if(indexI<0)
426  indexI=-indexI;
427  else if(indexI>=input.nrOfCol())
428  indexI=input.nrOfCol()-i;
429  if(y+j<0)
430  indexJ=-j;
431  else if(y+j>=input.nrOfRow())
432  indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
433  else
434  indexJ=(dimY-1)/2+j;
435  bool masked=false;
436  for(int imask=0;imask<m_noDataValues.size();++imask){
437  if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
438  masked=true;
439  break;
440  }
441  }
442  if(!masked){
443  std::vector<short>::const_iterator vit=m_class.begin();
444  if(!m_class.size())
445  ++occurrence[inBuffer[indexJ][indexI]];
446  else{
447  while(vit!=m_class.end()){
448  if(inBuffer[indexJ][indexI]==*(vit++))
449  ++occurrence[inBuffer[indexJ][indexI]];
450  }
451  }
452  windowBuffer.push_back(inBuffer[indexJ][indexI]);
453  }
454  }
455  }
456  switch(getFilterType(method)){
457  case(filter2d::nvalid):
458  outBuffer[x/down]=stat.nvalid(windowBuffer);
459  break;
460  case(filter2d::median):
461  if(windowBuffer.empty())
462  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
463  else
464  outBuffer[x/down]=stat.median(windowBuffer);
465  break;
466  case(filter2d::var):{
467  if(windowBuffer.empty())
468  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
469  else
470  outBuffer[x/down]=stat.var(windowBuffer);
471  break;
472  }
473  case(filter2d::stdev):{
474  if(windowBuffer.empty())
475  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
476  else
477  outBuffer[x/down]=sqrt(stat.var(windowBuffer));
478  break;
479  }
480  case(filter2d::mean):{
481  if(windowBuffer.empty())
482  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
483  else
484  outBuffer[x/down]=stat.mean(windowBuffer);
485  break;
486  }
487  case(filter2d::min):{
488  if(windowBuffer.empty())
489  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
490  else
491  outBuffer[x/down]=stat.mymin(windowBuffer);
492  break;
493  }
494  case(filter2d::ismin):{
495  if(windowBuffer.empty())
496  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
497  else
498  outBuffer[x/down]=(stat.mymin(windowBuffer)==windowBuffer[centre])? 1:0;
499  break;
500  }
501  case(filter2d::minmax):{//is the same as homog?
502  double min=0;
503  double max=0;
504  if(windowBuffer.empty())
505  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
506  else{
507  stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
508  if(min!=max)
509  outBuffer[x/down]=0;
510  else
511  outBuffer[x/down]=windowBuffer[centre];//centre pixels
512  }
513  break;
514  }
515  case(filter2d::max):{
516  if(windowBuffer.empty())
517  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
518  else
519  outBuffer[x/down]=stat.mymax(windowBuffer);
520  break;
521  }
522  case(filter2d::ismax):{
523  if(windowBuffer.empty())
524  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
525  else
526  outBuffer[x/down]=(stat.mymax(windowBuffer)==windowBuffer[centre])? 1:0;
527  break;
528  }
529  case(filter2d::order):{
530  if(windowBuffer.empty())
531  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
532  else{
533  double lbound=0;
534  double ubound=dimX*dimY;
535  double theMin=stat.mymin(windowBuffer);
536  double theMax=stat.mymax(windowBuffer);
537  double scale=(ubound-lbound)/(theMax-theMin);
538  outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[centre]-theMin)+lbound);
539  }
540  break;
541  }
542  case(filter2d::sum):{
543  outBuffer[x/down]=stat.sum(windowBuffer);
544  break;
545  }
546  case(filter2d::percentile):{
547  assert(m_threshold.size());
548  outBuffer[x/down]=stat.percentile(windowBuffer,windowBuffer.begin(),windowBuffer.end(),m_threshold[0]);
549  break;
550  }
551  case(filter2d::proportion):{
552  if(windowBuffer.size()){
553  double sum=stat.sum(windowBuffer);
554  if(sum)
555  outBuffer[x/down]=100.0*windowBuffer[centre]/stat.sum(windowBuffer);
556  else
557  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
558  }
559  else
560  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
561  break;
562  }
563  case(filter2d::homog):
564  if(occurrence.size()==1)//all values in window are the same
565  outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
566  else
567  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
568  break;
569  case(filter2d::heterog):{
570  if(occurrence.size()==windowBuffer.size())
571  outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
572  else
573  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
574  // if(occurrence.size()==1)//all values in window are the same
575  // outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
576  // else
577  // outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
578  // break;
579  // for(std::vector<double>::const_iterator wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
580  // if(wit==windowBuffer.begin()+windowBuffer.size()/2)
581  // continue;
582  // else if(*wit!=inBuffer[(dimY-1)/2][x]){
583  // outBuffer[x/down]=1;
584  // break;
585  // }
586  // else if(*wit==inBuffer[(dimY-1)/2][x]){//todo:wit mag niet central pixel zijn
587  // outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
588  // break;
589  // }
590  // }
591  // break;
592  }
593  case(filter2d::density):{
594  if(windowBuffer.size()){
595  std::vector<short>::const_iterator vit=m_class.begin();
596  while(vit!=m_class.end())
597  outBuffer[x/down]+=100.0*occurrence[*(vit++)]/windowBuffer.size();
598  }
599  else
600  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
601  break;
602  }
603  case(filter2d::countid):{
604  if(windowBuffer.size())
605  outBuffer[x/down]=occurrence.size();
606  else
607  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
608  break;
609  }
610  case(filter2d::mode):{
611  if(occurrence.size()){
612  std::map<long int,int>::const_iterator maxit=occurrence.begin();
613  for(std::map<long int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
614  if(mit->second>maxit->second)
615  maxit=mit;
616  }
617  if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
618  outBuffer[x/down]=maxit->first;
619  else//favorize original value in case of ties
620  outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
621  }
622  else
623  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
624  break;
625  }
626  case(filter2d::threshold):{
627  assert(m_class.size()==m_threshold.size());
628  if(windowBuffer.size()){
629  outBuffer[x/down]=inBuffer[(dimY-1)/2][x];//initialize with original value (in case thresholds not met)
630  for(int iclass=0;iclass<m_class.size();++iclass){
631  if(100.0*(occurrence[m_class[iclass]])/windowBuffer.size()>m_threshold[iclass])
632  outBuffer[x/down]=m_class[iclass];
633  }
634  }
635  else
636  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
637  break;
638  }
639  case(filter2d::scramble):{//could be done more efficiently window by window with random shuffling entire buffer and assigning entire buffer at once to output image...
640  if(windowBuffer.size()){
641  int randomIndex=std::rand()%windowBuffer.size();
642  if(randomIndex>=windowBuffer.size())
643  outBuffer[x/down]=windowBuffer.back();
644  else if(randomIndex<0)
645  outBuffer[x/down]=windowBuffer[0];
646  else
647  outBuffer[x/down]=windowBuffer[randomIndex];
648  }
649  else
650  outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
651  break;
652  }
653  case(filter2d::mixed):{
654  enum Type { BF=11, CF=12, MF=13, NF=20, W=30 };
655  double nBF=occurrence[BF];
656  double nCF=occurrence[CF];
657  double nMF=occurrence[MF];
658  double nNF=occurrence[NF];
659  double nW=occurrence[W];
660  if(windowBuffer.size()){
661  if((nBF+nCF+nMF)&&(nBF+nCF+nMF>=nNF+nW)){//forest
662  if(nBF/(nBF+nCF)>=0.75)
663  outBuffer[x/down]=BF;
664  else if(nCF/(nBF+nCF)>=0.75)
665  outBuffer[x/down]=CF;
666  else
667  outBuffer[x/down]=MF;
668  }
669  else{//non-forest
670  if(nW&&(nW>=nNF))
671  outBuffer[x/down]=W;
672  else
673  outBuffer[x/down]=NF;
674  }
675  }
676  else
677  outBuffer[x/down]=inBuffer[indexJ][indexI];
678  break;
679  }
680  default:{
681  std::ostringstream ess;
682  ess << "Error: filter method " << method << " not supported" << std::endl;
683  throw(ess.str());
684  break;
685  }
686  }
687  }
688  progress=(1.0+y/down);
689  progress+=(output.nrOfRow()*iband);
690  progress/=output.nrOfBand()*output.nrOfRow();
691  pfnProgress(progress,pszMessage,pProgressArg);
692  //write outBuffer to file
693  try{
694  output.writeData(outBuffer,y/down,iband);
695  }
696  catch(std::string errorstring){
697  std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
698  }
699  }
700  }
701  pfnProgress(1.0,pszMessage,pProgressArg);
702 }
703 
704 void filter2d::Filter2d::mrf(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, double beta, bool eightConnectivity, short down, bool verbose){
705  assert(m_class.size()>1);
706  Vector2d<double> fullBeta(m_class.size(),m_class.size());
707  for(int iclass1=0;iclass1<m_class.size();++iclass1)
708  for(int iclass2=0;iclass2<m_class.size();++iclass2)
709  fullBeta[iclass1][iclass2]=beta;
710  mrf(input,output,dimX,dimY,fullBeta,eightConnectivity,down,verbose);
711 }
712 
713 //beta[classTo][classFrom]
714 void filter2d::Filter2d::mrf(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, Vector2d<double> beta, bool eightConnectivity, short down, bool verbose)
715 {
716  const char* pszMessage;
717  void* pProgressArg=NULL;
718  GDALProgressFunc pfnProgress=GDALTermProgress;
719  double progress=0;
720  pfnProgress(progress,pszMessage,pProgressArg);
721 
722  assert(dimX);
723  assert(dimY);
724 
725  Vector2d<short> inBuffer(dimY,input.nrOfCol());
726  Vector2d<double> outBuffer(m_class.size(),(input.nrOfCol()+down-1)/down);
727  assert(input.nrOfBand()==1);
728  assert(output.nrOfBand()==m_class.size());
729  assert(m_class.size()>1);
730  assert(beta.size()==m_class.size());
731  int indexI=0;
732  int indexJ=0;
733  //initialize last half of inBuffer
734  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
735  try{
736  input.readData(inBuffer[indexJ],abs(j));
737  }
738  catch(std::string errorstring){
739  std::cerr << errorstring << "in line " << indexJ << std::endl;
740  }
741  ++indexJ;
742  }
743  for(int y=0;y<input.nrOfRow();++y){
744  if(y){//inBuffer already initialized for y=0
745  //erase first line from inBuffer
746  if(dimY>1)
747  inBuffer.erase(inBuffer.begin());
748  //read extra line and push back to inBuffer if not out of bounds
749  if(y+dimY/2<input.nrOfRow()){
750  //allocate buffer
751  if(dimY>1)
752  inBuffer.push_back(inBuffer.back());
753  try{
754  input.readData(inBuffer[inBuffer.size()-1],y+dimY/2);
755  }
756  catch(std::string errorstring){
757  std::cerr << errorstring << "in line " << y << std::endl;
758  }
759  }
760  else{
761  int over=y+dimY/2-input.nrOfRow();
762  int index=(inBuffer.size()-1)-over;
763  assert(index>=0);
764  assert(index<inBuffer.size());
765  inBuffer.push_back(inBuffer[index]);
766  }
767  }
768  if((y+1+down/2)%down)
769  continue;
770  for(int x=0;x<input.nrOfCol();++x){
771  if((x+1+down/2)%down)
772  continue;
773  std::vector<short> potential(m_class.size());
774  for(int iclass=0;iclass<m_class.size();++iclass){
775  potential[iclass]=0;
776  outBuffer[iclass][x/down]=0;
777  }
778  std::vector<double> windowBuffer;
779  int centre=dimX*(dimY-1)/2+(dimX-1)/2;
780  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
781  for(int i=-(dimX-1)/2;i<=dimX/2;++i){
782  if(i!=0&&j!=0&&!eightConnectivity)
783  continue;
784  if(i==0&&j==0)
785  continue;
786  indexI=x+i;
787  //check if out of bounds
788  if(indexI<0)
789  indexI=-indexI;
790  else if(indexI>=input.nrOfCol())
791  indexI=input.nrOfCol()-i;
792  if(y+j<0)
793  indexJ=-j;
794  else if(y+j>=input.nrOfRow())
795  indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
796  else
797  indexJ=(dimY-1)/2+j;
798  bool masked=false;
799  for(int imask=0;imask<m_noDataValues.size();++imask){
800  if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
801  masked=true;
802  break;
803  }
804  }
805  if(!masked){
806  for(int iclass=0;iclass<m_class.size();++iclass){
807  if(inBuffer[indexJ][indexI]==m_class[iclass])
808  potential[iclass]+=1;
809  }
810  }
811  }
812  }
813  double norm=0;
814  for(int iclass1=0;iclass1<m_class.size();++iclass1){
815  assert(beta[iclass1].size()==m_class.size());
816  double pot=0;
817  for(int iclass2=0;iclass2<m_class.size();++iclass2)
818  if(iclass2!=iclass1)
819  pot+=potential[iclass2]*beta[iclass1][iclass2];
820  double prior=exp(-pot);
821  outBuffer[iclass1][x/down]=prior;
822  norm+=prior;
823  }
824  if(norm){
825  for(int iclass1=0;iclass1<m_class.size();++iclass1)
826  outBuffer[iclass1][x/down]/=norm;
827  }
828  }
829  progress=(1.0+y/down)/output.nrOfRow();
830  pfnProgress(progress,pszMessage,pProgressArg);
831  //write outBuffer to file
832  assert(outBuffer.size()==m_class.size());
833  assert(y<output.nrOfRow());
834  for(int iclass=0;iclass<m_class.size();++iclass){
835  assert(outBuffer[iclass].size()==output.nrOfCol());
836  try{
837  output.writeData(outBuffer[iclass],y/down,iclass);
838  }
839  catch(std::string errorstring){
840  std::cerr << errorstring << "in class " << iclass << ", line " << y << std::endl;
841  }
842  }
843  }
844 }
845 
846 void filter2d::Filter2d::shift(ImgReaderGdal& input, ImgWriterGdal& output, double offsetX, double offsetY, double randomSigma, RESAMPLE resample, bool verbose)
847 {
848  assert(input.nrOfCol()==output.nrOfCol());
849  assert(input.nrOfRow()==output.nrOfRow());
850  assert(input.nrOfBand()==output.nrOfBand());
851  const char* pszMessage;
852  void* pProgressArg=NULL;
853  GDALProgressFunc pfnProgress=GDALTermProgress;
854  double progress=0;
855  pfnProgress(progress,pszMessage,pProgressArg);
856  //process band per band in memory
857  Vector2d<double> inBuffer(input.nrOfRow(),output.nrOfCol());
858  Vector2d<double> outBuffer(input.nrOfRow(),output.nrOfCol());
859  for(int iband=0;iband<input.nrOfBand();++iband){
860  input.readDataBlock(inBuffer,0,inBuffer.nCols()-1,0,inBuffer.nRows()-1,iband);
861  shift(inBuffer,outBuffer,offsetX,offsetY,randomSigma,resample,verbose);
862  output.writeDataBlock(outBuffer,0,outBuffer.nCols()-1,0,outBuffer.nRows()-1,iband);
863  }
864 }
865 
866 //todo: re-implement without dependency of CImg and reg libraries
867 // void filter2d::Filter2d::dwt_texture(const std::string& inputFilename, const std::string& outputFilename,int dim, int scale, int down, int iband, bool verbose)
868 // {
869 // ImgReaderGdal input;
870 // ImgWriterGdal output;
871 // if(verbose)
872 // std::cout << "opening file " << inputFilename << std::endl;
873 // input.open(inputFilename);
874 // double magicX=1,magicY=1;
875 // output.open(outputFilename,(input.nrOfCol()+down-1)/down,(input.nrOfRow()+down-1)/down,scale*3,GDT_Float32,input.getImageType());
876 // if(input.isGeoRef()){
877 // output.setProjection(input.getProjection());
878 // output.copyGeoTransform(input);
879 // }
880 // if(verbose)
881 // std::cout << "Dimension texture (row x col x band) = " << (input.nrOfCol()+down-1)/down << " x " << (input.nrOfRow()+down-1)/down << " x " << scale*3 << std::endl;
882 // assert(dim%2);
883 // int dimX=dim;
884 // int dimY=dim;
885 // Vector2d<float> inBuffer(dimY,input.nrOfCol());
886 // Vector2d<float> outBuffer(scale*3,(input.nrOfCol()+down-1)/down);
887 // //initialize last half of inBuffer
888 // int indexI=0;
889 // int indexJ=0;
890 // for(int j=-dimY/2;j<(dimY+1)/2;++j){
891 // try{
892 // if(verbose)
893 // cout << "reading input line " << abs(j) << std::endl;
894 // input.readData(inBuffer[indexJ],GDT_Float32,abs(j),iband);
895 // ++indexJ;
896 // }
897 // catch(std::string errorstring){
898 // std::cerr << errorstring << "in band " << iband << ", line " << indexJ << std::endl;
899 // }
900 // }
901 // const char* pszMessage;
902 // void* pProgressArg=NULL;
903 // GDALProgressFunc pfnProgress=GDALTermProgress;
904 // double progress=0;
905 // pfnProgress(progress,pszMessage,pProgressArg);
906 // for(int y=0;y<input.nrOfRow();y+=down){
907 // if(verbose)
908 // std::cout << "calculating line " << y/down << std::endl;
909 // if(y){//inBuffer already initialized for y=0
910 // //erase first line from inBuffer
911 // inBuffer.erase(inBuffer.begin());
912 // //read extra line and push back to inBuffer if not out of bounds
913 // if(y+dimY/2<input.nrOfRow()){
914 // //allocate buffer
915 // inBuffer.push_back(inBuffer.back());
916 // try{
917 // if(verbose)
918 // std::cout << "reading input line " << y+dimY/2 << std::endl;
919 // input.readData(inBuffer[inBuffer.size()-1],GDT_Float32,y+dimY/2,iband);
920 // }
921 // catch(std::string errorstring){
922 // std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
923 // }
924 // }
925 // }
926 // for(int x=0;x<input.nrOfCol();x+=down){
927 // Vector2d<double> texture_feature(scale,3);
928 // CImg<> texture_in(dimX,dimY);
929 // int r=0;//index for row of texture_in
930 // for(int j=-dimY/2;j<(dimY+1)/2;++j){
931 // int c=0;
932 // for(int i=-dimX/2;i<(dimX+1)/2;++i){
933 // indexI=x+i;
934 // //check if out of bounds
935 // if(indexI<0)
936 // indexI=-indexI;
937 // else if(indexI>=input.nrOfCol())
938 // indexI=input.nrOfCol()-i;
939 // if(y+j<0)
940 // indexJ=-j;
941 // else if(y+j>=input.nrOfRow())
942 // indexJ=dimY/2-j;//indexJ=inBuffer.size()-1-j;
943 // else
944 // indexJ=dimY/2+j;
945 // assert(indexJ<inBuffer.size());
946 // assert(indexI<inBuffer[indexJ].size());
947 // texture_in(r,c)=inBuffer[indexJ][indexI];
948 // c++;
949 // }
950 // ++r;
951 // }
952 // texture_in.dwt_texture(texture_feature,scale);
953 // for(int v=0;v<scale*3;++v)
954 // outBuffer[v][x/down]=texture_feature[v/3][v%3];
955 // }
956 // //write outBuffer to file
957 // try{
958 // if(verbose)
959 // std::cout << "writing line " << y/down << std::endl;
960 // for(int v=0;v<scale*3;++v)
961 // output.writeData(outBuffer[v],GDT_Float32,y/down,v);
962 // }
963 // catch(std::string errorstring){
964 // std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
965 // }
966 // progress=(1.0+y)/output.nrOfRow();
967 // pfnProgress(progress,pszMessage,pProgressArg);
968 // }
969 // input.close();
970 // output.close();
971 // }
972 
973 void filter2d::Filter2d::morphology(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, const std::vector<double> &angle, bool disc)
974 {
975  const char* pszMessage;
976  void* pProgressArg=NULL;
977  GDALProgressFunc pfnProgress=GDALTermProgress;
978  double progress=0;
979  pfnProgress(progress,pszMessage,pProgressArg);
980 
981  assert(dimX);
982  assert(dimY);
983 
985  for(int iband=0;iband<input.nrOfBand();++iband){
986  Vector2d<double> inBuffer(dimY,input.nrOfCol());
987  std::vector<double> outBuffer(input.nrOfCol());
988  int indexI=0;
989  int indexJ=0;
990  //initialize last half of inBuffer
991  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
992  try{
993  input.readData(inBuffer[indexJ],abs(j),iband);
994  ++indexJ;
995  }
996  catch(std::string errorstring){
997  std::cerr << errorstring << "in line " << indexJ << std::endl;
998  }
999  }
1000  for(int y=0;y<input.nrOfRow();++y){
1001  if(y){//inBuffer already initialized for y=0
1002  //erase first line from inBuffer
1003  if(dimY>1)
1004  inBuffer.erase(inBuffer.begin());
1005  //read extra line and push back to inBuffer if not out of bounds
1006  if(y+dimY/2<input.nrOfRow()){
1007  //allocate buffer
1008  if(dimY>1)
1009  inBuffer.push_back(inBuffer.back());
1010  try{
1011  input.readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
1012  }
1013  catch(std::string errorstring){
1014  std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
1015  }
1016  }
1017  else{
1018  int over=y+dimY/2-input.nrOfRow();
1019  int index=(inBuffer.size()-1)-over;
1020  assert(index>=0);
1021  assert(index<inBuffer.size());
1022  inBuffer.push_back(inBuffer[index]);
1023  }
1024  }
1025  for(int x=0;x<input.nrOfCol();++x){
1026  double currentValue=inBuffer[(dimY-1)/2][x];
1027  outBuffer[x]=currentValue;
1028  std::vector<double> statBuffer;
1029  bool currentMasked=false;
1030  int centre=dimX*(dimY-1)/2+(dimX-1)/2;
1031  for(int imask=0;imask<m_noDataValues.size();++imask){
1032  if(currentValue==m_noDataValues[imask]){
1033  currentMasked=true;
1034  break;
1035  }
1036  }
1037  if(currentMasked){
1038  outBuffer[x]=currentValue;
1039  }
1040  else{
1041  for(int j=-(dimY-1)/2;j<=dimY/2;++j){
1042  for(int i=-(dimX-1)/2;i<=dimX/2;++i){
1043  double d2=i*i+j*j;//square distance
1044  if(disc&&(d2>(dimX/2)*(dimY/2)))
1045  continue;
1046  if(angle.size()){
1047  double theta;
1048  //use polar coordinates in radians
1049  if(i>0){
1050  if(j<0)
1051  theta=atan(static_cast<double>(-j)/(static_cast<double>(i)));
1052  else
1053  theta=-atan(static_cast<double>(j)/(static_cast<double>(i)));
1054  }
1055  else if(i<0){
1056  if(j<0)
1057  theta=PI-atan(static_cast<double>(-j)/(static_cast<double>(-i)));
1058  else
1059  theta=PI+atan(static_cast<double>(j)/(static_cast<double>(-i)));
1060  }
1061  else if(j<0)
1062  theta=PI/2.0;
1063  else if(j>0)
1064  theta=3.0*PI/2.0;
1065  //convert to North (0), East (90), South (180), West (270) in degrees
1066  theta=360-(theta/PI*180)+90;
1067  if(theta<0)
1068  theta+=360;
1069  while(theta>360)
1070  theta-=360;
1071  bool alligned=false;
1072  for(int iangle=0;iangle<angle.size();++iangle){
1073  if(sqrt((theta-angle[iangle])*(theta-angle[iangle]))<10){
1074  alligned=true;
1075  break;
1076  }
1077  }
1078  if(!alligned)
1079  continue;
1080  }
1081  indexI=x+i;
1082  //check if out of bounds
1083  if(indexI<0)
1084  indexI=-indexI;
1085  else if(indexI>=input.nrOfCol())
1086  indexI=input.nrOfCol()-i;
1087  if(y+j<0)
1088  indexJ=-j;
1089  else if(y+j>=input.nrOfRow())
1090  indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
1091  else
1092  indexJ=(dimY-1)/2+j;
1093  //todo: introduce novalue as this: ?
1094  // if(inBuffer[indexJ][indexI]==(m_noDataValues.size())? m_noDataValues[0] : 0)
1095  // continue;
1096  bool masked=false;
1097  for(int imask=0;imask<m_noDataValues.size();++imask){
1098  if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
1099  masked=true;
1100  break;
1101  }
1102  }
1103  if(!masked){
1104  short binValue=0;
1105  for(int iclass=0;iclass<m_class.size();++iclass){
1106  if(inBuffer[indexJ][indexI]==m_class[iclass]){
1107  binValue=1;
1108  break;
1109  }
1110  }
1111  if(m_class.size())
1112  statBuffer.push_back(binValue);
1113  else
1114  statBuffer.push_back(inBuffer[indexJ][indexI]);
1115  }
1116  }
1117  }
1118  if(statBuffer.size()){
1119  switch(getFilterType(method)){
1120  case(filter2d::dilate):
1121  outBuffer[x]=stat.mymax(statBuffer);
1122  break;
1123  case(filter2d::erode):
1124  outBuffer[x]=stat.mymin(statBuffer);
1125  break;
1126  default:
1127  std::ostringstream ess;
1128  ess << "Error: morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << std::endl;
1129  throw(ess.str());
1130  break;
1131  }
1132  }
1133  if(outBuffer[x]&&m_class.size())
1134  outBuffer[x]=m_class[0];
1135  }
1136  }
1137  //write outBuffer to file
1138  try{
1139  output.writeData(outBuffer,y,iband);
1140  }
1141  catch(std::string errorstring){
1142  std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
1143  }
1144  progress=(1.0+y);
1145  progress+=(output.nrOfRow()*iband);
1146  progress/=output.nrOfBand()*output.nrOfRow();
1147  pfnProgress(progress,pszMessage,pProgressArg);
1148  }
1149  }
1150 }
1151 
1152 void filter2d::Filter2d::shadowDsm(ImgReaderGdal& input, ImgWriterGdal& output, double sza, double saa, double pixelSize, short shadowFlag){
1153  Vector2d<float> inputBuffer;
1154  Vector2d<float> outputBuffer;
1155  input.readDataBlock(inputBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, 0);
1156  shadowDsm(inputBuffer, outputBuffer, sza, saa, pixelSize, shadowFlag);
1157  output.writeDataBlock(outputBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,0);
1158 }
1159 
1160 void filter2d::Filter2d::dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
1161  Vector2d<float> theBuffer;
1162  for(int iband=0;iband<input.nrOfBand();++iband){
1163  input.readDataBlock(theBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
1164  std::cout << "filtering band " << iband << std::endl << std::flush;
1165  dwtForward(theBuffer, wavelet_type, family);
1166  output.writeDataBlock(theBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1167  }
1168 }
1169 
1170 void filter2d::Filter2d::dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
1171  Vector2d<float> theBuffer;
1172  for(int iband=0;iband<input.nrOfBand();++iband){
1173  input.readDataBlock(theBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
1174  std::cout << "filtering band " << iband << std::endl << std::flush;
1175  dwtInverse(theBuffer, wavelet_type, family);
1176  output.writeDataBlock(theBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1177  }
1178 }
1179 
1180 void filter2d::Filter2d::dwtCut(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut, bool verbose){
1181  Vector2d<float> theBuffer;
1182  for(int iband=0;iband<input.nrOfBand();++iband){
1183  input.readDataBlock(theBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
1184  std::cout << "filtering band " << iband << std::endl << std::flush;
1185  dwtCut(theBuffer, wavelet_type, family, cut);
1186  output.writeDataBlock(theBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1187  }
1188 }
1189 
1190 void filter2d::Filter2d::linearFeature(ImgReaderGdal& input, ImgWriterGdal& output, float angle, float angleStep, float maxDistance, float eps, bool l1, bool a1, bool l2, bool a2, int band, bool verbose){
1191  Vector2d<float> inputBuffer;
1192  std::vector< Vector2d<float> > outputBuffer;
1193  input.readDataBlock(inputBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, band);
1194  if(maxDistance<=0)
1195  maxDistance=sqrt(static_cast<float>(input.nrOfCol()*input.nrOfRow()));
1196  linearFeature(inputBuffer,outputBuffer,angle,angleStep,maxDistance,eps, l1, a1, l2, a2,verbose);
1197  for(int iband=0;iband<outputBuffer.size();++iband)
1198  output.writeDataBlock(outputBuffer[iband],0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1199 }
1200 
1201 void filter2d::Filter2d::linearFeature(const Vector2d<float>& input, std::vector< Vector2d<float> >& output, float angle, float angleStep, float maxDistance, float eps, bool l1, bool a1, bool l2, bool a2, bool verbose)
1202 {
1203  output.clear();
1204  int nband=0;//linear feature
1205  if(l1)
1206  ++nband;
1207  if(a1)
1208  ++nband;
1209  if(l2)
1210  ++nband;
1211  if(a2)
1212  ++nband;
1213  output.resize(nband);
1214  for(int iband=0;iband<output.size();++iband)
1215  output[iband].resize(input.nRows(),input.nCols());
1216  if(maxDistance<=0)
1217  maxDistance=sqrt(static_cast<float>(input.nRows()*input.nCols()));
1218  int indexI=0;
1219  int indexJ=0;
1220  const char* pszMessage;
1221  void* pProgressArg=NULL;
1222  GDALProgressFunc pfnProgress=GDALTermProgress;
1223  double progress=0;
1224  pfnProgress(progress,pszMessage,pProgressArg);
1225  for(int y=0;y<input.nRows();++y){
1226  for(int x=0;x<input.nCols();++x){
1227  float currentValue=input[y][x];
1228  //find values equal to current value with some error margin
1229  //todo: add distance for two opposite directions
1230  float lineDistance1=0;//longest line of object
1231  float lineDistance2=maxDistance;//shortest line of object
1232  float lineAngle1=0;//angle to longest line (North=0)
1233  float lineAngle2=0;//angle to shortest line (North=0)
1234  float northAngle=0;//rotating angle
1235  for(northAngle=0;northAngle<180;northAngle+=angleStep){
1236  if(angle<=360&&angle>=0&&angle!=northAngle)
1237  continue;
1238  //test
1239  if(verbose)
1240  std::cout << "northAngle: " << northAngle << std::endl;
1241  float currentDistance=0;
1242  float theDir=0;
1243  for(short side=0;side<=1;side+=1){
1244  theDir=PI/2.0-DEG2RAD(northAngle)+side*PI;//in radians
1245  //test
1246  if(verbose)
1247  std::cout << "theDir in deg: " << RAD2DEG(theDir) << std::endl;
1248  if(theDir<0)
1249  theDir+=2*PI;
1250  //test
1251  if(verbose)
1252  std::cout << "theDir in deg: " << RAD2DEG(theDir) << std::endl;
1253  float nextValue=currentValue;
1254  for(float currentRay=1;currentRay<maxDistance;++currentRay){
1255  indexI=x+currentRay*cos(theDir);
1256  indexJ=y-currentRay*sin(theDir);
1257  if(indexJ<0||indexJ>=input.size())
1258  break;
1259  if(indexI<0||indexI>=input[indexJ].size())
1260  break;
1261  nextValue=input[indexJ][indexI];
1262  if(verbose){
1263  std::cout << "x: " << x << std::endl;
1264  std::cout << "y: " << y << std::endl;
1265  std::cout << "currentValue: " << currentValue << std::endl;
1266  std::cout << "theDir in degrees: " << RAD2DEG(theDir) << std::endl;
1267  std::cout << "cos(theDir): " << cos(theDir) << std::endl;
1268  std::cout << "sin(theDir): " << sin(theDir) << std::endl;
1269  std::cout << "currentRay: " << currentRay << std::endl;
1270  std::cout << "currentDistance: " << currentDistance << std::endl;
1271  std::cout << "indexI: " << indexI << std::endl;
1272  std::cout << "indexJ: " << indexJ << std::endl;
1273  std::cout << "nextValue: " << nextValue << std::endl;
1274  }
1275  if(fabs(currentValue-nextValue)<=eps){
1276  ++currentDistance;
1277  //test
1278  if(verbose)
1279  std::cout << "currentDistance: " << currentDistance << ", continue" << std::endl;
1280  }
1281  else{
1282  if(verbose)
1283  std::cout << "currentDistance: " << currentDistance << ", break" << std::endl;
1284  break;
1285  }
1286  }
1287  }
1288  if(lineDistance1<currentDistance){
1289  lineDistance1=currentDistance;
1290  lineAngle1=northAngle;
1291  }
1292  if(lineDistance2>currentDistance){
1293  lineDistance2=currentDistance;
1294  lineAngle2=northAngle;
1295  }
1296  if(verbose){
1297  std::cout << "lineDistance1: " << lineDistance1 << std::endl;
1298  std::cout << "lineAngle1: " << lineAngle1 << std::endl;
1299  std::cout << "lineDistance2: " << lineDistance2 << std::endl;
1300  std::cout << "lineAngle2: " << lineAngle2 << std::endl;
1301  }
1302  }
1303  int iband=0;
1304  if(l1)
1305  output[iband++][y][x]=lineDistance1;
1306  if(a1)
1307  output[iband++][y][x]=lineAngle1;
1308  if(l2)
1309  output[iband++][y][x]=lineDistance2;
1310  if(a2)
1311  output[iband++][y][x]=lineAngle2;
1312  assert(iband==nband);
1313  }
1314  progress=(1.0+y);
1315  progress/=input.nRows();
1316  pfnProgress(progress,pszMessage,pProgressArg);
1317  }
1318 }
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: Filter.h:33
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96
bool writeDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Write pixel cell values for a range of columns and rows for a specific band (all indices start counti...
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
void close(void)
Close the image.
void readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98