biovision_d.h

Go to the documentation of this file.

/* Copyright (C) 2000-2016 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */
 
#ifndef BRAINRAT_CLASSIFICATION_BIOVISION_D_H
#define BRAINRAT_CLASSIFICATION_BIOVISION_D_H
 
//--- brainrat ----------------------------------------------------------
#include <brainrat/classification/biovision.h>      // class declaration
#include <brainrat/data/classes.h>                      // bio::BVclasses
//--- aims -------------------------------------------------------------------
#include <aims/data/data_g.h>                                      // AimsData
#include <aims/io/writer.h>                                    // aims::Writer
#include <aims/io/reader.h>                                    // aims::Reader
#include <aims/utility/converter_g.h>
#include <aims/utility/progress.h>                      // aims::ProgressInfo
#include <aims/threshold/emthreshold.h>                     // AimsEMThreshold
#include <aims/roi/roiIterator.h>                         // aims::RoiIterator
#include <aims/roi/maskIterator.h>                       // aims::MaskIterator
#include <aims/connectivity/connectivity.h>              // aims::Connectivity
#include <aims/connectivity/component.h>// aims::AimsLabeledConnectedComponent
#include <aims/interpolation/interpolatedvolume.h>
//--- cartodata --------------------------------------------------------------
#include <cartodata/volume/volume.h>
//--- soma-io ----------------------------------------------------------------
#include <soma-io/allocator/allocator.h>
//--- cartobase --------------------------------------------------------------
#include <cartobase/stream/fileutil.h>
#include <cartobase/smart/rcptr.h>
#include <cartobase/type/string_conversion.h>
//--- system -----------------------------------------------------------------
#include <vector>
//----------------------------------------------------------------------------
 
 
namespace biovision {
    
//============================================================================
//   L E A R N I N G
//============================================================================
 
  template <typename T> int32_t
  LearningAlgorithm::processOneImage( carto::VolumeRef<T>      & image,
                                      const std::string        & inputGraph,
                                      const std::string        & fileOutMeasure,
                                      std::vector<carto::Object >                     & classDesc,
                                      std::vector<std::vector<std::vector<double> > > & classVectors )
  {
    // Get image processor instance
    carto::rc_ptr<aims::MaskIterator> mask = aims::getMaskIterator( inputGraph );
 
    // Get interpolator instances
    std::vector< aims::InterpolatedVolume > interpolators;
    interpolators.reserve(_measuresize);
    
    for( int j=0; j<_measuresize; ++j) {
      carto::VolumeRef<T> view = image.view( Point4di( 0, 0, 0, j), Point4di(image.getSizeX(), image.getSizeY(), image.getSizeZ(), 1) );
      aims::InterpolatedVolume interpolated;
      interpolated.setVolume( view.copy(), 1 ); 
      interpolators.push_back( interpolated );
    }
 
    // Read ROI graph
    const carto::rc_ptr<aims::RoiIterator> roiIterator = aims::getRoiIterator( inputGraph );
 
    while( roiIterator->isValid() ) {   // Classes
      int k = 0;
      bool newregion = true;
 
      // Test if ROI already exists
      for( k=0; k<classDesc.size(); ++k )
      {
        // macos 10.4 copmiler needs this carto::GenericObject:: thing
        const std::string & s = classDesc[k]->getProperty("n")->carto::GenericObject::value<std::string>();
        if( s == roiIterator->regionName() ) {
          // Add measure vectors to existing region
          if( carto::verbose )
            std::cout << "Add measure vectors to existing region "
                      << roiIterator->regionName() << std::endl;
          newregion = false;
          break;
        }
      }
 
      if( newregion ) {
        // Create a new region
        if( carto::verbose )
          std::cout << "Create new region " << roiIterator->regionName()
                    << std::endl;
 
        carto::Object c = carto::Object::value(carto::PropertySet());
        c->setProperty("n", roiIterator->regionName() ); // Set region prior
        c->setProperty("p", 0.);
        classDesc.push_back(c);
 
        std::vector<std::vector<double> > cv; // Add vector entry for the class vectors
        classVectors.push_back( cv );
      }
 
      carto::rc_ptr<aims::MaskIterator> roiMaskIterator = roiIterator->maskIterator();
      size_t point_count = 0;
 
      // Process each pixel contained in the mask
      while( roiMaskIterator->isValid() ) {
        // Create measure vector v and set values
        std::vector<double> v(_measuresize);
        Point3df coord_f = roiMaskIterator->valueMillimeters();
 
        for( int32_t j=0; j<_measuresize; ++j )
          v[j] = interpolators[j].atMm( coord_f[0], coord_f[1], coord_f[2] );
 
        classVectors[k].push_back(v); // add to existing measure vector
        ++point_count;
 
        roiMaskIterator->next();
      }
 
      if( carto::verbose )
        std::cout << carto::toString(point_count)
                  << " measure vectors added to "
                  << roiIterator->regionName()
                  << std::endl;
 
      roiIterator->next();
    }
 
    return EXIT_SUCCESS;
  }
 
