interesttiling.cc

/*
* This file is part of rasdaman community.
*
* Rasdaman community is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Rasdaman community is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with rasdaman community.  If not, see <http://www.gnu.org/licenses/>.
*
* Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Peter Baumann /
rasdaman GmbH.
*
* For more information please see <http://www.rasdaman.org>
* or contact Peter Baumann via <baumann@rasdaman.com>.
/
/**
 * SOURCE: interesttiling.cc
 *
 * MODULE: rasodmg
 * CLASS:  r_Interest_Tiling
 *
 * COMMENTS:
 *      None
*/

#include "rasodmg/interesttiling.hh"
#include "rasodmg/alignedtiling.hh"
#include "rasodmg/dirdecompose.hh"
#include "rasodmg/dirtiling.hh"
#include "raslib/rminit.hh"

#include <assert.h>
#include <string.h>
#include <math.h>
#include <cstdlib>

#ifdef AIX
#include <strings.h>
#endif

// Uncoment the _VISUALIZE_2D_DECOMP_ line to generate ppm files the
// visualization of the domain decomposition done by the algoritm
// #define _VISUALIZE_2D_DECOMP_

// Uncoment the folling line for some debug printfs on this class
// #define _DEBUG_INTERESTTILING_

#ifdef _VISUALIZE_2D_DECOMP_
#include "tools/visualtiling2d.hh"
#endif

// This is a data structure for internal use within interesttiling.cc.
// It is defined here because of problems with ptrepository.
// -------------
// This structure consists of a r_Minterval and a counter. Also, it has
// some auxiliary constructors.

class Classified_Block
{
public:
    int intersection_count;                         // Intersection count
    r_Minterval block;                              // The actual block

    // Default constructor
    Classified_Block(int count = 0);

    // Constructor with the actual block and the counter
    Classified_Block(const r_Minterval& b, int count = 0);

    // Same data structure (this operator is needed because of dlist)
    bool operator==(const Classified_Block& other) const;

    // Different data structure (this operator is needed because of dlist)
    bool operator!=(const Classified_Block& other) const;

    // Friend of std::ostream (this operator is needed because of dlist)
    friend std::ostream& operator<<(std::ostream& os, const Classified_Block& block);
};

const char*
r_Interest_Tiling::description = "dimensions, areas of interest, tile size (in bytes) and tile size limit [NOLIMIT|REGROUP|SUBTILING|REGROUPSUBTILING] (ex: \"2;[0:9,0:9];[100:109,0:9];100;REGROUPSUBTILING\")";

const char* r_Interest_Tiling::tilesizelimit_name_nolimit       ="NOLIMIT";
const char* r_Interest_Tiling::tilesizelimit_name_regroup       ="REGROUP";
const char* r_Interest_Tiling::tilesizelimit_name_subtiling     ="SUBTILING";
const char* r_Interest_Tiling::tilesizelimit_name_regroupandsubtiling   ="REGROUPSUBTILING";
const char* r_Interest_Tiling::all_tilesizelimit_names[r_Interest_Tiling::NUMBER]=
{
    tilesizelimit_name_nolimit,
    tilesizelimit_name_regroup,
    tilesizelimit_name_subtiling,
    tilesizelimit_name_regroupandsubtiling
};

r_Interest_Tiling::Tilesize_Limit
r_Interest_Tiling::get_tilesize_limit_from_name(const char* name)
{

    if(!name)
    {
        RMInit::logOut << "r_Interest_Tiling::get_tilesize_limit_from_name(" << (name?name: "NULL") << ")" << endl;
        return r_Interest_Tiling::NUMBER;
    }

    unsigned int i=r_Interest_Tiling::NUMBER;

    for (i=0; i<(unsigned int)r_Interest_Tiling::NUMBER; i++)
    {
        if (strcasecmp(name, all_tilesizelimit_names[i]) == 0)
            break;
    }
    return (r_Interest_Tiling::Tilesize_Limit)i;
}

