NICE_Segmentation/RSMarkovCluster.cpp

#include "RSMarkovCluster.h"

/*
 Alle Koordinatenangaben orientieren sich an der Matrixschreibweise, dh. (y,x)
 Zeile y und Spalte x. Gespeichert werden Koordianten (soweit nicht anders angegeben)
 im typedef Coord. Dieses baut auf std::pair<int,int> auf. Damit ist coord.first die
 Zeile und coord.second die Spalte.
*/


using namespace OBJREC;
using namespace NICE;
using namespace std;

//#define DEBUGMODE
#ifdef DEBUGMODE
#include <boost/date_time/posix_time/posix_time.hpp>
using namespace boost::posix_time;
#endif

///////////////////// ///////////////////// /////////////////////
//                   CONSTRUCTORS / DESTRUCTORS
///////////////////// ///////////////////// /////////////////////

RSMarkovCluster::RSMarkovCluster () : RegionSegmentationMethod()
{
  this->r = 4.5;
  this->mu = 10.0;
  this->p = 1.4;
  this->chaosThreshold = 0.001;
  this->iterations = 23;
}

RSMarkovCluster::RSMarkovCluster ( const NICE::Config * _conf )
{
  this->initFromConfig ( _conf );
}

RSMarkovCluster::~RSMarkovCluster()
{
}

void RSMarkovCluster::initFromConfig( const NICE::Config* _conf, const std::string & _confSection )
{
  //NOTE previously, the confSection defaulted to "MarkovCluster" !!!!!

  this->r               = _conf->gD ( _confSection, "clusterradius",        4.5 );
  this->mu              = _conf->gD ( _confSection, "edgeweight_parameter", 10.0 );
  this->p               = _conf->gD ( _confSection, "inflation_parameter",  1.4 );
  this->chaosThreshold  = _conf->gD ( _confSection, "chaos_threshold",      0.001 );
  this->iterations      = _conf->gI ( _confSection, "iterations",           23 );
}

///////////////////// ///////////////////// /////////////////////
//                      SEGMENTATION STUFF
///////////////////// ///////////////////// //////////////////


/*
 offsets = {(y,x) | norm((y,x),2)<=r}
 Menge aller moeglichen Offsets, dh. offsets enthaelt alle Koordinatenverschiebungen,
 welche vom Ursprung (0,0) gesehen kleiner als der gegebene compact-pruning-Parameter
 ist. Dadurch wird die Nachbarschaft definiert und die Groesse der Cluster begrenzt.
*/
void
RSMarkovCluster::precalcOffsets ( RSMarkovCluster::Offsets& offsets ) const
{
  // Fuege in die Offsets alle Pixel der 8er Nachbarschaft und sich selber ein.
  offsets.push_back ( Coord ( 0, 0 ) );
  offsets.push_back ( Coord ( -1, -1 ) );
  offsets.push_back ( Coord ( -1, 0 ) );
  offsets.push_back ( Coord ( -1, 1 ) );
  offsets.push_back ( Coord ( 0, -1 ) );
  offsets.push_back ( Coord ( 0, 1 ) );
  offsets.push_back ( Coord ( 1, -1 ) );
  offsets.push_back ( Coord ( 1, 0 ) );
  offsets.push_back ( Coord ( 1, 1 ) );


  int bound = round ( this->r );


  for ( int m = -bound; m <= bound; m++ )
  {
    for ( int n = -bound; n <= bound; n++ )
    {
      Coord newCoord ( m, n );

      if ( ( newCoord.pnorm ( 2 ) <= ( this->r ) )
           && ! ( ( m == -1 ) && ( n == -1 ) )
           && ! ( ( m == -1 ) && ( n == 0 ) )
           && ! ( ( m == -1 ) && ( n == 1 ) )
           && ! ( ( m == 0 ) && ( n == -1 ) )
           && ! ( ( m == 0 ) && ( n == 0 ) )
           && ! ( ( m == 0 ) && ( n == 1 ) )
           && ! ( ( m == 1 ) && ( n == -1 ) )
           && ! ( ( m == 1 ) && ( n == 0 ) )
           && ! ( ( m == 1 ) && ( n == 1 ) ) )

      {
        offsets.push_back ( newCoord );
      }
    }
  }
}

