260{
261 const char* pszMessage;
262 void* pProgressArg=NULL;
263 GDALProgressFunc pfnProgress=GDALTermProgress;
264 double progress=0;
265 pfnProgress(progress,pszMessage,pProgressArg);
266
268 double noDataValue=0;
269 if(m_noDataValues.size()){
270 stat.setNoDataValues(m_noDataValues);
271 noDataValue=m_noDataValues[0];
272 }
273 assert(dimX);
274 assert(dimY);
275
276 outputVector.resize((inputVector.size()+down-1)/down);
278 std::vector<T2> outBuffer((inputVector[0].size()+down-1)/down);
279
280 int indexI=0;
281 int indexJ=0;
282
283 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
284 inBuffer[indexJ]=inputVector[abs(j)];
285 ++indexJ;
286 }
287 for(int y=0;y<inputVector.size();++y){
288
289 if(y){
290
291 inBuffer.erase(inBuffer.begin());
292
293 if(y+dimY/2<inputVector.size())
294 inBuffer.push_back(inputVector[y+dimY/2]);
295 else{
296 int over=y+dimY/2-inputVector.size();
297 int index=(inBuffer.size()-1)-over;
298 assert(index>=0);
299 assert(index<inBuffer.size());
300 inBuffer.push_back(inBuffer[index]);
301 }
302 }
303 if((y+1+down/2)%down)
304 continue;
305 for(int x=0;x<inputVector[0].size();++x){
306 if((x+1+down/2)%down)
307 continue;
308 outBuffer[x/down]=0;
309
310 std::vector<T1> windowBuffer;
311 std::map<int,int> occurrence;
312 int centre=dimX*(dimY-1)/2+(dimX-1)/2;
313 bool centreMasked=false;
314 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
315 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
316 indexI=x+i;
317
318 if(indexI<0)
319 indexI=-indexI;
320 else if(indexI>=inputVector[0].size())
321 indexI=inputVector[0].size()-i;
322 if(y+j<0)
323 indexJ=-j;
324 else if(y+j>=inputVector.size())
325 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
326 else
327 indexJ=(dimY-1)/2+j;
328 windowBuffer.push_back(inBuffer[indexJ][indexI]);
329 if(!stat.isNoData(inBuffer[indexJ][indexI])){
330 std::vector<short>::const_iterator vit=m_class.begin();
331
332 if(!m_class.size())
333 ++occurrence[inBuffer[indexJ][indexI]];
334 else{
335 while(vit!=m_class.end()){
336 if(inBuffer[indexJ][indexI]==*(vit++))
337 ++occurrence[inBuffer[indexJ][indexI]];
338 }
339 }
340 }
341 }
342 }
343 switch(getFilterType(method)){
344 case(filter2d::nvalid):
345 outBuffer[x/down]=stat.nvalid(windowBuffer);
346 break;
347 case(filter2d::median):
348 outBuffer[x/down]=stat.median(windowBuffer);
349 break;
350 case(filter2d::var):{
351 outBuffer[x/down]=stat.var(windowBuffer);
352 break;
353 }
354 case(filter2d::stdev):{
355 T2 varValue=stat.var(windowBuffer);
356 if(stat.isNoData(varValue))
357 outBuffer[x/down]=noDataValue;
358 else
359 outBuffer[x/down]=sqrt(varValue);
360 break;
361 }
362 case(filter2d::mean):{
363 if(windowBuffer.empty())
364 outBuffer[x/down]=noDataValue;
365 else
366 outBuffer[x/down]=stat.mean(windowBuffer);
367 break;
368 }
369 case(filter2d::min):{
370 outBuffer[x/down]=stat.mymin(windowBuffer);
371 break;
372 }
373 case(filter2d::ismin):{
374 T1 minValue=stat.mymin(windowBuffer);
375 if(stat.isNoData(minValue))
376 outBuffer[x/down]=noDataValue;
377 else
378 outBuffer[x/down]=(windowBuffer[centre]==minValue)? 1:0;
379 break;
380 }
381 case(filter2d::minmax):{
382 T1 min=0;
383 T1 max=0;
384 stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
385 if(min!=max)
386 outBuffer[x/down]=0;
387 else
388 outBuffer[x/down]=windowBuffer[centre];
389 break;
390 }
391 case(filter2d::max):{
392 outBuffer[x/down]=stat.mymax(windowBuffer);
393 break;
394 }
395 case(filter2d::ismax):{
396 T1 maxValue=stat.mymax(windowBuffer);
397 if(stat.isNoData(maxValue))
398 outBuffer[x/down]=noDataValue;
399 else
400 outBuffer[x/down]=(windowBuffer[centre]==maxValue)? 1:0;
401 break;
402 }
403 case(filter2d::order):{
404 stat.eraseNoData(windowBuffer);
405 if(windowBuffer.empty())
406 outBuffer[x/down]=noDataValue;
407 else{
408 double lbound=0;
409 double ubound=dimX*dimY;
410 double theMin=stat.mymin(windowBuffer);
411 double theMax=stat.mymax(windowBuffer);
412 double scale=(ubound-lbound)/(theMax-theMin);
413 outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[centre]-theMin)+lbound);
414 }
415 break;
416 }
417 case(filter2d::sum):{
418 outBuffer[x/down]=stat.sum(windowBuffer);
419 break;
420 }
421 case(filter2d::percentile):{
422 assert(m_threshold.size());
423 outBuffer[x/down]=stat.percentile(windowBuffer,windowBuffer.begin(),windowBuffer.end(),m_threshold[0]);
424 break;
425 }
426 case(filter2d::proportion):{
427 stat.eraseNoData(windowBuffer);
428 T2 sum=stat.sum(windowBuffer);
429 if(sum)
430 outBuffer[x/down]=windowBuffer[centre]/sum;
431 else
432 outBuffer[x/down]=noDataValue;
433 break;
434 }
435 case(filter2d::homog):{
436 T1 centreValue=inBuffer[(dimY-1)/2][x];
437 bool isHomog=true;
438 stat.eraseNoData(windowBuffer);
439 typename std::vector<T1>::const_iterator wit;
440 for(wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
441 if(*wit==centreValue)
442 continue;
443 else{
444 isHomog=false;
445 break;
446 }
447 }
448 if(isHomog)
449 outBuffer[x/down]=1;
450 else
451 outBuffer[x/down]=noDataValue;
452 break;
453 }
454 case(filter2d::heterog):{
455 T1 centreValue=inBuffer[(dimY-1)/2][x];
456 bool isHeterog=true;
457 stat.eraseNoData(windowBuffer);
458 typename std::vector<T1>::const_iterator wit;
459 for(wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
460 if(*wit!=centreValue)
461 continue;
462 else{
463 isHeterog=false;
464 break;
465 }
466 }
467 if(isHeterog)
468 outBuffer[x/down]=1;
469 else
470 outBuffer[x/down]=noDataValue;
471 break;
472 }
473 case(filter2d::sauvola):{
474 try{
475 double theMean=0;
476 double theStdev=0;
477 bool invalid=false;
478 T1 centreValue=inBuffer[(dimY-1)/2][x];
479 if(windowBuffer.empty()||stat.isNoData(centreValue)){
480 invalid=true;
481 throw(invalid);
482 }
483 stat.meanVar(windowBuffer,theMean,theStdev);
484 theStdev=sqrt(theStdev);
485 double kValue=0.5;
486 double rValue=128;
487 if(m_threshold.size()==2){
488 kValue=m_threshold[0];
489 rValue=m_threshold[1];
490 }
491
492
493 double theThreshold=theMean*(1+kValue*(theStdev/rValue - 1));
494
495 outBuffer[x/down]=(centreValue>theThreshold) ? 1 : noDataValue;
496 }
497 catch(bool invalid){
498 outBuffer[x/down]=noDataValue;
499 }
500 break;
501 }
502 case(filter2d::density):{
503 int nvalid=stat.nvalid(windowBuffer);
504 if(nvalid){
505 std::vector<short>::const_iterator vit=m_class.begin();
506 while(vit!=m_class.end())
507 outBuffer[x/down]+=100.0*occurrence[*(vit++)]/nvalid;
508 }
509 else
510 outBuffer[x/down]=noDataValue;
511 break;
512 }
513 case(filter2d::countid):{
514 if(occurrence.size())
515 outBuffer[x/down]=occurrence.size();
516 else
517 outBuffer[x/down]=noDataValue;
518 break;
519 }
520 case(filter2d::mode):{
521 if(occurrence.size()){
522 std::map<int,int>::const_iterator maxit=occurrence.begin();
523 for(std::map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
524 if(mit->second>maxit->second)
525 maxit=mit;
526 }
527 if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)
528 outBuffer[x/down]=maxit->first;
529 else
530 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
531 }
532 else
533 outBuffer[x/down]=noDataValue;
534 break;
535 }
536 case(filter2d::threshold):{
537 assert(m_class.size()==m_threshold.size());
538 int nvalid=stat.nvalid(windowBuffer);
539 if(nvalid>0){
540 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
541 for(int iclass=0;iclass<m_class.size();++iclass){
542 if(100.0*(occurrence[m_class[iclass]])/nvalid>m_threshold[iclass])
543 outBuffer[x/down]=m_class[iclass];
544 }
545 }
546 else
547 outBuffer[x/down]=noDataValue;
548 break;
549 }
550 case(filter2d::scramble):{
551 if(windowBuffer.size()){
552 int randomIndex=std::rand()%windowBuffer.size();
553 if(randomIndex>=windowBuffer.size())
554 outBuffer[x/down]=windowBuffer.back();
555 else if(randomIndex<0)
556 outBuffer[x/down]=windowBuffer[0];
557 else
558 outBuffer[x/down]=windowBuffer[randomIndex];
559 }
560 else
561 outBuffer[x/down]=noDataValue;
562 break;
563 }
564 case(filter2d::mixed):{
565 enum MixType { BF=11, CF=12, MF=13, NF=20, W=30 };
566 double nBF=occurrence[BF];
567 double nCF=occurrence[CF];
568 double nMF=occurrence[MF];
569 double nNF=occurrence[NF];
570 double nW=occurrence[W];
571 if(windowBuffer.size()){
572 if((nBF+nCF+nMF)&&(nBF+nCF+nMF>=nNF+nW)){
573 if(nBF/(nBF+nCF)>=0.75)
574 outBuffer[x/down]=BF;
575 else if(nCF/(nBF+nCF)>=0.75)
576 outBuffer[x/down]=CF;
577 else
578 outBuffer[x/down]=MF;
579 }
580 else{
581 if(nW&&(nW>=nNF))
582 outBuffer[x/down]=W;
583 else
584 outBuffer[x/down]=NF;
585 }
586 }
587 else
588 outBuffer[x/down]=inBuffer[indexJ][indexI];
589 break;
590 }
591 default:
592 break;
593 }
594 }
595 progress=(1.0+y/down);
596 progress/=outputVector.size();
597 pfnProgress(progress,pszMessage,pProgressArg);
598
599 outputVector[y/down]=outBuffer;
600 }
601}