slices.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_DATA_SLICES_H
#define BRAINRAT_DATA_SLICES_H
 
#include <iomanip>
#include <cartobase/type/string_conversion.h>
#include <aims/io/process.h>
#include <aims/io/io_g.h>
#include <aims/data/data_g.h>
#include <aims/math/mathelem.h>
#include <aims/morphology/morphology_g.h>
#include <aims/utility/converter_g.h>
 
#define INCH             25.49     // inch->mm
#define NB_MODES         2         // Mode de number back + tissue approximation
#define MIN_HISTO        0
#define MAX_HISTO        64        // TODO test bug value 255...
#define MAX_INT16_T      32767
#define MAX_ITER         500
#define MID_WINDOW       4         // Median processed on window = 4-C-4
//#define MID_WINDOW       0         // Global median on all components
#define FRACTION         3.0       // Variation 1/N fraction 
                                   // tolerance / surf or dim [default = 3] 
#define BORDER_LOCAL     10        // Around component 
#define BORDER_GLOBAL    50        // Around component 
#define PERCENT_SURFACE  .2        // Test over mode surface detected / reject 
#define EROSION          1         // Size of erosion / separation
 
//------------------------------------------------------------------------------
// Process to detect histological slices in a scanned image
//------------------------------------------------------------------------------
class SlicesDetector : public aims::Process {
    public:
      std::string fileIn;
      std::string fileOut;
      //std::vector<int32_t> dimensionsizes, dimensionorders;
      float fraction, min_size;
      int32_t n, res;
      int32_t lt, ht;
      uint8_t c;
      bool disable_checks;
    
      SlicesDetector();
};
 
//------------------------------------------------------------------------------
// Processes to detect, extract and stack histological slices in a scanned image
//------------------------------------------------------------------------------
class SlicesStacker : public aims::Process {
    public:
      std::string fileIn;
      std::string fileInP;
      std::string fileOut;
      std::string outName;
      std::vector<int32_t> dimensionsizes, dimensionorders;
      std::string format;
      std::string background;
      bool nostack;
      float fraction, min_size;
      uint8_t sn;
      int32_t n, res;
      int32_t dx, dy, lt, ht;
      uint8_t c;
      float z;
    
    protected:
      SlicesStacker();
};
 
 
class RGBSlicesStacker : public SlicesStacker {
    public:
      RGBSlicesStacker();
};
 
class GreyLevelSlicesStacker : public SlicesStacker {
    public:
      GreyLevelSlicesStacker();
};
 
//------------------------------------------------------------------------------
// Enumerate type parameter to be tested = continuity
//------------------------------------------------------------------------------
enum ContinuityParameterTested
{
    ComponentSurface,
    ComponentSize_hz,
    ComponentSize_vt
};
 
enum Direction
{
    Horizontal,
    Vertical
};
 
enum HistoAnalysisFunction
{
    HistogramMode,
    ProjectionMode
};
 
//------------------------------------------------------------------------------
// Component structure
//------------------------------------------------------------------------------
typedef struct {
 
  int32_t label_i;   // i=init; m=inter; f=final
  int32_t label_m;   // value = -1
  int32_t label_f;
 
  int32_t abscissa_x;
  int32_t abscissa_X;
 
  int32_t ordinate_y;
  int32_t ordinate_Y;
  
  int32_t coordinate_z ;
  int32_t coordinate_Z ;
  
 
  int32_t gravity_center_x;
  int32_t gravity_center_y;
  int32_t gravity_center_z;
 
  int32_t bounding_hz;   // Extraction bounding box GC + max X
  int32_t ref_median_bounding_hz;
  int32_t test_median_bounding_hz;
 
  int32_t bounding_vt;
  int32_t ref_median_bounding_vt;
  int32_t test_median_bounding_vt;
 
  int32_t surface;
  int32_t surfaceMM;
  int32_t ref_median_surface;
  int32_t test_median_surface;
  
  float   diff_median_bounding_hz;
  float   diff_median_bounding_vt;
  float   diff_median_surface;
 
} Component, *ComponentPtr;
 