/*
 detours[e] = {(i,j) | offsets(i)+offsets(j)=offsets(e)}
 Menger aller moeglichen 2er Pfade von einem gegebenen Pixel zum offset[e].
*/
void
RSMarkovCluster::precalcDetours ( RSMarkovCluster::Detours& detours, const RSMarkovCluster::Offsets& offsets ) const
{
  // Anzahl der offsets (Nachbarn)
  const uint N_e = offsets.size();
  // Groesse von Detours anpassen
  detours.resize ( N_e );

  // diese Tour muss fuer alle Nachbarn berechnet werden
  for ( uint e = 0; e < N_e; e++ )
  {
    // Touren fuer den Offset e
    vector< pair< uint, uint > > detours_e;

    for ( uint i = 0; i < N_e; i++ )
    {
      for ( uint j = 0; j < N_e; j++ )
      {
        if ( ( offsets[i].first + offsets[j].first == offsets[e].first ) && ( offsets[i].second + offsets[j].second == offsets[e].second ) )
        {
          detours_e.push_back ( pair< uint, uint> ( i, j ) );
        }
      }
    }
    detours[e] = detours_e;
  }
}

int
RSMarkovCluster::segRegions ( const Image& img, Matrix& mask ) const
{
  RSMarkovCluster::Offsets offsets;
  RSMarkovCluster::Detours detours;

  // Vorausberechnungen der fuer das Clustering erforderlichen Offsets und Detours.
  precalcOffsets ( offsets );
  const uint n_e = offsets.size();

  precalcDetours ( detours, offsets );

  // Bilddaten
  const uchar* imageData    = img.getPixelPointerXY ( 0, 0 );
  const uint   imageHeight  = img.height();
  const uint   imageWidth   = img.width();
  const uint   channelCount = img.channels();

  // MarkovMatrix
  NodeCentricRepMatrix L ( imageHeight, imageWidth, n_e );

  // Initialisierung der MarkovMatrix.
  initMarkovMatrix ( imageData, imageHeight, imageWidth, channelCount, offsets, L );


#ifdef DEBUGMODE
  // Anfangszeitpunkt
  ptime start ( microsec_clock::local_time() ); //
#endif
  // Start des MarkovClustering Prozesses
  runClustering ( offsets, detours, L );
#ifdef DEBUGMODE
  // Endzeitpunkt
  ptime end ( microsec_clock::local_time() );   //
  printf ( "[log] Laufzeit Clustering: %ld.%3ldsek\n", ( end - start ).total_milliseconds() / 1000, ( end - start ).total_milliseconds() % 1000 );
#endif

  // Finde die Cluster innerhalb der Matrix
  const uint clusterCount = findCluster ( offsets, L, mask );

  return clusterCount;
}

int
RSMarkovCluster::segRegions ( const ColorImage &cimg, Matrix &mask ) const
{
  RSMarkovCluster::Offsets offsets;
  RSMarkovCluster::Detours detours;

  // Vorausberechnungen der fuer das Clustering erforderlichen Offsets und Detours.
  precalcOffsets ( offsets );
  const uint n_e = offsets.size();

  precalcDetours ( detours, offsets );

  // Bilddaten
  const uchar* imageData    = cimg.getPixelPointerXY ( 0, 0 );
  const uint   imageHeight  = cimg.height();
  const uint   imageWidth   = cimg.width();
  const uint   channelCount = cimg.channels();

  // MarkovMatrix
  NodeCentricRepMatrix L ( imageHeight, imageWidth, n_e );

  // Initialisierung der MarkovMatrix.
  initMarkovMatrix ( imageData, imageHeight, imageWidth, channelCount, offsets, L );

#ifdef DEBUGMODE
  // Anfangszeitpunkt
  ptime start ( microsec_clock::local_time() ); //
#endif
  // Start des MarkovClustering Prozesses
  runClustering ( offsets, detours, L );
#ifdef DEBUGMODE
  // Endzeitpunkt
  ptime end ( microsec_clock::local_time() );   //
  printf ( "[log] Laufzeit Clustering: %ld.%3ldsek\n", ( end - start ).total_milliseconds() / 1000, ( end - start ).total_milliseconds() % 1000 );
#endif

  // Finde die Cluster innerhalb der Matrix
  const uint clusterCount = findCluster ( offsets, L, mask );

  return clusterCount;
}

