slices_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
 */
 
#include <cartobase/object/object.h>
#include <cartobase/stream/directory.h>
#include <cartobase/stream/fileutil.h>
#include <aims/io/channelR.h>
#include <brainrat/data/slices.h>
#include <brainrat/exception/stackerexcept.h>
 
 
//------------------------------------------------------------------------------
// Utility functions
//------------------------------------------------------------------------------
template<class T> bool 
detectSlices(aims::Process & p, const std::string &, aims::Finder &)
{
    SlicesDetector & proc = (SlicesDetector &)p;
    
    // Instanciate detector object
    HistologyScanSlicesDetector detector(proc.min_size,
                                         1. / proc.fraction,
                                         NB_MODES,
                                         MAX_HISTO,
                                         !proc.disable_checks);
 
//     std::vector<Component> components;
    carto::Object bounds;
    
//     float default_value = (proc.res? INCH / proc.res > 0 : 1.);
//     std::vector<float> vs(default_value, 4);
 
    carto::VolumeRef<int16_t> vol;
    //
    // Read data to use for processing components
    // 
    aims::ChannelReader<carto::VolumeRef<int16_t> > reader(proc.fileIn);
    reader.read(vol, proc.c);
    std::cout << "Original image loaded" << std::endl;
    
    // Step 1: Detect components
    std::vector<Component> components = detector.analyze(vol, proc.n, proc.res,
                                                         proc.lt, proc.ht);
    
    // Step 2: Convert components to generic object that describes bounding
    // boxes
    bounds = detector.boundingBoxes(vol, components);
 
//     // Step 2: sort components
//     // It is necessary to first sort using 2 dimensions because images are 2d
//     if (proc.dimensionsizes.size() == 0) {
//         proc.dimensionsizes.push_back(1);
//     }
// 
//     if (proc.dimensionorders.size() == 0) {
//         proc.dimensionorders.push_back(0);
//     }
// 
//     std::vector<int32_t> columnsizes;
//     columnsizes.push_back(proc.dimensionsizes[0]);
//     columnsizes.push_back(proc.dimensionsizes[1] * proc.dimensionsizes[2]);
//     columnsizes.push_back(1);
//     // std::cout << "First sorting process." << std::endl;
//     components = sortComponent(components, columnsizes);
//     // std::cout << "Second sorting process." << std::endl;
//     components = changeComponentsIndex(components, 
//                                        proc.dimensionsizes, 
//                                        proc.dimensionorders, 
//                                        1, 1);
//         
//     std::cout << "Last sorting process (according to slice taking)." 
//               << std::endl;
//     std::cout << std::endl;
//     writeComponents(components);
//     std::cout << std::endl;
    
    // Write components in JSON format
//     carto::Object bounds = carto::Object::value(carto::PropertySet());
//     std::vector<carto::Object> b;
//     
//     std::vector<Component>::const_iterator end = components.end(), it;
//     for (it = components.begin(); it != end; ++it) {
//         carto::Object bound = carto::Object::value(carto::PropertySet());
//         
//         std::vector<std::vector<float> > v;
//         std::vector<float> lower_bound(3, 0.);
//         std::vector<float> upper_bound(3, 0.);
//         const Component & c = *it;
//         
//         // Get bounding box coordinates
//         lower_bound[0] = (vs[0] * c.abscissa_x);
//         lower_bound[1] = (vs[1] * c.ordinate_y);
//         lower_bound[2] = (vs[2] * c.coordinate_z);
//         upper_bound[0] = (vs[0] * c.abscissa_X);
//         upper_bound[1] = (vs[1] * c.ordinate_Y);
//         upper_bound[2] = (vs[2] * c.coordinate_Z);
//         v.push_back(lower_bound);
//         v.push_back(upper_bound);
//         
//         bound->setProperty("bounds", v);
//         b.push_back(bound);
//     }
//     bounds->setProperty("bounding_boxes", b);
    
    aims::Writer<carto::Object> w(proc.fileOut);
    const std::string format("JSON");
    return w.write(bounds, false, &format);
}
 