const char*
r_Interest_Tiling::get_name_from_tilesize_limit(Tilesize_Limit tsl)
{
    static const char* unknown="UNKNOWN";
    unsigned int idx = (unsigned int)tsl;

    if (idx >= (unsigned int)r_Interest_Tiling::NUMBER)
        return unknown;

    return all_tilesizelimit_names[idx];
}

r_Interest_Tiling::r_Interest_Tiling(const char* encoded) throw (r_Error)
    :   r_Dimension_Tiling(0,0)
{
    if(!encoded)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << (encoded?encoded: "NULL") << ")" << endl;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

    std::vector<r_Minterval> vectInterestAreas;
    r_Interest_Tiling::Tilesize_Limit tileSizeLimit;
    r_Dimension tileD=0;
    r_Bytes tileS=0, lenToConvert=0;
    const char *pStart=NULL, *pRes=NULL, *pTemp=NULL, *pEnd=NULL;
    char *pToConvert=NULL;


    pStart=encoded;
    pEnd=pStart+strlen(pStart);
    pTemp=pStart;
    pRes=strstr(pTemp,TCOLON);

    if(!pRes)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile dimension from \"" << pTemp << "\"." << endl;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

//deal with dimension
    lenToConvert=pRes-pTemp;
    pToConvert=new char[lenToConvert+1];
    memcpy(pToConvert,pTemp, lenToConvert);
    pToConvert[lenToConvert]='\0';

    tileD=strtol(pToConvert, (char**)NULL, DefaultBase);
    if (!tileD)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile dimension from \"" << pToConvert << "\"." << endl;
        delete[] pToConvert;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }
    if (tileD < 0)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile dimension from \"" << pToConvert << "\", is negative." << endl;
        delete[] pToConvert;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

//skip COLON && free buffer
    delete[] pToConvert;
    if(pRes != (pEnd-1))
        pRes++;
    else
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas, end of stream." << endl;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

//parse interest areas
    pTemp=pRes;
    pRes=strstr(pTemp,TCOLON);

    if(!pRes)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas." << endl;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

    while(pRes)
    {
        //is interest areas?
        if(*pTemp != *LSQRBRA)
            break;

        //copy parsed interest area
        lenToConvert=pRes-pTemp;
        pToConvert=new char[lenToConvert+1];
        memcpy(pToConvert, pTemp, lenToConvert);
        pToConvert[lenToConvert]='\0';

        //add to vector
        try
        {
            r_Minterval a(pToConvert);
            vectInterestAreas.push_back(a);
        }
        catch(r_Error& err)
        {
            RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest area from \"" << pToConvert << "\"." << endl;
            RMInit::logOut << "Error " << err.get_errorno() << " : " << err.what() << endl;
            delete [] pToConvert;
            throw r_Error(TILINGPARAMETERNOTCORRECT);
        }

        //skip COLON
        delete[] pToConvert;
        if(pRes != (pEnd-1))
            pRes++;
        else
        {
            RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas, end of stream." << endl;
            throw r_Error(TILINGPARAMETERNOTCORRECT);
        }

        //try next item
        pTemp=pRes;
        pRes=strstr(pTemp, TCOLON);
    }

    if(vectInterestAreas.empty())
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding interest areas, no interest areas specified." << endl;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

//deal with tile size
    lenToConvert=pRes-pTemp;
    pToConvert=new char[lenToConvert+1];
    memcpy(pToConvert, pTemp, lenToConvert);
    pToConvert[lenToConvert]='\0';

    tileS=strtol(pToConvert, (char**)NULL, DefaultBase);
    if (!tileS)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size from \"" << pToConvert << "\"." << endl;
        delete[] pToConvert;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }
    if (tileS < 0)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size from \"" << pToConvert << "\", is negative." << endl;
        delete[] pToConvert;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