int
RSMarkovCluster::findCluster ( const RSMarkovCluster::Offsets & offsets, NodeCentricRepMatrix& L, Matrix& mask ) const
{
  // Bild zum Anzeigen der Attracktoren
  Matrix att ( L.getM(), L.getN(), 0.0 );
  // L' Matrix: dort werden die Veraenderungen der L-Matrix gespeichert.
  NodeCentricRepMatrix M ( L.getM(), L.getN(), L.getE() );


  Image attImage ( L.getN(), L.getM() );
  for ( uint y = 0;y < L.getM();y++ )
  {
    for ( uint x = 0;x < L.getN();x++ )
    {
      attImage.setPixel ( x, y, 0 );
    }
  }


  printf ( "schwache Kanten entfernen\n" );

  // Delta um den Kantenschwellwert zu berechnen. Dieser wird von der Summe der abgehenden Kanten abgezogen
  const double delta = 0.01;
  // Loeschen von schwachten Kanten & Finden der Attracktoren.
  for ( uint y = 0;y < L.getM();y++ )
  {
    for ( uint x = 0;x < L.getN();x++ )
    {
      // Pointer auf das die "Nachbarschaft"
      double* efirst = &L ( y, x, 0 );
      // Alle Kanten, die vom Knoten (y,x) abgehen
      valarray<double> outgoingEdges = valarray<double> ( efirst, L.getE() );
      // die Summe der quadratischen abgehenden Kantengewichte
      double edgeThreshold = pow ( outgoingEdges, 2.0 ).sum() - delta;

      for ( uint e = 0;e < L.getE();e++ )
      {
        if ( L ( y, x, e ) >= edgeThreshold )
        {
          M ( y, x, e ) = L ( y, x, e );
        }
        else
        {
          M ( y, x, e ) = 0.0;
        }
        if ( L ( y, x, e ) > L ( y, x, 0 ) )
        {
          M ( y, x, e ) = L ( y, x, e );
        }
        else
        {
          M ( y, x, e ) = 0.0;
        }

      }

      // Pointer auf das die "Nachbarschaft"
      efirst = &M ( y, x, 0 );
      // Alle Kanten, die vom Knoten (y,x) abgehen
      outgoingEdges = valarray<double> ( efirst, M.getE() );
      // die Summe der quadratischen abgehenden Kantengewichte
      double edgeSum = outgoingEdges.sum();

      if ( edgeSum == 0.0 )
        //if(M(y,x,0)!=0.0)
      {
        att ( y, x ) = 1.0;
        if ( L ( y, x, 0 ) > 0.75 )
          attImage.setPixel ( x, y, 255 );
        else if ( L ( y, x, 0 ) > 0.5 )
          attImage.setPixel ( x, y, 180 );
        else if ( L ( y, x, 0 ) > 0.25 )
          attImage.setPixel ( x, y, 100 );
        else
          attImage.setPixel ( x, y, 50 );
      }
    }
  }

  showImage ( attImage );

  printf ( "aquivalenzklassen finden\n" );
  for ( uint y = 0;y < att.rows();y++ )
  {
    for ( uint x = 0;x < att.cols();x++ )
    {
      // aktueller Pixel ist ein Attraktor
      if ( att ( y, x ) == 1.0 )
      {
        uint strongestAttIndex = 0;
        double strongestAttValue = 0.0;

        for ( uint e = 1;e < L.getE();e++ )
        {
          // "-->" Relation
          if ( L ( y, x, e ) > 0 )
          {
            Coord st = Coord ( y, x ) + offsets[e];
            if ( att ( st.first, st.second ) == 1.0 )
            {
              // benachbarter Attraktor gefunden

              // Welcher der beiden ist der "staerkere" ?
              if ( ( L ( y, x, 0 ) <= L ( st.first, st.second, 0 ) ) && ( L ( st.first, st.second, 0 ) > strongestAttValue ) )
              {
                strongestAttValue = L ( st.first, st.second, 0 );
                strongestAttIndex = e;
              }
            }
          }
        }

        if ( strongestAttIndex != 0 )
        {
          M ( y, x, strongestAttIndex ) = strongestAttValue;
          att ( y, x ) = 0.0;
          attImage.setPixel ( x, y, 0 );
        }
      }
    }
  }
  showImage ( attImage );

  //return 0;
  // Erste Operation auf der Ausgabemaske: Markieren der Attraktoren

  uint clusterCount = 0;
  printf ( "Attraktoren markieren\n" );
//#pragma omp parallel for
  for ( uint y = 0;y < att.rows();y++ )
  {
    for ( uint x = 0;x < att.cols();x++ )
    {
      if ( att ( y, x ) == 1.0 )
      {
        // Jeder Attraktor bekommt eine eindeutige Identifikation & wird in der Maske markiert
        mask ( y, x ) = ++clusterCount;
      }
    }
  }
  printf ( "Attraktoren markiert\n" );

  printf ( "bla\n" );
  printf ( "finde cluster\n" );

  for ( uint y = 0;y < mask.rows();y++ )
  {
    for ( uint x = 0;x < mask.cols();x++ )
    {
      if ( mask ( y, x ) == 0.0 )
      {
        // (x,y) Pixel wurde noch keinem Cluster zugewiesen
        Coord st = Coord ( y, x );
        do
        {
          // Pointer auf das die "Nachbarschaft"
          const double* efirst = &M.at ( st.first, st.second, 0 );
          // Alle Kanten, die vom Knoten (y,x) abgehen
          const valarray<double> outgoingEdges = valarray<double> ( efirst, M.getE() );
          // maximaler Kantenwert, der von (y,x) abgeht
          const double maxEdgeWeight = outgoingEdges.max();

          // (s,t) Nachbarpixel zu dem die staerkste Kante fuerhte
          uint e = 0;
          for ( ;e < M.getE();e++ ) {
            if ( M ( st.first, st.second, e ) == maxEdgeWeight ) break;
          };
          st += offsets[e];
          //printf("(%d,%d)\n",st.second,st.first);
        } while ( att ( st.first, st.second ) == 0.0 );

        mask ( y, x ) = mask ( st.first, st.second );
      }
    }
  }



  //showImage(att);

  return clusterCount;


  //printf("[log] start clustering\n");
  // erzeuge den DAG fuer die L-Matrix. Dh. es werden nur Kanten behalten, welche folgendes Kriterium erfuellen:
  // seien A,B Knoten aus V und sei w(*) das Gewicht der Kante der beiden Knoten. Dann enthaelt der DAG nur die
  // Kanten, fuer die gilt: w(A->B)>w(A->A). Dabei entstehen senken, welche welche von den Attracktoren gebildet
  // werden, da diese diese Bedingung nicht erfüllen. Dieses Vorgehen erzeugt ein ueberlappendes Clustering, dass
  // in der Bildsegmentierung nicht gewuenscht ist. Darum wird der DAG wie folgt erweitert:
  // Def. siehe oben. Zusaetzliche Bedingung ist, dass w(A->B) die staerkste abgehende Kante von A ist.

  // Maske zum Speichern der max. Kante. Dabei wird der entsprechende Offsetindex gespeichert.
  //Matrix maxEdgeIndexMatrix(L.getM(),L.getN(),-1);

  /*
  #pragma omp parallel for
  for(uint y=0;y<L.getM();y++)
  {
   for(uint x=0;x<L.getN();x++)
   {
    // Pointer auf das die "Nachbarschaft"
    double* efirst=&L.at(y,x,0);
    // Alle Kanten, die vom Knoten (y,x) abgehen
    valarray<double> outgoingEdges=valarray<double>(efirst,L.getE());
    // maximaler Kantenwert, der von (y,x) abgeht
    double maxEdgeWeight=outgoingEdges.max();


    if(maxEdgeWeight==0.5) printf("[debug] max. edgeweight=%.2f\n",maxEdgeWeight);

    // Alle Kanten loeschen, die kleineres Gewicht als
    for(uint e=0;e<L.getE();e++)
    {
     if(L.at(y,x,e)!=maxEdgeWeight)
     {
      L.at(y,x,e)=0.0;
     }
     else
     {
      // Index der maximalen Kante wird gespeichert.
      maxEdgeIndexMatrix(y,x)=e;
     }

    }
   }
  }
  */

  /*
  for(uint y=0;y<L.getM();y++)
  {
   for(uint x=0;x<L.getN();x++)
   {
    // wenn es sich um einen Attracktor handelt, dann ist der Eintrag in der maxEdgeIndexMatrix=0
    if(maxEdgeIndexMatrix(y,x)==0)
    {
     clusterCount++;
     mask(y,x)=indexOfPixel(Coord(y,x),L.getN());
     continue;
    }

    Coord st(y,x);
    do
    {
     // Update st auf die Koordinaten des Knotens zu dem die staerkste Kante geht.
     st+=offsets[maxEdgeIndexMatrix(st.first,st.second)];

     // Wir sind auf einen Pixel gestossen, der bereits markiert ist.
     if(mask(st.first,st.second)!=(-1))
     {
      mask(y,x)=mask(st.first,st.second);
      break;
     }
    }while(maxEdgeIndexMatrix(st.first,st.second)!=0);

    if(mask(y,x)==(-1)) mask(y,x)=indexOfPixel(st,L.getN());
   }
  }*/


  //return clusterCount;
}