//------------------------------------------------------------------------------
// Validation_Test structure 
//------------------------------------------------------------------------------
typedef struct {
  // Mandatory
  bool test_nb_component;
 
  // Optional
  bool test_size_continuity;
  bool test_horizontal_continuity;
  bool test_vertical_continuity;
 
} Tests, *TestsPtr;
 
//------------------------------------------------------------------------------
// Maxima histogram structure
//------------------------------------------------------------------------------
typedef struct {
 
  int32_t iteration;
  int32_t numero;
  int32_t value;
  int32_t position;
  int32_t surface;
  int32_t test_surface;
  int32_t extrema;
 
} Maxima, *MaximaPtr;
 
//------------------------------------------------------------------------------
// HistologyScanSlicesDetector
//------------------------------------------------------------------------------
class HistologyScanSlicesDetector {   
public:
    
    // Constructor fo class HistologyScanSlicesDetector
    HistologyScanSlicesDetector(
        float   slice_min_surface = 1.,
        float   slice_median_variation = 1. / FRACTION,
        int32_t histogram_modes = NB_MODES,
        int32_t histogram_bins = MAX_HISTO,
        bool    check_slices = true);
    
    // Analyzes histology scan and return a generic object containing
    // bounding boxes about detected slices
    carto::Object boundingBoxes(const carto::VolumeRef<int16_t> & vol, 
                                const std::vector<Component> & components);
 
    std::vector<Component> analyze(carto::VolumeRef<int16_t> & vol,
                                   const int32_t slices_number = 6,
                                   const int32_t scan_resolution = 0,
                                   const int32_t threshold_lowerbound = 0, 
                                   const int32_t threshold_upperbound = 0);
    
protected:
    
    bool detectHistogramModes(const carto::VolumeRef<int32_t> & histogram,
                              const int32_t slices_number,
                              const int32_t slices_min_voxels,
                              const HistoAnalysisFunction analysis_function,
                              int32_t & background,
                              int32_t & slices_mode_lowerbound,
                              int32_t & slices_mode_upperbound);
    
    std::vector<Component> detectComponents(
                              carto::VolumeRef<int16_t> & vol, 
                              const int32_t slices_number,
                              const int32_t threshold_lowerbound, 
                              const int32_t threshold_upperbound);
    
    std::vector<Component> checkComponents(
        carto::VolumeRef<int16_t> & vol,
        std::vector<Component> &components,
        const int32_t slices_number);
    
    std::vector<Component> processOutlierComponents(
        std::vector<Component> &components,
        carto::VolumeRef<int16_t> &vol,
        const int32_t slices_number);
    
    std::vector<Component> openComponent(
        const std::vector<Component> & components, 
        const int32_t component, 
        carto::VolumeRef<int16_t> &vol);
 
    std::vector<Component> separateComponent(
        const std::vector<Component> & components, 
        const int32_t component, 
        carto::VolumeRef<int16_t> &vol,
        const Direction dir);
    
    void writeTestsWarnings(const Tests &t);
    
    //int64_t getSurfaceVariationInMm();
    
private:
    
//    carto::VolumeRef<T> _vol;
//    int32_t _slices_number;
    float   _slice_min_surface;
    float   _slice_median_variation;
    int32_t _histogram_modes; 
    int32_t _histogram_bins; 
    bool    _checks_slices;
//     std::vector<int32_t> _dimension_sizes, 
//                          _dimension_orders;
};
 
//------------------------------------------------------------------------------
// Iterator type definitions MAXIMA / COMPONENT 
//------------------------------------------------------------------------------
typedef std::vector<Component>::iterator CompIterator;
typedef std::vector<Maxima>::iterator MaxIterator;
 
Maxima affectMaxima(int32_t iter, int32_t num, int32_t val, 
                    int32_t pos, int32_t extr);
 
//------------------------------------------------------------------------------
// Write Component + vector<Point3d)
//------------------------------------------------------------------------------
void writeComponents(std::vector<Component> &comp);
void writeVector3int32_t(std::vector<AimsVector<int32_t, 3> > &vect);
void writeVectorInt32_t(std::vector<int32_t> &vect);
void writeAimsInt32_t(carto::VolumeRef<int32_t> &vol);
void writeAimsFloat(carto::VolumeRef<float> &vol);
void writeMaxima(std::vector<Maxima> &mxm);
void writeMaximaIter(std::vector<Maxima> &mxm,
                      int32_t iter);
 
