utility_d.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 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_UTILITY_UTILITY_D_H
#define BRAINRAT_UTILITY_UTILITY_D_H

//--- brainrat ----------------------------------------------------------
#include <brainrat/utility/utility.h>       // class declaration
#include <brainrat/utility/imageprocessor.h>       // bio::ImageProcessor
#include <brainrat/data/measurevector.h>            // bio::MeasureVector
#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/channel.h>                           // ChannelSelector
#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/resampling/linearInterpolator.h>          // aims::Inteprolator
//--- 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 obso_temp {

//============================================================================
//   U T I L I T I E S
//============================================================================

  // Return the brainmask from input image
  template <typename T>
  AimsData<int16_t> getBrainMask(const AimsData<T> &image, byte channel)
  {
    int dx = image.dimX();
    int dy = image.dimY();
    int dz = image.dimZ();

    //  Get channel component
    ChannelSelector< AimsData<T>, AimsData<double> > selector;
    AimsData<double> G = selector.select(image, channel);

    //  EM Algorithm  TODO : remplacer par un seuillage automatique
    AimsEMThreshold< double, int16_t> em(AIMS_BETWEEN, 0);
    AimsData<int16_t> thImage = em.bin(G);

    // Hypothesis : background = 1 connected component including borders
    // Get connected components and label image
    AimsBucket<Void> component; // Components = {background , tissue, ...}
    AimsData<int16_t> C2(dx, dy, dz); // Label image
    C2 = aims::AimsLabeledConnectedComponent(
                    component,
                    thImage,                         // ne modifie pas l'image
                    aims::Connectivity::CONNECTIVITY_6_XYZ,
                    int16_t(-1),                       // background  MODIFIER
                    false,                                      // binary mode
                    (size_t) 0,                                      //minSize
                    (size_t) 0,                      //number of component Max
                    false                                            //verbose
         );

    // Seek for component which label is appearing the most within borders
    int bordure = 2;
    AimsData <int> fond (component.size());
    fond = 0;
    for (int z = 0; ((z < bordure) && (z < dz)) ; z++) // Parcours de l'image des labels
      for (int y = 0; ((y < bordure) && (y < dy)); y++)
        for (int x = 0; ((x < bordure) && (x < dx)); x++)
        {
          fond(C2(x, y, z)) ++ ;
          fond(C2(dx - 1 - x, dy - 1 - y, dz - 1 - z)) ++ ;
        }

    int ifond = 0; // indice de la composante du fond
    for (int icomp = 0; icomp < component.size(); icomp++){  // Parcours des composantes
      ifond = (fond (ifond) > fond (icomp)) ? ifond : icomp; // On garde max occurences
    }

    // Get max size of connected components (excluding background)
    int max_taille;
    bool init = true;
    for (int icomp = 0; icomp < component.size(); icomp++){  // Parcours des composantes
      std::list<AimsBucketItem<Void> > CC = component[icomp];
      if (icomp != ifond){
        if (init){
          max_taille =  CC.size();
          init = false;
        }else{
          max_taille = (max_taille > CC.size()) ? max_taille : CC.size();
        }
      }
    }

    // Add non-background components to the mask
    AimsData<int16_t> brainMask (dx, dy, dz);
    brainMask.volume()->header() = image.volume()->header();
    brainMask = 0;

    for (int icomp = 0; icomp < component.size(); icomp++){
      std::list<AimsBucketItem<Void> > CC = component[icomp];
      if ( (icomp != ifond) && (CC.size() >= 0.001 * max_taille) ) { // filtre petites composantes
        std::list<AimsBucketItem<Void> >::iterator itel, itele = CC .end();
        for (itel = CC.begin(); itel != itele; ++itel){ //parcours des buckets
          AimsBucketItem<Void> pixel = *itel;
          brainMask(pixel.location()) = 1;
        }
      }
    }

    return brainMask;
  }

} // namespace obso_temp

#endif