  int32_t LearningAlgorithm::learnClustering( std::list<double>                               & inputListOfPrior,
                                              std::vector<carto::Object >                     & classDesc,
                                              std::vector<std::vector<std::vector<double> > > & classVectors )
  {
    // Clustering
    if (carto::verbose)
      std::cout << "Clustering measure vectors ..." << std::endl
                << "Clustering parameters [ threshold = " << _threshold
                << ", cluster number maximum = " << _nbclustermax
                << ", sigma minimum = " << _sigmamin
                << ", sigma ratio = " << _sigmaratio << " ]" << std::endl;
 
    if( !inputListOfPrior.empty() && inputListOfPrior.size() != classDesc.size() )
        throw ( "number of priors and classes must match! " ) ;
 
    std::list<double>::iterator itp = inputListOfPrior.begin(); // prior iterator
 
    int nbtotal = 0; // Information only : number of pixels in learning set
    std::vector<carto::Array<carto::Array<carto::Array<int> > > > histo(classDesc.size());
 
    /* Compute Density Parameters :
     * For each class (contained in classDesc) is given a set
     * Q=classVectors.at( p ) of vectors. We want to find the number of
     * gaussians in the mixture and thus clusterize the set Q.
     */
    for( int p = 0; p < classDesc.size(); ++p )
    {
      // We get all measures of given class p
      std::vector<std::vector<double> > & mv = classVectors[p];
      std::vector<double> d; // Clusters density parameters
 
      int nbvectors = mv.size(); // Number of measures in class
      nbtotal += nbvectors;
 
      if( carto::verbose ) {
        std::cout << "Clustering " << carto::toString(nbvectors)
                  << " measure vectors of " << classDesc[p]->getProperty("n")->value<std::string>()
                  << "..." << std::endl;
      }
 
      carto::VolumeRef<int>  C( nbvectors ); // Cluster labeling
      C = 0;
      std::vector<std::vector<double> > centroids
        = Clustering( mv, _threshold, _sigmamin, _sigmaratio,
                      _nbclustermax, _newclusterpolicy, _clustersizemin, C ) ; // Centroids of each cluster
 
      // For each cluster, compute cluster probability and covariance matrix
      int nbclusters = centroids.size(); // Named 'nq' in Biovision article
      d.push_back(nbclusters);
 
      if (carto::verbose)
        std::cout << "Processing density parameters for "
                  << carto::toString(nbclusters)
                  << " clusters..." << std::endl;
 
      // ro probability of each cluster
      carto::VolumeRef<double> ro( nbclusters );
      ro = 0;
      for( int i = 0 ; i < nbvectors; ++i ) {
        if( C(i) >= 0 )
          ro( C(i) )++;
      }
 
      // For each cluster Qi
      for( int k = 0 ; k < nbclusters ; ++k ) {
 
        // Cluster probability
        ro(k) /= nbvectors;
        d.push_back(ro(k));
 
        if (carto::verbose)
          std::cout << "ro(" << carto::toString(k) << ") : "
                    << carto::toString(ro(k)) << std::endl;
 
        // Centroid
        std::vector<double> mu = centroids[k];
 
        for (int i = 0; i < _measuresize; ++i) {
          d.push_back(mu[i]);
        }
 
        // All the measure vectors in cluster
        std::vector<std::vector<double> > M;
        std::vector<std::vector<double> > qk;
        for (int i = 0 ; i < nbvectors; ++i) {
          if (C(i) == k)
            qk.push_back(mv[i]);
        }
 
        // Compute cluster histograms
        carto::Array<carto::Array<int> > mh;
        for (int i = 0; i < _measuresize; ++i) {
          aims::RegularBinnedHistogram<double> h;
 
          // Convert measure vectors to VolumeRef
          carto::VolumeRef<double> hdata(qk.size());
          for(int32_t j = 0; j < hdata.getSizeX(); ++j) {
            hdata(j) = ((std::vector<double>)qk[j])[i];
          }
 
          h.doit(hdata, 0, 255);
 
          // Convert histogram to array
          carto::Array<int> v;
          for (int32_t j = 0; j < h.data().getSizeX(); ++j) {
            v.push_back(h.data()(j));
          }
          mh.push_back(v);
        }
        histo[p].push_back(mh);
 
        int nvcluster = qk.size(); // Number of measure vectors in cluster
        d.push_back(nvcluster);
 
        // Compute covariance matrix :  cov =  1/n ( M*M' )
        for (int i = 0 ; i < nvcluster; ++i) {
          std::vector<double> diff(qk[i].size());
          for(int j = 0;j<qk[i].size(); ++j){
              diff[j] = qk[i][j] - mu[j];
          }
          M.push_back(diff);
        }
 
        carto::VolumeRef<double> cov(_measuresize, _measuresize), I(_measuresize, _measuresize);
        I = 0, cov = 0;
        for (int i = 0; (i< _measuresize); ++i) {
          I(i,i) = 1.0;
          for (int j = 0 ; (j < _measuresize); ++j) {
            for (int k = 0 ; k < nvcluster; ++k) {
              cov (i,j) += M[k][i] * M[k][j] ; // M(i,k) * M'(k,j)
            }
            cov(i,j) /= nvcluster; // Normalize covariance matrix
          }
        }
 
        // Inverse of covariance matrix and determinant
        double detcov = det(cov, I);
        d.push_back(detcov);
        for (int i = 0; i < _measuresize; ++i) {
          for (int j = 0 ; j < _measuresize; ++j) {
            d.push_back(I(i,j));
          }
        }
      }
 
      classDesc[p]->setProperty("d", d);        
      classDesc[p]->setProperty("h", histo[p]);
 
      // Change a priori probability of each class
      if (inputListOfPrior.empty ())
        classDesc[p]->getProperty("p")->value<double>() = 1.0 / ((double)classDesc.size());  // Uniform
      else
        classDesc[p]->getProperty("p")->value<double>() = (*itp) ; // Input proba
 
      itp++;
 
      if (carto::verbose){
        std::cout << "prior : " << classDesc[p]->getProperty("p")->value<double>() << std::endl;
        std::cout << "density : [" << d[0];
        for(int i=1; i<d.size(); ++i)
            std::cout << "," << d[i];
        std::cout << "]" << std::endl;
      }
    }
 
    std::cout << carto::toString(nbtotal) << " pixels clustered." << std::endl;
    
    _learningmodel = carto::Object::value(carto::PropertySet());
    _learningmodel->setProperty( "data", classDesc );
 
    return EXIT_SUCCESS;
  }
 
//============================================================================
//   C L A S S I F I C A T I O N
//============================================================================
 