//skip COLON && free buffer
    delete[] pToConvert;
    if(pRes != (pEnd-1))
        pRes++;
    else
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size limit, end of stream." << endl;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

    pTemp=pRes;

    tileSizeLimit=r_Interest_Tiling::get_tilesize_limit_from_name(pTemp);
    if(tileSizeLimit==r_Interest_Tiling::NUMBER)
    {
        RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << encoded << "): Error decoding tile size limit from \"" << pTemp << "\"." << endl;
        throw r_Error(TILINGPARAMETERNOTCORRECT);
    }

    iareas = vectInterestAreas;
    ts_strat = tileSizeLimit;
    dimension = tileD;
    tile_size = tileS;
}

r_Interest_Tiling::r_Interest_Tiling(r_Dimension dim, const std::vector<r_Minterval>& interest_areas, r_Bytes ts, Tilesize_Limit strat) throw (r_Error)
    :   r_Dimension_Tiling(dim, ts),
        iareas(interest_areas),
        ts_strat(strat)
{
    for (std::vector<r_Minterval>::iterator it = iareas.begin(); it != iareas.end(); it++)
        if (it->dimension() != dimension)
        {
            RMInit::logOut << "r_Interest_Tiling::r_Interest_Tiling(" << dim << ", " << interest_areas << ", " << ts << ", " << strat << ") the interest area domain " << *it << " does not match the dimension of this tiling scheme (" << dimension << ")" << endl;
            throw r_Edim_mismatch(dimension, it->dimension());
        }
}

r_Interest_Tiling::~r_Interest_Tiling()
{
}

r_Tiling* r_Interest_Tiling::clone() const
{
    r_Tiling* copy = new r_Interest_Tiling(dimension, iareas, tile_size, ts_strat);
    return copy;
}

void r_Interest_Tiling::print_status(std::ostream& os) const
{
    os << "r_Interest_Tiling[ ";
    r_Dimension_Tiling::print_status(os);
    os << " interest areas = " << iareas << ", tiling strategy = " << ts_strat << " ]";
}

r_Tiling_Scheme
r_Interest_Tiling::get_tiling_scheme() const
{
    return r_InterestTiling;
}

static int r_Range_comp(const void *elem1, const void *elem2)
{
    r_Range e1 = *((r_Range*) elem1);
    r_Range e2 = *((r_Range*) elem2);

    if (e1 == e2)
        return 0;

    if (e1 < e2)
        return -1;
    else
        return +1;
}

std::vector<r_Dir_Decompose>*
r_Interest_Tiling::make_partition(const r_Minterval& domain) const
{
    r_Dimension dim = domain.dimension();
    int total = 2 * iareas.size();

    // We need one decomp from each dimension
    std::vector<r_Dir_Decompose>* part = new std::vector<r_Dir_Decompose>(dim);

    // We have at most (number of interest areas + 2) intervals
    r_Range* intervals = new r_Range[total + 2];

    // Create iterator for interest areas
    std::vector<r_Minterval>::const_iterator it = iareas.begin();
    int total_iareas = iareas.size();


    // For all dimensions
    for (r_Dimension i = 0; i < dim; i++)
    {
        it = iareas.begin();                              // Reset iterator
        intervals[0] = domain[i].low();               // Input lower domain limit
        intervals[total+1] = domain[i].high();        // Input higher domain limit


        for (int j = 1; j < total + 1; j += 2, ++it)          // For all possible intervals
        {
            if ((*it)[i].low()-1 <= domain[i].low())    // Input low iarea limit
                intervals[j] = domain[i].low();
            else
                intervals[j] = (*it)[i].low()-1;

            intervals[j+1] = (*it)[i].high();           // Input higher iarea limit
        }

        // Sort the table
        qsort((void*) intervals, total+2, sizeof(r_Range), r_Range_comp);

        // Create partition using the limits table
        for (int k=0; k<total+2; k++)               // all limits must be checked
        {
            if (k == total+1)                         // if on the last limit...
                ((*part)[i]) << intervals[k];                //   input it
            else                                      // else
                if (intervals[k] != intervals[k+1])     //   if it is unique
                    ((*part)[i]) << intervals[k];              //     input it
        }
    }

    // Free memory
    delete [] intervals;

    // Return result
    return part;
}