template<class T> bool 
individStackSlices(aims::Process & p, const std::string &, aims::Finder &)
{
    SlicesStacker & proc = (SlicesStacker &)p;
 
// #ifdef INDIVID_STACK_SLICES
    carto::VolumeRef<T> vol;
    
    std::vector<Component> components;
 
    if (proc.dimensionsizes.size() == 0) {
        proc.dimensionsizes.push_back(1);
    }
 
    if (proc.dimensionorders.size() == 0) {
        proc.dimensionorders.push_back(0);
    }
    
    
    {
        // Only allocate data_analysis during analysis step
        carto::VolumeRef<int16_t> volgrey;
        //
        // Read data to use for processing components
        // 
        if (!proc.fileInP.empty())
        {
            aims::ChannelReader<carto::VolumeRef<int16_t> > reader(proc.fileInP);
            reader.read(volgrey, proc.c);
            std::cout << "Pre-processed image loaded" << std::endl;
        }
        else
        {
            aims::ChannelReader<carto::VolumeRef<int16_t> > reader(proc.fileIn);
            reader.read(volgrey, proc.c);
            std::cout << "Original image loaded" << std::endl;
        }
        // Instanciate detector object
        HistologyScanSlicesDetector detector(proc.min_size,
                                             1. / proc.fraction,
                                             NB_MODES,
                                             MAX_HISTO,
                                             true);
            
        // Step 1: Detect components
        components = detector.analyze(volgrey, proc.n, proc.res, 
                                      proc.lt, proc.ht);
    }
 
    // Step 2: sort components
    // It is necessary to first sort using 2 dimensions because images are 2d
    std::vector<int32_t> columnsizes;
    columnsizes.push_back(proc.dimensionsizes[0]);
    columnsizes.push_back(proc.dimensionsizes[1] * proc.dimensionsizes[2]);
    columnsizes.push_back(1);
    // std::cout << "First sorting process." << std::endl;
    components = sortComponent(components, columnsizes);
    // std::cout << "Second sorting process." << std::endl;
    components = changeComponentsIndex(components, 
                                       proc.dimensionsizes, 
                                       proc.dimensionorders, 
                                       1, 1);
        
    std::cout << "Last sort processed (according to slice taking) ..." 
              << std::endl;
              
    std::cout << std::endl;
    writeComponents(components);
    std::cout << std::endl;
        
    // Step 3: reconstruct volume
    aims::Reader<carto::VolumeRef<T> > fileIn(proc.fileIn);
    fileIn.read(vol);
 
    T bg = background<T>(components, vol);
    std::cout << "Background : " << toString(bg) << std::endl;
    
    bool result = false;
    if (proc.nostack)
    {
        // Use input file name by default
        std::string outname = (proc.outName.empty() ? FileUtil::removeExtension(
                                                    FileUtil::basename(
                                                        proc.fileIn))
                                               : proc.outName);
 
        result = writeNumberedSlices(proc.dx, proc.dy, 
                                     proc.res, proc.z,
                                     bg, 
                                     components, 
                                     vol,
                                     proc.fileOut,
                                     outname,
                                     proc.format,
                                     proc.sn);
        std::cout << "Volume written to output directory. Terminated."
                  << std::endl;
    }
    else
    {
        result = writeStackedSlices(proc.dx, proc.dy, 
                                    proc.res, proc.z,
                                    bg, 
                                    components, 
                                    vol,
                                    proc.fileOut,
                                    proc.format);
        std::cout << "Volume written to output file. Terminated." << std::endl;
    }
 
    return result;
// #else
//     std::cout << "!!!WARNING!!! indiviualize and stack processing has been "
//                  "disabled at compile time." << std::endl;    
// #endif
}
 
 
template<class T>
void rescale(carto::VolumeRef<T> &in, RescalerInfo & info)
{
    if (std::isnan(info.vmin)) {
        if (! info.usevtypelimits) {
            info.vmin = (double) in.min();
        }
    }
 
    if (std::isnan(info.vmax)) {
        if (! info.usevtypelimits) {
            info.vmax = (double) in.max();
        }
    }
 
    DefaultedRescalerInfo<T, T> defaultedinfo(info);
 
    int x, y, z, t, 
        dx = in.getSizeX(), 
        dy = in.getSizeY(), 
        dz = in.getSizeZ(), 
        dt = in.getSizeT();
 
    for (t = 0; t < dt; ++t)
    {
        for (z = 0; z < dz; ++z)
        {
            for (y = 0; y < dy; ++y)
            {
                for (x = 0; x < dx; ++x)
                    in(x, y, z, t) = defaultedinfo.getScaledValue(
                                        in(x, y, z, t));
            }
        }
    }
}
 
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> &data,
                        std::string outputfile,
                        std::string format)
{
    int32_t x, y, l;
 
    const std::string *wf = 0;
    if (!format.empty())
        wf = &format;
 
    Object options = Object::value(Dictionary());
    //options->setProperty("exact_format", true);
    
    if (dx == 0)
        dx = extractMaxima(
             extractVectorFromComponent(cpnt, ComponentSize_hz)) + 2 
             * (BORDER_LOCAL + BORDER_GLOBAL);
    else if (dx < extractMaxima(
                  extractVectorFromComponent(cpnt, ComponentSize_hz)) + 2 
                  * BORDER_LOCAL)
    {
        throw dimension_x_error();
    }
    
    if (dy == 0)
        dy = extractMaxima(
             extractVectorFromComponent(cpnt, ComponentSize_vt)) + 2 
             * (BORDER_LOCAL + BORDER_GLOBAL);
    else if (dy < extractMaxima(
                  extractVectorFromComponent(cpnt, ComponentSize_vt)) + 2 
                  * BORDER_LOCAL) 
    {
        throw dimension_y_error();
    }
    
    carto::VolumeRef<T> recons_volume(dx, dy, cpnt.size());
    std::vector<float> vs(4, 1.);
    vs[0] = INCH / res;
    vs[1] = INCH / res;
    vs[2] = z;
    
    recons_volume.header().setProperty("voxel_size", vs);
    
    //   Background optimized from histogram analysis
    recons_volume = background;
 
    std::cout << "Reconstructing volume ..." << std::endl << std::endl;  
    for (l = 0; l < int32_t(cpnt.size()); ++l)
    {
        Point2d shift, size_component;
        
        size_component[0] = cpnt[l].bounding_hz + 2 * BORDER_LOCAL;
        size_component[1] = cpnt[l].bounding_vt + 2 * BORDER_LOCAL;
        
        shift[0] = (dx - size_component[0]) / 2;
        shift[1] = (dy - size_component[1]) / 2;
        
        for (y = 0; y < size_component[1] ; ++y) 
            for (x = 0; x < size_component[0] ; ++x)
            {
                if (((x + cpnt[l].abscissa_x - BORDER_LOCAL) >= 0)
                 && ((x + cpnt[l].abscissa_x - BORDER_LOCAL) < data.getSizeX())
                 && ((y + cpnt[l].ordinate_y - BORDER_LOCAL) >= 0)
                 && ((y + cpnt[l].ordinate_y - BORDER_LOCAL) < data.getSizeY()))
                {
                    recons_volume(x + shift[0], 
                                  y + shift[1], 
                                  cpnt[l].label_f - 1) = data(
                                      x + cpnt[l].abscissa_x - BORDER_LOCAL,
                                      y + cpnt[l].ordinate_y - BORDER_LOCAL);
                }
            }
    }
 
    std::cout << "Stacking process achieved." << std::endl;
 
    //
    // Ecriture du fichier
    //
    aims::Writer<carto::VolumeRef<T> > volumeWriter(outputfile, options);
    return volumeWriter.write(recons_volume, false, wf);
}
 
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> &data,
                         std::string outputdir,
                         std::string outputname,
                         std::string format,
                         uint8_t slicemin)
{
    int32_t x, y, l;
    
    Directory dir = Directory(outputdir);
    dir.makedirs();
 
    std::string dirname = dir.dirname();
    std::string extension = FileUtil::extension(outputname);
    std::string wf = (format.empty() && extension.empty() ? "GIS" : format);
 
    Object options = Object::value(Dictionary());
    options->setProperty("exact_format", true);
 
    if (dx == 0)
        dx = extractMaxima(
             extractVectorFromComponent(cpnt, ComponentSize_hz)) + 2 
             * (BORDER_LOCAL + BORDER_GLOBAL);
    else if (dx < extractMaxima(
                  extractVectorFromComponent(cpnt, ComponentSize_hz)) + 2
                  * BORDER_LOCAL) 
    {
        throw dimension_x_error();
    }
    
    if (dy == 0)
        dy = extractMaxima(
             extractVectorFromComponent(cpnt, ComponentSize_vt)) + 2
             * (BORDER_LOCAL + BORDER_GLOBAL);
    else if (dy < extractMaxima(
                  extractVectorFromComponent(cpnt, ComponentSize_vt)) + 2
                  * BORDER_LOCAL)
    {
        throw dimension_y_error();
    }
    
    carto::VolumeRef<T> slice(dx, dy);
    std::vector<float> vs(4, 1.);
    vs[0] = INCH / res;
    vs[1] = INCH / res;
    vs[2] = z;
    
    slice.header().setProperty("voxel_size", vs);
    
    std::cout << "Writing volume slices ..." << std::endl << std::endl;
 
    for (l=0; l<int32_t(cpnt.size()); ++l)
    {
        // Background optimized from histogram analysis
        slice = background;
 
        Point2d shift, size_component;
        
        size_component[0] = cpnt[l].bounding_hz + 2 * BORDER_LOCAL;
        size_component[1] = cpnt[l].bounding_vt + 2 * BORDER_LOCAL;
        
        shift[0] = (dx - size_component[0]) / 2;
        shift[1] = (dy - size_component[1]) / 2;
        
        for (y = 0; y < size_component[1] ; ++y) 
            for (x = 0; x < size_component[0] ; ++x)
            {
                if (((x + cpnt[l].abscissa_x - BORDER_LOCAL) >= 0)
                 && ((x + cpnt[l].abscissa_x - BORDER_LOCAL) < data.getSizeX())
                 && ((y + cpnt[l].ordinate_y - BORDER_LOCAL) >= 0)
                 && ((y + cpnt[l].ordinate_y - BORDER_LOCAL) < data.getSizeY()))
                {
                    slice(x + shift[0], 
                          y + shift[1]) = data(x + cpnt[l].abscissa_x
                                               - BORDER_LOCAL,
                                               y + cpnt[l].ordinate_y
                                               - BORDER_LOCAL);
                }
            }
 
        std::stringstream ss;
        ss << outputdir << FileUtil::separator() 
           << FileUtil::removeExtension(FileUtil::basename(outputname));
        ss << "_" << setfill('0') << setw(4) << slicemin + cpnt[l].label_f - 1;
 
        if (! extension.empty())
            ss << "." << extension;
 
        // std::cout << "Writing file : " << ss.str() 
        //           << std::endl << std::flush;
        
        aims::Writer<carto::VolumeRef<T> > sliceWriter(ss.str(), options);
        if (! sliceWriter.write(slice, false, &wf))
            return false;
    }
 
    return true;
}