  template <typename T>
  int32_t ClassificationAlgorithm::classify( const Point3dl      &t,
                                             carto::VolumeRef<T> &image )
  {
    std::vector<double>  w(_measuresize);
    for ( int32_t j = 0; j < _measuresize; ++j) {
      w[j] = image(t[0],t[1],t[2],j);
    }
 
    return Classifier(w, t, false, false, true);
  }
 
  template <typename T> int32_t
  ClassificationAlgorithm::classify( carto::VolumeRef<T>            &image,
                                     const std::vector<int32_t>     &labels,
                                     const std::vector<std::string> &classesorder,
                                     const std::string              &maskFile,
                                     const bool                      keepClasses,
                                     const bool                      keepProba,
                                     const bool                      keepDensity)
  {
    // Read or Generate labels from class indexes
    std::vector<int32_t> l = labels;
    if (l.empty()) {
      // Generate labels using class indexes
      for (int32_t i = 0; i < _nbclasses; ++i)
        l.push_back(i + 1);
    }
    else if (l.size() != _nbclasses) {
      throw std::runtime_error( "Number of labels specified does not match the number of classes in model." );
    }
    
    std::map<std::string, int32_t> orders;
    if (classesorder.empty()) {
      // Generate classes order dictionary using class indexes
      for (int32_t i = 0; i < _nbclasses; ++i)
        orders[_classdesc[i]->getProperty("n")->carto::GenericObject::value<std::string>()] = i;
    }
    else if (classesorder.size() != _nbclasses) {
      throw std::runtime_error( "Number of classes specified in classes order does not match the number of classes in model." );
    }
    else {
      // Generate classes order dictionary using specified classesorder
      int32_t i = 0;
      bool found;
      std::vector<std::string>::const_iterator k, ek = classesorder.end();
      for ( k = classesorder.begin(); k != ek; ++k, ++i) {        
        // Checks that k exists in model
        found = false;
        for (int32_t j = 0; j < _classdesc.size(); ++j) {
          if (_classdesc[j]->getProperty("n")->carto::GenericObject::value<std::string>() == *k) {
            found = true;
            break;
          }
        }
        
        if (found) {
          // Sets order for class
          orders[*k] = i;
        }
        else
          throw std::runtime_error( "Class " + *k + " specified in classes order was not found in model." );
      }
    }
    
    // Image dimensions
    int dx = image.getSizeX(), dy = image.getSizeY(), dz = image.getSizeZ();
 
    // Mask image
    carto::rc_ptr<aims::MaskIterator> mask;
 
    // Read mask
    if (!maskFile.empty()) {
      if (carto::verbose)
        std::cout << "Reading mask image ..." << std::endl;
      mask = aims::getMaskIterator( maskFile );
    }
 
    // Classes label image
    if( keepClasses ) {
      if (carto::verbose)
        std::cout << "Allocating classification image ..." << std::endl;
      _Cimage = carto::VolumeRef<int16_t>(dx, dy, dz);
      _Cimage.copyHeaderFrom(image.header());
      _Cimage.fill(0);
    } 
 
    // Probabilities image
    if( keepProba ) {
      if (carto::verbose)
        std::cout << "Allocating probability image ..." << std::endl;
      _Pimage = carto::VolumeRef<double>( dx, dy, dz, _nbclasses );
      _Pimage.copyHeaderFrom(image.header());
      _Pimage.fill(0.);
    }
    
    // Density image
    if( keepDensity ) {
      if (carto::verbose)
        std::cout << "Allocating density image ..." << std::endl;
      _Dimage = carto::VolumeRef<double>( dx, dy, dz, _nbclasses );
      _Dimage.copyHeaderFrom(image.header());
      _Dimage.fill(0.);
    }
    
    // Classify
    if (carto::verbose) {
      std::cout << "Starting segmentation ..." << std::endl;
      for (int32_t i = 0; i < _nbclasses; ++i)
      {
        const std::string & s = _classdesc[i]->getProperty("n")->carto::GenericObject::value<std::string>();
        std::cout << "Class " << s
                  << " label " << carto::toString(l[i]) << std::endl;
      }
    }
 
    std::vector<double> w(_measuresize);
    int ic;
    Point3dl c;
 
    if (!mask.isNull())
      mask->restart();
 
 
    // Progression tracking tool
    aims::Progression progress(dz - 1);
    if (dz > 1)
      std::cout << "Progress: ";
      
    for( int32_t z = 0; z < dz; ++z, ++progress ) {
      if(dz>1) std::cout << progress << std::flush;
      for( int32_t y = 0; y < dy; ++y )
        for( int32_t x = 0; x < dx; ++x )
          if( mask.isNull() || mask->contains( Point3df(
              image.getVoxelSize()[0] * x, image.getVoxelSize()[1] * y, image.getVoxelSize()[2] * z ) ) )
          {
            c = Point3dl( x, y, z );
            for ( int32_t j = 0; j < _measuresize; j++)
              w[j] = image(x,y,z,j);
            ic = Classifier( w, c, keepProba, keepDensity );
            // Classify only if the voxel has a probability to belong 
            // to a class of the model
            if( (ic > -1) && keepClasses ) {
              // Get the class order for the found class
              int32_t o = orders[_classdesc[ic]->getProperty("n")->carto::GenericObject::value<std::string>()];
              _Cimage(x,y,z) = l[ o ];
            }
          }
    }
 
    return EXIT_SUCCESS;
  }
 