inline uint
RSMarkovCluster::indexOfPixel ( const Coord& yx, const uint xSize ) const
{
  return yx.first*xSize + yx.second;
}


void
RSMarkovCluster::normalize ( NodeCentricRepMatrix& L ) const
{
#pragma omp parallel for
  for ( uint y = 0;y < L.getM();y++ )
  {
    for ( uint x = 0;x < L.getN();x++ )
    {
      // Pointer auf das die "Nachbarschaft"
      double* efirst = &L ( y, x, 0 );
      // Alle Kanten, die vom Knoten (y,x) abgehen
      valarray<double> outgoingEdges = valarray<double> ( efirst, L.getE() );
      // maximaler Kantenwert, der von (y,x) abgeht
      double sum = outgoingEdges.sum();

      for ( uint e = 0;e < L.getE();e++ )
      {
        L ( y, x, e ) = checkForNAN ( L ( y, x, e ) / sum );
      }

    }
  }

}

void
RSMarkovCluster::inflation ( const NodeCentricRepMatrix& Lin, NodeCentricRepMatrix& Lout ) const
{
#pragma omp parallel for
  for ( uint y = 0;y < Lin.getM();y++ )
  {
    for ( uint x = 0;x < Lin.getN();x++ )
    {
      for ( uint e = 0;e < Lin.getE();e++ )
      {
        Lout ( y, x, e ) = checkForNAN ( pow ( Lin ( y, x, e ), this->p ) );
      }
    }
  }
}