std::vector<r_Minterval>*
r_Interest_Tiling::group(std::vector<r_Minterval>& blocks, r_Bytes typelen, Blocks_Type btype) const
{
    r_Bytes tilesize = get_tile_size();
    int joins = 0;
    bool group_blocks = true;

    // The list of threated blocks
    std::vector<r_Minterval>* treated = new std::vector<r_Minterval>;

    // An iterator for the blocks list
    std::vector<r_Minterval>::iterator blocks_it = blocks.begin();

    // The current block
    r_Minterval current_block;

    // For all the blocks in list
    while (!blocks.empty())
    {
        // Get first block from list
        current_block = blocks.back();
        blocks.pop_back();

        //this is neccessary when the compiler optimizes out the .end() check
        unsigned int numberOfLevels = blocks.size();
        for (blocks_it = blocks.begin(); blocks_it != blocks.end(); blocks_it++)
        {
            if (numberOfLevels == 0)
            {
                RMInit::logOut << "r_Interest_Tiling::group() the for loop was incorrectly optimized.  breaking the loop." << endl;
                break;
            }
            r_Minterval aux = *blocks_it;

            // In principle two blocks can't be merged
            group_blocks = false;

            // If they can be merged
            if (current_block.is_mergeable(aux))
            {
                std::vector<r_Minterval>::iterator ia_it = blocks.begin();

                switch (btype)
                {
                case BLOCKS_A:

                    group_blocks = true;

                    // Check if the two blocks belong exaclty to the same iareas
                    for (; ia_it != blocks.end(); ia_it++)
                    {
                        if (aux.intersects_with(*ia_it) != current_block.intersects_with(*ia_it))
                        {
                            group_blocks = false;
                            break;
                        }
                    }

                    break;

                case BLOCKS_B:

                    for (; ia_it != blocks.end(); ia_it++) // For all iareas
                    {
                        // Find the one this block intersects
                        if (current_block.intersects_with(*ia_it))
                        {
                            // Check if the other area intersects it
                            if (!aux.intersects_with(*ia_it))
                                group_blocks = false;
                            else
                                group_blocks = true;

                            break;
                        }
                    }

                    if (!group_blocks)
                        break;

                    group_blocks = false;

                case BLOCKS_C: // Falls in (this is, also applies to B);

                    // Only on this two strategies, tilesize should be looked at
                    if ((ts_strat == REGROUP) || (ts_strat == REGROUP_AND_SUBTILING))
                    {
                        // If the resulting size isn't larger than tilesize
                        if ((current_block.cell_count()+aux.cell_count()) * typelen
                                < get_tile_size())
                            group_blocks = true;
                    }
                    else
                        group_blocks = true;

                    break;
                }
            }

            // take care of the iterator advance, if is possible
            if (group_blocks)
            {
                // Merge them
                blocks.erase(blocks_it); // Automatically advances to next block
                // take care of the size of the blocks
                numberOfLevels--;
                current_block.closure_with(aux);
                ++joins;
            }
            else
                numberOfLevels--;
        }

        // Update the treated list with the current block
        treated->push_back(current_block);
    }

    std::vector<r_Minterval>* result = 0;

    // If there was joins, the algoritm must be repeted
    if (joins != 0)
    {
        result = group(*treated, typelen, btype);
        delete treated;
    }
    else
    {
        result = treated;
    }

    // Return the result
    return result;
}