//------------------------------------------------------------------------------
// Functions to integrate to StackSlices
//------------------------------------------------------------------------------
 
//------------------------------------------------------------------------------
// 0) Analyse histogram to extract component mode  
//------------------------------------------------------------------------------
template <class T>
float computeHistogramRatio(const carto::VolumeRef<T> &vol, 
                            const int32_t bins);
 
carto::VolumeRef<int32_t> generateHistogram(
    const carto::VolumeRef<int16_t> &vol,
    const int32_t bins = MAX_HISTO);
 
std::vector<AimsVector<int32_t, 3> > updateVDM(
    std::vector<AimsVector<int32_t, 3> > &vdm);
 
carto::VolumeRef<int16_t> applyThreshold(
    const carto::VolumeRef<int16_t> &vol, 
    const int16_t threshold_lower_bound, 
    const int16_t threshold_upper_bound,
    const int16_t new_value = MAX_INT16_T);
 
//------------------------------------------------------------------------------
// 1) Connex component extraction
//------------------------------------------------------------------------------
carto::VolumeRef<int16_t> computeConnexComponent(
    carto::VolumeRef<int16_t> &vol);
 
//------------------------------------------------------------------------------
// 2) Extract information from component
//------------------------------------------------------------------------------
std::vector<Component> removeComponent(std::vector<Component> cpnt, int32_t l);
 
std::vector<Component> removeComponents(
    const std::vector<Component> & cpnt, 
    const std::vector<Component> & to_remove);
 
std::vector<Component> insertComponent(std::vector<Component> cpnt, 
                                       int32_t comp,
                                       std::vector<Component> cpnt2slices, 
                                       int32_t l);
 
std::vector<int32_t> getComponentLabels(const std::vector<Component> & cpnt);
 
 
std::vector<Component> extractInformationFromComponent(
    const carto::VolumeRef<int16_t> &vol);
 
//------------------------------------------------------------------------------
// 3) Test function: 1) init and 2) update
//------------------------------------------------------------------------------
 
Tests testComponents(const std::vector<Component> &cpnt, int32_t number_cpnt);
 
void displayComponentsStrategy(const std::vector<Component> &cpnt);
 
void displayStatus(const Tests &t);
 
bool isValid(const Tests &t);
 
//------------------------------------------------------------------------------
// 4) Test size && dimension X/Y continuity
//------------------------------------------------------------------------------
// Point3d return=value, mean, std
std::vector<int32_t> extractVectorFromComponent(
    std::vector<Component> cpnt, 
    ContinuityParameterTested param);
 
std::vector<Component> processMedianConsistencyParameters(
    std::vector<Component> &cpnt, float variation, int32_t mid_window = 0);
 
//------------------------------------------------------------------------------
// 4 bis) Process data with morphomath 
//------------------------------------------------------------------------------ 
carto::VolumeRef<int32_t> projectData(carto::VolumeRef<int16_t> &X);
 
carto::VolumeRef<int16_t> transposeVolume(
    carto::VolumeRef<int16_t> &vol);
 
bool removeComponentsInVolume(carto::VolumeRef<int16_t> & vol, 
                              const std::vector<Component> & to_remove);
 
carto::VolumeRef<int16_t> extractFromVolume(
    const carto::VolumeRef<int16_t> &vol,
    const int32_t x,
    const int32_t X,
    const int32_t y,
    const int32_t Y,
    const int32_t old_label,
    const int32_t new_label);
 
bool copyToVolume(const carto::VolumeRef<int16_t> &vol,
                  carto::VolumeRef<int16_t> &destvol,
                  const int32_t x,
                  const int32_t X,
                  const int32_t y,
                  const int32_t Y,
                  const int32_t old_label,
                  const int32_t new_label);
 
std::vector<Component> openingComponent(std::vector<Component> &cpnt, 
                                        int32_t comp, 
                                        carto::VolumeRef<int16_t> vol,
                                        float fraction);
 