void
RSMarkovCluster::expansion ( const NodeCentricRepMatrix& Lin, NodeCentricRepMatrix& Lout, const Offsets& offsets, const Detours& detours ) const
{
#pragma omp parallel for
  for ( uint y = 0;y < Lin.getM();y++ )
  {
    for ( uint x = 0;x < Lin.getN();x++ )
    {
      // neues Kantengewicht fuer alle "adjazenten" Knoten berechnen
      for ( uint e = 0;e < Lin.getE();e++ )
      {
        double w_xy_mn = 0.0;

        for ( uint d = 0;d < detours[e].size();d++ )
        {
          uint efirst = detours[e][d].first;
          uint esecond = detours[e][d].second;

          Coord st = Coord ( y, x ) + offsets[efirst];

          double w_xy_st = 0.0;
          double w_st_mn = 0.0;

          if ( ( st.first >= 0 ) && ( st.first < (int)Lin.getM() ) && ( st.second >= 0 ) && ( st.second < (int)Lin.getN() ) )
          {
            w_xy_st = Lin ( y, x, efirst );
            w_st_mn = Lin ( st.first, st.second, esecond );
          }

          w_xy_mn += w_xy_st * w_st_mn;
        }
        Lout ( y, x, e ) = checkForNAN ( w_xy_mn );
      }
    }
  }
}