std::vector<r_Minterval>*
r_Interest_Tiling::compute_tiles(const r_Minterval& domain, r_Bytes typelen) const throw (r_Error)
{
    r_Dimension num_dims = domain.dimension();                   // Dimensionality of dom
    if (domain.dimension() != dimension)
    {
        RMInit::logOut << "r_Interest_Tiling::compute_tiles(" << domain << ", " << typelen << ") dimension (" << dimension << ") does not match dimension of object to tile (" << num_dims << ")" << endl;
        throw r_Edim_mismatch(dimension, num_dims);
    }
    if (typelen > tile_size)
    {
        RMInit::logOut << "r_Interest_Tiling::compute_tiles(" << domain << ", " << typelen << ") tile size (" << tile_size << ") is smaller than type length (" << typelen << ")" << endl;
        throw r_Error(TILESIZETOOSMALL);
    }

#ifdef _VISUALIZE_2D_DECOMP_                           // User wants a visual 
    static int count;                                    // Number of decomps
    ++count;                                             // Update num decomps
    // of the 2D decomp.
    Visual_Tiling_2D* vis;
    if (domain.dimension() == 2)
    {
        // Create an object for visualization
        char fname[80];
        sprintf(fname, "2D_decomp_int_%d.ppm", count);
        vis = new Visual_Tiling_2D(domain, fname);
    }

#endif


    // *** Main algoritm ***

    // The result
    std::vector<r_Minterval>* result = new std::vector<r_Minterval>;

    // Create a partition for dir tiling
    std::vector<r_Dir_Decompose>* part = make_partition(domain);

    // Perform dirtiling
    r_Dir_Tiling *dir_tiling= NULL;
    std::vector<r_Minterval>* dir_domain=NULL;

    try
    {
        dir_tiling=new r_Dir_Tiling(num_dims, *part, tile_size, r_Dir_Tiling::WITHOUT_SUBTILING);

        // Compute all tiles (directional tiles)
        dir_domain = dir_tiling->compute_tiles(domain, typelen);

        delete dir_tiling;
    }
    catch(r_Error& err)
    {
        delete result;
        delete part;
        if(dir_tiling)
            delete dir_tiling;
        throw;
    }

    // Get an interator for interest areas
    //std::vector<r_Minterval>::iterator interest_area;

    // Create a list for holding the classifed blocks
    std::vector<Classified_Block> part_domain;

    // Finds how many intersections exist between a block an the interest areas
    for (std::vector<r_Minterval>::iterator dir_block = dir_domain->begin(); dir_block != dir_domain->end(); dir_block++)
    {
        Classified_Block b(*dir_block, 0);

        for (std::vector<r_Minterval>::const_iterator interest_area = iareas.begin(); interest_area != iareas.end(); interest_area++)
        {
            if (b.block.intersects_with(*interest_area))
                ++b.intersection_count;
        }

        part_domain.push_back(b);
    }

    // Lists used for grouping blocks
    std::vector<r_Minterval> Out;
    std::vector<r_Minterval> In_Unique;
    std::vector<r_Minterval> In_Common;

    // Divide blocks into lists according to their number of intersections
    for (std::vector<Classified_Block>::iterator class_block = part_domain.begin(); class_block != part_domain.end(); class_block++)
    {
        switch ((*class_block).intersection_count)
        {
        case 0:
            Out.push_back((*class_block).block);
            break;
        case 1:
            In_Unique.push_back((*class_block).block);
            break;
        default:
            In_Common.push_back((*class_block).block);
            break;
        }
    }


    // Group blocks
    std::vector<r_Minterval>* Blocks_A = group(In_Common, typelen, BLOCKS_A);
    std::vector<r_Minterval>* Blocks_B = group(In_Unique, typelen, BLOCKS_B);
    std::vector<r_Minterval>* Blocks_C = group(Out, typelen, BLOCKS_C);


    std::vector<r_Minterval>::iterator it_A = Blocks_A->begin();
    std::vector<r_Minterval>::iterator it_B = Blocks_B->begin();
    std::vector<r_Minterval>::iterator it_C = Blocks_C->begin();

    // If may be necessary to perform sub-tiling
    if ((ts_strat == SUB_TILING) || (ts_strat == REGROUP_AND_SUBTILING))
    {
        // Variable to hold result of sub-tiling
        std::vector<r_Minterval>* subtiles = 0;

        // We need an array to hold the 3 iterators of the resulting blocks
        std::vector<r_Minterval>::iterator* blocks_it[3];
        std::vector<r_Minterval>* blocks_vec[3];

        // Put iterators on the list
        blocks_it[0] = &it_A;
        blocks_it[1] = &it_B;
        blocks_it[2] = &it_C;

        blocks_vec[0] = Blocks_A;
        blocks_vec[1] = Blocks_B;
        blocks_vec[2] = Blocks_C;

        // For all the lists (Blocs_A, Blocks_B and Blocks_C)
        for (int j=0; j<3; j++)
        {
            // Reset iterator to the begin
            (*blocks_it[j]) = blocks_vec[j]->begin();

            // Tile each block if necessary
            for (; *blocks_it[j] != blocks_vec[j]->end(); (*blocks_it[j])++)
            {
                if ((**blocks_it[j]).cell_count()*typelen > get_tile_size())
                {
                    // Create a specification of a regular n-dim cube grid
                    r_Minterval specs(num_dims);
                    for (r_Dimension i=0; i<num_dims; i++)
                    {
                        specs << r_Sinterval((r_Range)0, (**blocks_it[j])[i].high() - (**blocks_it[j])[i].low());
                    }

                    // Class for performing sub-tiling
                    r_Aligned_Tiling subtiling(specs, get_tile_size());

                    subtiles = subtiling.compute_tiles(**blocks_it[j], typelen);
                    for (std::vector<r_Minterval>::iterator it = subtiles->begin(); it != subtiles->end(); it++)
                    {
                        result->push_back(*it);

#ifdef _VISUALIZE_2D_DECOMP_
                        if (domain.dimension() == 2)
                            (*vis) << (*it);
#endif
                    }

                    delete subtiles;
                }
                else    // No subtiling needed
                {
                    // Insert block as it is
                    result->push_back(**blocks_it[j]);

#ifdef _VISUALIZE_2D_DECOMP_
                    if (domain.dimension() == 2)
                        (*vis) << (**blocks_it[j]);
#endif
                }
            }
        }

#ifdef _VISUALIZE_2D_DECOMP_
        if (domain.dimension() == 2)
        {
            vis->set_pen(255, 255, 0);

            interest_area.seek_begin();
            while (interest_area.not_done())
            {
                (*vis) << (*interest_area);
                ++interest_area;
            }
        }
#endif

    }
    else
    {
        // The result is just the sum of all blocks

        result->insert(result->end(), Blocks_A->begin(), Blocks_A->end());
        result->insert(result->end(), Blocks_B->begin(), Blocks_B->end());
        result->insert(result->end(), Blocks_C->begin(), Blocks_C->end());

        // Visualization

#ifdef _VISUALIZE_2D_DECOMP_
        if (domain.dimension() == 2)
        {
            vis->set_pen(255, 255, 0);
            for (it_A.seek_begin(); it_A.not_done(); it_A++)
                (*vis) << *it_A;

            vis->set_pen(0, 255, 0);
            for (it_B.seek_begin(); it_B.not_done(); it_B++)
                (*vis) << *it_B;

            vis->set_pen(0, 0, 255);
            for (it_C.seek_begin(); it_C.not_done(); it_C++)
                (*vis) << *it_C;
        }
#endif
    }

    // *** End of the main algorithm ***

    // *** Clean up ***

    delete part;
    delete dir_domain;
    delete Blocks_A;
    delete Blocks_B;
    delete Blocks_C;

#ifdef _VISUALIZE_2D_DECOMP_
    if (domain.dimension() == 2)
        delete vis;
#endif

    // Return result

    return result;
}



std::ostream& operator<<(std::ostream& os, const Classified_Block block);


Classified_Block::Classified_Block(int count)
    :   intersection_count(count)
{
}

Classified_Block::Classified_Block(const r_Minterval& b, int count)
    :   intersection_count(count),
        block(b)
{
}


bool Classified_Block::operator==(const Classified_Block& other) const
{
    return (this->block == other.block)
           && (this->intersection_count == other.intersection_count);
}

bool Classified_Block::operator!=(const Classified_Block& other) const
{
    return !(*this == other);
}


std::ostream& operator<<(std::ostream& os, const Classified_Block block)
{
    os << "CBlock(" << block.intersection_count << "x: " << block.block << ")";

    return os;
}