std::vector<Component> separateComponent(std::vector<Component> &cpnt, 
                                         int32_t comp, 
                                         carto::VolumeRef<int16_t> vol,
                                         int32_t min_surface,
                                         Direction dir,
                                         float fraction);
 
std::vector<Component> getOutliersBySurface(const std::vector<Component> & cpnt, 
                                            const float min_surface);
 
std::vector<Component> getOutliersByWeight(const std::vector<Component> & cpnt, 
                                           const int32_t slices_number);
 
std::vector<Component> filterComponentsBySurface(
    const std::vector<Component> & cpnt, 
    const float min_surface);
 
std::vector<Component> weightedRemoveComponents(
                                         std::vector<Component> &cpnt, 
                                         int32_t slices_number);
 
std::vector<Component> processOutlierComponent(std::vector<Component> &cpnt, 
                                               carto::VolumeRef<int16_t> &vol,
                                               int32_t min_surface,
                                               int32_t slices_number,
                                               float fraction);
 
//------------------------------------------------------------------------------
// 5) Sort component according to Right->Left / Top->Down
//------------------------------------------------------------------------------
int32_t getDimensionCount(std::vector<int32_t> dimensionsizes);
 
std::vector<Component> changeComponentsIndex(
    std::vector<Component>& components,
    std::vector<int32_t> dimensionsizes, 
    std::vector<int32_t> dimensionorders, 
    int32_t offset = 1, int32_t step = 1);
 
bool compareComponentsByLabelM(const Component & c1, const Component & c2);
 
std::vector<Component> sortComponent(std::vector<Component> & cpnt,
                                     std::vector<int32_t> dimensionsizes);
 
//------------------------------------------------------------------------------
// 6) Apply new labelling
//------------------------------------------------------------------------------ 
carto::VolumeRef<int16_t> applyNewLabelling(std::vector<Component> cpnt,
                                    carto::VolumeRef<int16_t> vol);
 
//------------------------------------------------------------------------------
// 7) Dectect slices
//------------------------------------------------------------------------------ 
template<class T> bool 
detectSlices(aims::Process & p, const std::string &, aims::Finder &);
 
//------------------------------------------------------------------------------
// 8) Reconstruct volume
//------------------------------------------------------------------------------
int32_t extractMaxima(std::vector<int32_t> vctr);
 
template<class T>
void rescale(carto::VolumeRef<T> &in, carto::RescalerInfo & info);
 
template<class T>
bool writeStackedSlices(int32_t dx, int32_t dy, int32_t res, 
                         float z,
                         T background,
                         std::vector<Component> &cpnt,
                         carto::VolumeRef<T> &vol,
                         std::string outputfile,
                         std::string format);
 
template<class T>
bool writeNumberedSlices(int32_t dx, int32_t dy, int32_t res, 
                          float z,
                          T background,
                          std::vector<Component> &cpnt,
                          carto::VolumeRef<T> &vol,
                          std::string outputdir,
                          std::string outputname,
                          std::string format,
                          uint8_t slicemin);
 
template<class T>
static bool individStackSlices(aims::Process & p, 
                               const std::string &, 
                               aims::Finder &);
 
 
//------------------------------------------------------------------------------
// Utilities
//------------------------------------------------------------------------------
template <class T>
float getVoxelLebesgue(const carto::VolumeRef<T> & vol);
 
template <class T>
int32_t getVoxelCount(const carto::VolumeRef<T> & vol, float value);
 
template<class T>
inline T background(std::vector<Component> & cpnt, carto::VolumeRef<T> &vol)
{
    double bg = 0;
    int32_t count = 0;
 
    for (int32_t i = 0; i <(int32_t)cpnt.size(); i++) 
    {
        // Uses pixels around the component bounding box to process background
        // mean value
        Component & c = cpnt[ i ];
 
        if (((c.abscissa_x - BORDER_LOCAL)> 0) 
         && ((c.abscissa_X + BORDER_LOCAL) <= vol.getSizeX()))
        {
            for (int32_t y = c.ordinate_y; y <c.ordinate_Y; y++)
            {
                bg += vol(c.abscissa_x - 1, y);
                bg += vol(c.abscissa_X, y);
            }
 
            count += (c.bounding_vt * 2);
        }
 
        if (((c.ordinate_y - BORDER_LOCAL)> 0)
         && ((c.ordinate_Y  + BORDER_LOCAL) <= vol.getSizeY()))
        {
            for (int32_t x = c.abscissa_x; x <c.abscissa_X; x++)
            {
                bg += vol(x, c.ordinate_y - 1);
                bg += vol(x, c.ordinate_Y);
            }
            count += (c.bounding_hz * 2);
        }
    }
    
    return (T)(bg / count);
}
 