inline double
RSMarkovCluster::checkForNAN ( const double val ) const
{
  if ( ! ( val == val ) || ( val == -std::numeric_limits < double >::quiet_NaN () ) )
    return 0.0;
  else
    return val;
}

void
RSMarkovCluster::runClustering ( const RSMarkovCluster::Offsets& offsets, const RSMarkovCluster::Detours& detours, NodeCentricRepMatrix& L1 ) const
{
  //printMatrix(L1,offsets,"Linit");

  // Initialisierte Markovmatrix muss normalisiert werden.
  normalize ( L1 );

  // Mass fuer die Aenderungen innerhalb der Matrix.
  //double chaos = 1.0;
  //uint sameChaosCounter = 0; // Wenn sich Chaos eine bestimmte Anzahl von Iterationen nicht mehr veraendert, dann wird abgebrochen.
  // 2te Markovmatrix zum Vergleichen der Veraenderungen
  NodeCentricRepMatrix L2 ( L1.getM(), L1.getN(), L1.getE() );

  uint iterationCounter = 0;

#ifdef DEBUGMODE
  time_duration expAvgTime ( 0, 0, 0, 0 );
  time_duration infAvgTime ( 0, 0, 0, 0 );
  time_duration norAvgTime ( 0, 0, 0, 0 );
#endif
  printf ( "%d\n", iterations );

  //while(chaos>chaosThreshold)
  while ( iterationCounter < (uint)iterations )
  {
#ifdef DEBUGMODE
    // Anfangszeitpunkt
    ptime start ( microsec_clock::local_time() ); //
#endif
    // rufe Matrix L1 und speichere in Matrix L2 (L2<--exp(L1))
    expansion ( L1, L2, offsets, detours );
#ifdef DEBUGMODE
    // Endzeitpunkt
    ptime end ( microsec_clock::local_time() );   //

    // Zeit, die der Expansionschritt gebraucht hat.
    time_duration expDuration ( end - start );
    printf ( "[log] expansion: %ld.%3ldsek\n", expDuration.total_milliseconds() / 1000, expDuration.total_milliseconds() % 1000 );

    // Gesamtlaufzeit des Expansionschritts
    expAvgTime += expDuration;
#endif
#ifdef DEBUGMODE
    // Anfangszeitpunkt
    start = microsec_clock::local_time(); //
#endif
    // rufe Matrix L2 und speichere in Matrix L1 (L1<--inf(L2))
    inflation ( L2, L1 );
#ifdef DEBUGMODE
    // Endzeitpunkt
    end = microsec_clock::local_time(); //

    // Zeit, die der Expansionschritt gebraucht hat.
    time_duration infDuration ( end - start );
    printf ( "[log] inflation: %ld.%3ldsek\n", infDuration.total_milliseconds() / 1000, infDuration.total_milliseconds() % 1000 );

    // Gesamtlaufzeit des Expansionschritts
    infAvgTime += infDuration;
#endif
#ifdef DEBUGMODE
    // Anfangszeitpunkt
    start = microsec_clock::local_time(); //
#endif
    // vor dem Ende jeder Iteration muss normalisiert werden.
    normalize ( L1 );
#ifdef DEBUGMODE
    // Endzeitpunkt
    end = microsec_clock::local_time(); //

    // Zeit, die der Expansionschritt gebraucht hat.
    time_duration norDuration ( end - start );
    printf ( "[log] normalize: %ld.%3ldsek\n", norDuration.total_milliseconds() / 1000, norDuration.total_milliseconds() % 1000 );

    // Gesamtlaufzeit des Expansionschritts
    norAvgTime += norDuration;
#endif

    //double lastChaos=chaos;
    // berechne die Veraenderungen der Matrix nach dem Schritt t
    //double maxDiffValue=NodeCentricRepMatrix::maxValue(L1-L2);
    //chaos=min(maxDiffValue,chaos);

    //printf("[log] %d iteration: chaos=%.5f\n",iterationCounter,chaos);
    //printf("[debug] %.5f\n",maxDiffValue);

    iterationCounter++;
    printf ( "[log] %d. Iteration\n", iterationCounter );

    /*
    if(lastChaos==chaos)
    {
     sameChaosCounter++;

     if(sameChaosCounter==maxConstantInteration)
     {
      break;
     }
    }
    else
    {
     sameChaosCounter=0;
    }
    */
  }

#ifdef DEBUGMODE
  long int expAvgTimeUSec = expAvgTime.total_milliseconds() / ( iterationCounter - 1 );
  long int infAvgTimeUSec = infAvgTime.total_milliseconds() / ( iterationCounter - 1 );
  long int norAvgTimeUSec = norAvgTime.total_milliseconds() / ( iterationCounter - 1 );

  printf ( "[log] Durchschnittszeit fuer Expansion: %ld.%3ldsek\n", expAvgTimeUSec / 1000, expAvgTimeUSec % 1000 );
  printf ( "[log] Durchschnittszeit fuer Inflation: %ld.%3ldsek\n", infAvgTimeUSec / 1000, infAvgTimeUSec % 1000 );
  printf ( "[log] Durchschnittszeit fuer Normalize: %ld.%3ldsek\n", norAvgTimeUSec / 1000, norAvgTimeUSec % 1000 );
#endif
}


