TestMedian.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <Base/Common/CompareFloatingPoint.hh>
#include <Base/Debug/Error.hh>
#include <Base/Math/Random.hh>
#include <MathAlgo/Median1D.hh>
#include <iostream>
#include <vector>

using namespace BIAS;
using namespace std;

/**
   @file TestMedian.cpp
   @brief Test median computation for vectors.
   @relates Median1D
   @ingroup g_tests
   @author MIP
*/

int main(int argc, char *argv[])
{
  Median1D<double> medDouble;
  Median1D<float>  medFloat;
  Median1D<int>    medInt;

  vector<double>   vecDouble;
  vector<float>    vecFloat;
  vector<int>      vecInt;

  const unsigned int numRuns =  1000;
  const double       minVal  = -100.0;
  const double       maxVal  =  100.0;
  const unsigned int minSize =  3;
  const unsigned int maxSize =  150;

  Random randomizer;
  unsigned int size, medIdx, countLT, countGT, numLT, numGT;
  double median;
  bool addLT;
  int result = 0;
  
  for (unsigned int run = 1; run <= numRuns; run++)
  {
    // create random vectors
    size   = randomizer.GetUniformDistributedInt(minSize, maxSize);
    median = randomizer.GetUniformDistributed(minVal, maxVal);
    medIdx = randomizer.GetUniformDistributedInt(0, (int)size-1);
    vecDouble.resize(size);
    vecFloat.resize(size);
    vecInt.resize(size);
    vecDouble[medIdx] = median;
    countLT = countGT = 0;
    numGT = size / 2;
    numLT = size - 1 - numGT;
    for (unsigned int i = 0; i < size; i++) {
      if (i != medIdx) {
        if (countLT == numLT)
          addLT = false;
        else if (countGT == numGT)
          addLT = true;
        else
          addLT = randomizer.GetUniformDistributedInt(0, 100) < 50;
        if (addLT) {
          vecDouble[i] = randomizer.GetUniformDistributed(minVal, median);
          countLT++;
        } else {
          vecDouble[i] = randomizer.GetUniformDistributed(median, maxVal);
          countGT++;
        }
      }
      vecFloat[i] = (float)vecDouble[i];
      vecInt[i] = (int)vecDouble[i];
    }
    BIASASSERT(countLT + countGT + 1 == size);
    BIASASSERT(vecDouble.size() == size && vecFloat.size() == size && vecInt.size() == size);
    if (countLT != numLT || countGT != numGT) {
      cout << "Error in run " << run << "/" << numRuns << " : " << endl
           << "  found " << countLT << " values below median, expected "
           << numLT << endl << "  found " << countGT << " values above median, expected "
           << numGT << endl << ", vector size is " << size << endl << flush;
      result = -1;
    }

    // compute and compare median
    medDouble.Compute(vecDouble);
    medFloat.Compute(vecFloat);
    medInt.Compute(vecInt);
    if (!Equal(medDouble.GetMedian(), median, 1e-4) ||
        !Equal((double)medFloat.GetMedian(), median, 1e-4) ||
        medInt.GetMedian() != (int)median) {
      cout << "Error in run " << run << "/" << numRuns << " : " << endl
           << "  [double] computed " << medDouble.GetMedian()
           << ",\t expected " << median << endl
           << "  [float]  computed " << medFloat.GetMedian()
           << ",\t expected " << (float)median << endl
           << "  [int]    computed " << medInt.GetMedian()
           << ",\t expected " << (int)median << endl
           << "  vector size is " << size << endl << flush;
      result = -1;
    }
  }
  
  return result;
}