template<>
inline AimsRGB background<AimsRGB>(std::vector<Component> & cpnt, 
                                   carto::VolumeRef<AimsRGB> &vol)
{
    AimsRGB bg(0,0,0);
    AimsVector<double, 3> s(0,0,0);
    int32_t count = 0;
 
    for (int32_t i = 0; i <(int32_t)cpnt.size(); i++) 
    {
        // Uses pixels around the component bounding box to process background
        // mean value
        Component & c = cpnt[ i ];
 
        if (((c.abscissa_x - BORDER_LOCAL)> 0)
         && ((c.abscissa_X + BORDER_LOCAL) <= vol.getSizeX()))
        {
            for (int32_t y = c.ordinate_y; y <c.ordinate_Y; y++)
            {
                for (int32_t v = 0; v <3; v++)
                {
                    s[v] += vol(c.abscissa_x - 1, y)[v];
                    s[v] += vol(c.abscissa_X, y)[v];
                }
            }
 
            count += (c.bounding_vt * 2);
        }
 
        if (((c.ordinate_y - BORDER_LOCAL)> 0)
         && ((c.ordinate_Y  + BORDER_LOCAL) <= vol.getSizeY()))
        {
            for (int32_t x = c.abscissa_x; x <c.abscissa_X; x++)
            {
                for (int32_t v = 0; v <3; v++)
                {
                    s[v] += vol(x, c.ordinate_y - 1)[v];
                    s[v] += vol(x, c.ordinate_Y)[v];
                }
            }
            count += (c.bounding_hz * 2);
        }
    }
    
    bg[0] = (uint8_t)(s[0] / count);
    bg[1] = (uint8_t)(s[1] / count);
    bg[2] = (uint8_t)(s[2] / count);
    
    return (AimsRGB)(bg);
}
 
template<>
inline AimsRGBA background<AimsRGBA>(std::vector<Component> & cpnt,
                                     carto::VolumeRef<AimsRGBA> &vol)
{
    AimsRGBA bg(0,0,0,0);
    AimsVector<double, 4> s(0,0,0,0);
    int32_t count = 0;
 
    for (int32_t i = 0; i <(int32_t)cpnt.size(); i++) 
    {
        // Uses pixels around the component bounding box to process background
        // mean value
        Component & c = cpnt[ i ];
 
        if (((c.abscissa_x - BORDER_LOCAL)> 0)
         && ((c.abscissa_X + BORDER_LOCAL) <= vol.getSizeX()))
        {
            for (int32_t y = c.ordinate_y; y <c.ordinate_Y; y++)
            {
                for (int32_t v = 0; v <4; v++)
                {
                    s[v] += vol(c.abscissa_x - 1, y)[v];
                    s[v] += vol(c.abscissa_X, y)[v];
                }
            }
            count += (c.bounding_vt * 2);
        }
 
        if (((c.ordinate_y - BORDER_LOCAL)> 0)
         && ((c.ordinate_Y  + BORDER_LOCAL) <= vol.getSizeY()))
        {
            for (int32_t x = c.abscissa_x; x <c.abscissa_X; x++)
            {
                for (int32_t v = 0; v <4; v++)
                {
                    s[v] += vol(x, c.ordinate_y - 1)[v];
                    s[v] += vol(x, c.ordinate_Y)[v];
                }
            }
            count += (c.bounding_hz * 2);
        }
    }
    
    bg[0] = (uint8_t)(bg[0] / count);
    bg[1] = (uint8_t)(bg[1] / count);
    bg[2] = (uint8_t)(bg[2] / count);
    bg[3] = (uint8_t)(bg[3] / count);
    
    return (AimsRGBA)(bg);
}
 
#endif