void
RSMarkovCluster::initMarkovMatrix ( const uchar* imageData, const uint imageHeight, const uint imageWidth, const uint channelCount, const RSMarkovCluster::Offsets& offsets, NodeCentricRepMatrix& L ) const
{
  // Initialisierung erfolgt fuer jeden Pixel und die 8er Nachbarschaft. Dabei wird im Loops initialisiert.
#pragma omp parallel for
  for ( uint y = 0;y < imageHeight;y++ )
  {
    for ( uint x = 0;x < imageWidth;x++ )
    {
      for ( uint e = 0; e < 9; e++ )
      {
        // Koordinaten des Nachbarpixels
        RSMarkovCluster::Coord yxN = RSMarkovCluster::Coord ( y, x ) + offsets[e];

        // Eine Verbindung zu einem Nachbarn ausserhalb des Bildes hat das Gewicht 0.
        if ( ( yxN.first < 0 ) || ( yxN.first >= (int)imageHeight ) || ( yxN.second < 0 ) || ( yxN.second >= (int)imageWidth ) )
          continue;

        Vector imageValue ( channelCount );

        for ( uint c = 0;c < channelCount;c++ )
        {
          uint indexP = y * imageWidth * channelCount + x * channelCount + c;
          uint indexN = yxN.first * imageWidth * channelCount + yxN.second * channelCount + c;

          imageValue[c] = ( double ( imageData[indexP] ) - double ( imageData[indexN] ) ) / 255.0;
        }

        L.at ( y, x, e ) = edgeWeight ( imageValue );
      }
    }
  }
}

inline double
RSMarkovCluster::edgeWeight ( const Vector& imageValueDiff ) const
{
  return exp ( - ( this->mu ) * pow ( imageValueDiff.normL2 (), 2 ) );
}

void
RSMarkovCluster::printMatrix ( const NodeCentricRepMatrix& L, const RSMarkovCluster::Offsets& offsets , const string& filename ) const
{
  // Bemerkung in Offset ist [0]=x && [1]=y
  const uint n_x = L.getN ();
  const uint n_y = L.getM ();
  const uint n_e = L.getE ();
  Matrix oMatrix ( ( n_y*n_x ), ( n_y*n_x ) );

  Matrix::iterator it = oMatrix.begin();
  for ( ;it != oMatrix.end();++it ) {
    ( *it ) = 0;
  }

  for ( uint y = 0;y < n_y;y++ )
  {
    for ( uint x = 0;x < n_x;x++ )
    {
      for ( uint e = 0;e < n_e;e++ )
      {
        double value = L.at ( y, x, e );
        int n = offsets[e].second;
        int m = offsets[e].first;

        int indexP = y * n_x + x;
        int indexN = ( y + m ) * n_x + ( x + n );

        if ( ( indexP >= 0 ) && ( indexP < (int)( n_y*n_x ) ) && ( indexN >= 0 ) && ( indexN < (int)( n_y*n_x ) ) )
          oMatrix ( indexP, indexN ) = value;
      }
    }
  }

  ofstream fout;
  fout.open ( filename.c_str() );
  for ( int y = 0;y < (int)( n_y*n_x );y++ )
  {
    string line = "";
    for ( int x = 0;x < (int)( n_y*n_x );x++ )
    {
      if ( oMatrix ( x, y ) == 0.0 ) line += "----\t";
      else
      {
        char entry[16];
        sprintf ( entry, "%.3f\t", oMatrix ( x, y ) );
        line += string ( entry );
      }

    }
    fout << line << endl;
  }
  fout.close();
}