  int ClassificationAlgorithm::writeResults( const std::string &outputFile,
                                             const std::string &outputProbaFile,
                                             const std::string &outputDensityFile )
  {
    int out = EXIT_SUCCESS;
    
    // Write proba file
    if( !outputProbaFile.empty() && !_Pimage.get() )
    {
      std::cout << "Proba image was not saved ! "
                << std::endl
                << "Run classification again with keepProba parameter "
                << "set at true."
                << std::endl;
      out = EXIT_FAILURE;
    }
    else if( !outputProbaFile.empty() )
    {
      if( carto::verbose )
        std::cout << "Writing proba image : " << outputProbaFile << std::endl;
 
      aims::Writer<carto::Volume<double> > outputProbaWriter( outputProbaFile );
      outputProbaWriter << *(_Pimage);
    }
 
    // Write density file
    if( !outputDensityFile.empty() && !_Dimage.get() )
    {
      std::cout << "Density image was not saved ! "
                << std::endl
                << "Run classification again with keepDensity parameter "
                << "set at true."
                << std::endl;
      out = EXIT_FAILURE;
    }
    else if( !outputDensityFile.empty() )
    {
      if( carto::verbose )
        std::cout << "Writing density image : " << outputDensityFile << std::endl;
 
      aims::Writer<carto::Volume<double> > outputDensityWriter( outputDensityFile );
      outputDensityWriter << *(_Dimage);
    }
 
    // Write image file
    if( !outputFile.empty() && !_Cimage.get() )
    {
      std::cout << "Classification image was not saved ! "
                << std::endl
                << "Run classification again with keepClasses parameter "
                << "set at true."
                << std::endl;
      out = EXIT_FAILURE;
    }
    else if( !outputFile.empty() )
    {
      if (carto::verbose)
        std::cout << "Writing output image : " << outputFile << std::endl;
 
      aims::Writer<carto::Volume<int16_t> > outputWriter(outputFile);
      outputWriter << *(_Cimage);
    }
 
    return out;
  }
 
}// namespace biovision
 
#endif // BRAINRAT_CLASSIFICATION_BIOVISION_D_H