///////////////////// INTERFACE PERSISTENT /////////////////////
// interface specific methods for store and restore
///////////////////// INTERFACE PERSISTENT ///////////////////// 

void RSMarkovCluster::restore ( std::istream & is, int format )
{
  //delete everything we knew so far...
  this->clear();
  
  
  if ( is.good() )
  {
    
    std::string tmp;
    is >> tmp; //class name 
    
    if ( ! this->isStartTag( tmp, "RSMarkovCluster" ) )
    {
      std::cerr << " WARNING - attempt to restore RSMarkovCluster, but start flag " << tmp << " does not match! Aborting... " << std::endl;
      throw;
    }   
    
    bool b_endOfBlock ( false ) ;
    
    while ( !b_endOfBlock )
    {
      is >> tmp; // start of block 
      
      if ( this->isEndTag( tmp, "RSMarkovCluster" ) )
      {
        b_endOfBlock = true;
        continue;
      }      
      
      tmp = this->removeStartTag ( tmp );
     
      if ( tmp.compare("r") == 0 )
      {
        is >> this->r;
        is >> tmp; // end of block 
        tmp = this->removeEndTag ( tmp );
      }
      else if ( tmp.compare("mu") == 0 )
      {
        is >> this->mu;
        is >> tmp; // end of block 
        tmp = this->removeEndTag ( tmp );
      }      
      else if ( tmp.compare("p") == 0 )
      {
        is >> this->p;
        is >> tmp; // end of block 
        tmp = this->removeEndTag ( tmp );
      }
      else if ( tmp.compare("chaosThreshold") == 0 )
      {
        is >> this->chaosThreshold;
        is >> tmp; // end of block 
        tmp = this->removeEndTag ( tmp );
      }
      else if ( tmp.compare("iterations") == 0 )
      {
        is >> this->iterations;
        is >> tmp; // end of block 
        tmp = this->removeEndTag ( tmp );
      }  
      else
      {
      std::cerr << "WARNING -- unexpected RSMarkovCluster object -- " << tmp << " -- for restoration... aborting" << std::endl;
      throw;
      }
    }
  }
  else
  {
    std::cerr << "RSMarkovCluster::restore -- InStream not initialized - restoring not possible!" << std::endl;
    throw;
  }
}

void RSMarkovCluster::store ( std::ostream & os, int format ) const
{ 
  if (os.good())
  {
    // show starting point
    os << this->createStartTag( "RSMarkovCluster" ) << std::endl;    

    os << this->createStartTag( "r" ) << std::endl;
    os << this->r << std::endl;
    os << this->createEndTag( "r" ) << std::endl;
    
    os << this->createStartTag( "mu" ) << std::endl;
    os << this->mu << std::endl;
    os << this->createEndTag( "mu" ) << std::endl;
    
    os << this->createStartTag( "p" ) << std::endl;
    os << this->p << std::endl;
    os << this->createEndTag( "p" ) << std::endl;  

    os << this->createStartTag( "chaosThreshold" ) << std::endl;
    os << this->chaosThreshold << std::endl;
    os << this->createEndTag( "chaosThreshold" ) << std::endl;  

    os << this->createStartTag( "iterations" ) << std::endl;
    os << this->iterations << std::endl;
    os << this->createEndTag( "iterations" ) << std::endl;  
    
    // done
    os << this->createEndTag( "RSMarkovCluster" ) << std::endl;    
  }
  else
  {
    std::cerr << "OutStream not initialized - storing not possible!" << std::endl;
  }
}

void RSMarkovCluster::clear ()
{
}