Commit 3e38dc0c9afbb4e34287b2a71d47b468135ffbee

Authored by Pierre Lassalle
1 parent 3f4c3298
Exists in master

New version of the GRM library

Applications/CMakeLists.txt 0 → 100644
... ... @@ -0,0 +1,19 @@
  1 +#=========================================================================
  2 +
  3 +# Program: Generic Region Merging Library (GRM)
  4 +# Language: C++
  5 +# author: Lassalle Pierre
  6 +
  7 +
  8 +
  9 +# Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
  10 +
  11 +# See grmlib-copyright.txt for details.
  12 +
  13 +# This software is distributed WITHOUT ANY WARRANTY; without even
  14 +# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15 +# PURPOSE. See the above copyright notices for more information.
  16 +
  17 +#=========================================================================
  18 +add_executable(RegionMergingSegmentation RegionMergingSegmentation.cxx)
  19 +target_link_libraries(RegionMergingSegmentation OTBGRM)
... ...
Applications/RegionMergingSegmentation.cxx 0 → 100644
... ... @@ -0,0 +1,52 @@
  1 +#include <iostream>
  2 +#include <otbImage.h>
  3 +#include <otbImageFileReader.h>
  4 +#include <otbVectorImage.h>
  5 +#include "lsrmBaatzSegmenter.h"
  6 +
  7 +int main(int argc, char *argv[])
  8 +{
  9 + if(argc != 7)
  10 + {
  11 + std::cerr << "Usage: ./" << argv[0] << "\n"
  12 + << "\t[input image path] : (.jpg, .png, .tif)\n"
  13 + << "\t[output clustered image] : (.jpg, .png, .tif)\n"
  14 + << "\t[output label image] : (.tif)\n"
  15 + << "\t[spectral weight] : range from 0 to 1\n"
  16 + << "\t[shape weight] : range from 0 to 1\n"
  17 + << "\t[scale threshold] : unlimited positive value"
  18 + << std::endl;
  19 + return 1;
  20 + }
  21 +
  22 + lsrm::BaatzParam params;
  23 + const char * input_image = argv[1];
  24 + const char * clustered_image = argv[2];
  25 + const char * label_image = argv[3];
  26 + params.m_SpectralWeight = atof(argv[4]);
  27 + params.m_ShapeWeight = atof(argv[5]);
  28 + float sqrt_scale = atof(argv[6]);
  29 + const float scale = sqrt_scale * sqrt_scale;
  30 +
  31 + typedef float PixelType;
  32 + typedef otb::VectorImage<PixelType, 2> ImageType;
  33 + typedef lsrm::BaatzSegmenter<ImageType> SegmenterType;
  34 +
  35 + SegmenterType segmenter;
  36 + segmenter.SetParam(params);
  37 + segmenter.SetThreshold(scale);
  38 + segmenter.SetInputFileName(input_image);
  39 + segmenter.SetClusteredImageFileName(clustered_image);
  40 + segmenter.SetLabelImageFileName(label_image);
  41 +
  42 + segmenter.RunSegmentation();
  43 +
  44 + return 0;
  45 +}
  46 +
  47 +
  48 +
  49 +
  50 +
  51 +
  52 +
... ...
CMakeLists.txt 0 → 100644
... ... @@ -0,0 +1,33 @@
  1 +#=========================================================================
  2 +
  3 +# Program: Generic Region Merging Library (GRM)
  4 +# Language: C++
  5 +# author: Lassalle Pierre
  6 +
  7 +
  8 +
  9 +# Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
  10 +
  11 +# See grmlib-copyright.txt for details.
  12 +
  13 +# This software is distributed WITHOUT ANY WARRANTY; without even
  14 +# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15 +# PURPOSE. See the above copyright notices for more information.
  16 +
  17 +#=========================================================================
  18 +project(GRM)
  19 +
  20 +cmake_minimum_required(VERSION 2.8)
  21 +
  22 +find_package(OTB)
  23 +IF(OTB_FOUND)
  24 + include(${OTB_USE_FILE})
  25 +ELSE(OTB_FOUND)
  26 + message(FATAL_ERROR
  27 + "Cannot build OTB project without OTB. Please set OTB_DIR.")
  28 +ENDIF(OTB_FOUND)
  29 +
  30 +include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Code)
  31 +
  32 +add_subdirectory(Code)
  33 +add_subdirectory(Applications)
... ...
Code/CMakeLists.txt 0 → 100644
... ... @@ -0,0 +1,22 @@
  1 +#=========================================================================
  2 +
  3 +# Program: Generic Region Merging Library (GRM)
  4 +# Language: C++
  5 +# author: Lassalle Pierre
  6 +
  7 +
  8 +
  9 +# Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
  10 +
  11 +# See grmlib-copyright.txt for details.
  12 +
  13 +# This software is distributed WITHOUT ANY WARRANTY; without even
  14 +# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15 +# PURPOSE. See the above copyright notices for more information.
  16 +
  17 +#=========================================================================
  18 +add_library(OTBGRM
  19 + lsrmContourOperations.cpp
  20 + lsrmNeighborhood.cpp)
  21 +
  22 +target_link_libraries(OTBGRM OTBCommon OTBIO)
... ...
Code/lsrmBaatzSegmenter.h 0 → 100644
... ... @@ -0,0 +1,102 @@
  1 +#ifndef __LSRM_BAATZ_SEGMENTER_H
  2 +#define __LSRM_BAATZ_SEGMENTER_H
  3 +#include "lsrmSegmenter.h"
  4 +#include "lsrmGraphToOtbImage.h"
  5 +/*
  6 + * Tutorial : Implementation of the Baatz & Schape criterion
  7 + *
  8 + * Details about the criterion can be found in the following publication:
  9 + *
  10 + * Martin Baatz and Arno Schape. Multiresolution segmentation: an optimization approach for high quality multi-scale image segmentation.
  11 + * Angewandte Geographische Informationsverarbeitung XII, pages 12–23, 2000.
  12 + *
  13 + * The steps are ordered has to be followed in a chronological way for a full understanding.
  14 + * This tutorial does not aim at explaining all the details of the large scale region merging
  15 + * but gives all the elements to use it properly.
  16 + */
  17 +
  18 +namespace lsrm
  19 +{
  20 + /*
  21 + * Step 1 :
  22 + * We define the specific attributes required for the Baatz & Schape criterion.
  23 + * Regions are represented by nodes in a graph and the lsrm library has an internal
  24 + * representation of the nodes with the class lsrm::Node.
  25 + * --> A node contains a unique number to be indentified, it corresponds to the vectorized
  26 + * coordinates of its first pixel and its name is m_Id.
  27 + * To retrieve the 2D coordinates from the value of m_Id, it is necessary to know the dimension
  28 + * of the image (the number of rows and columns), then :
  29 + * x = m_Id % columns and y = m_Id / columns.
  30 + *
  31 + * --> A node contains the perimeter (m_Perimeter) of its corresponding region, it is the length of the external
  32 + * boundary of the region. For example a region of one pixel has a perimeter of 4.
  33 + *
  34 + * --> A node contains the area (m_Area) of its corresponding region, it is the number of pixels within the
  35 + * region.
  36 + *
  37 + * --> A node contains the minimum rectangular bounding box (parrallel to the image axis) of the region (m_Bbox).
  38 + * the boundix box is determined by the structure lsrm::BoundingBox with the coordinates of the upper left
  39 + * corner pixel and the width and the height (m_UX, m_UY, m_W, m_H).
  40 + *
  41 + * After Reading the paper about the Baatz & Schape criterion, we can conclude that we have to define additional
  42 + * spectral attributes for a node: its mean vector, the square mean vector, the sum of all the pixels
  43 + * and the standard deviation vector.
  44 + *
  45 + * We define a new class BaatzNode which inherits from the template structure Node. Notice that the template
  46 + * is the derived node type, then the template type is BaatzNode.
  47 + */
  48 + struct BaatzNode : Node<BaatzNode>
  49 + {
  50 + std::vector<float> m_Means;
  51 + std::vector<float> m_SquareMeans;
  52 + std::vector<float> m_SpectralSum;
  53 + std::vector<float> m_Std;
  54 + };
  55 +
  56 + /*
  57 + * Step 2
  58 + *
  59 + * The Baatz & Schape criterion has user-defined parameters, the spectral weight
  60 + * which is the relative importance given to the spectral components and the shape
  61 + * weight for the shape components.
  62 + *
  63 + * We define then a structure wrapping these two parameters we called BaatzParam.
  64 + */
  65 + struct BaatzParam
  66 + {
  67 + float m_SpectralWeight;
  68 + float m_ShapeWeight;
  69 + };
  70 +
  71 + /*
  72 + * Step 3 :
  73 + *
  74 + */
  75 +
  76 + template<class TImage>
  77 + class BaatzSegmenter : public Segmenter< TImage, BaatzNode, BaatzParam>
  78 + {
  79 + public:
  80 +
  81 + /* Some convenient typedefs */
  82 + typedef Segmenter<TImage, BaatzNode, BaatzParam> Superclass;
  83 + typedef TImage ImageType;
  84 + typedef typename Superclass::GraphType GraphType;
  85 + typedef typename Superclass::NodePointerType NodePointerType;
  86 + typedef typename Superclass::GraphOperatorType GraphOperatorType;
  87 + typedef GraphToOtbImage<GraphType> IOType;
  88 +
  89 + BaatzSegmenter();
  90 +
  91 + void RunSegmentation();
  92 + float ComputeMergingCost(NodePointerType n1, NodePointerType n2);
  93 + void UpdateSpecificAttributes(NodePointerType n1, NodePointerType n2);
  94 + void InitFromImage();
  95 + };
  96 +
  97 +} // end of namespace lsrm
  98 +#include "lsrmBaatzSegmenter.txx"
  99 +#endif
  100 +
  101 +
  102 +
... ...
Code/lsrmBaatzSegmenter.txx 0 → 100644
... ... @@ -0,0 +1,159 @@
  1 +#ifndef __LSRM_BAATZ_SEGMENTER_TXX
  2 +#define __LSRM_BAATZ_SEGMENTER_TXX
  3 +#include <otbImageFileReader.h>
  4 +#include <itkImageRegionIterator.h>
  5 +
  6 +namespace lsrm
  7 +{
  8 +
  9 + template<class TImage>
  10 + BaatzSegmenter<TImage>::BaatzSegmenter() : Superclass()
  11 + {
  12 + this->m_DoBFSegmentation = true;
  13 + this->m_NumberOfIterations = 75;
  14 + }
  15 +
  16 + template<class TImage>
  17 + void
  18 + BaatzSegmenter<TImage>::InitFromImage()
  19 + {
  20 + typedef otb::ImageFileReader<TImage> ImageReaderType;
  21 + typedef itk::ImageRegionIterator<TImage> ImageIterator;
  22 +
  23 + typename ImageReaderType::Pointer reader = ImageReaderType::New();
  24 + reader->SetFileName(this->m_InputFileName);
  25 + reader->Update();
  26 +
  27 + this->m_ImageWidth = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
  28 + this->m_ImageHeight = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[1];
  29 + this->m_NumberOfComponentsPerPixel = reader->GetOutput()->GetNumberOfComponentsPerPixel();
  30 +
  31 + std::size_t idx = 0;
  32 + ImageIterator it(reader->GetOutput(), reader->GetOutput()->GetLargestPossibleRegion());
  33 + for(it.GoToBegin(); !it.IsAtEnd(); ++it)
  34 + {
  35 + this->m_Graph.m_Nodes[idx]->m_Means.reserve(this->m_NumberOfComponentsPerPixel);
  36 + this->m_Graph.m_Nodes[idx]->m_SquareMeans.reserve(this->m_NumberOfComponentsPerPixel);
  37 + this->m_Graph.m_Nodes[idx]->m_SpectralSum.reserve(this->m_NumberOfComponentsPerPixel);
  38 + this->m_Graph.m_Nodes[idx]->m_Std.assign(this->m_NumberOfComponentsPerPixel, 0.0f);
  39 +
  40 + for(std::size_t b = 0; b < this->m_NumberOfComponentsPerPixel; ++b)
  41 + {
  42 + this->m_Graph.m_Nodes[idx]->m_Means.push_back(it.Get()[b]);
  43 + this->m_Graph.m_Nodes[idx]->m_SquareMeans.push_back((it.Get()[b])*(it.Get()[b]));
  44 + this->m_Graph.m_Nodes[idx]->m_SpectralSum.push_back(it.Get()[b]);
  45 + }
  46 + ++idx;
  47 + }
  48 + }
  49 +
  50 + template<class TImage>
  51 + float
  52 + BaatzSegmenter<TImage>::ComputeMergingCost(NodePointerType n1, NodePointerType n2)
  53 + {
  54 + const std::size_t bands = n1->m_Means.size();
  55 + const unsigned int a1 = n1->m_Area, a2 = n2->m_Area, a_sum = a1 + a2;
  56 +
  57 + float spect_cost = 0.0f;
  58 + float mean, square_mean, sum, std;
  59 +
  60 + for (unsigned int b = 0; b < this->m_NumberOfComponentsPerPixel; b++)
  61 + {
  62 + mean = ((a1 * n1->m_Means[b]) + (a2 * n2->m_Means[b])) / a_sum;
  63 + square_mean = n1->m_SquareMeans[b] + n2->m_SquareMeans[b];
  64 + sum = n1->m_SpectralSum[b] + n2->m_SpectralSum[b];
  65 + std = std::sqrt((square_mean - 2*mean*sum + a_sum * mean* mean) / a_sum);
  66 + spect_cost += (a_sum * std - a1 * n1->m_Std[b] - a2 * n2->m_Std[b]);
  67 + }
  68 + spect_cost *= this->m_Param.m_ShapeWeight;
  69 +
  70 + if(spect_cost < this->m_Threshold)
  71 + {
  72 + float shape_cost, smooth_f, compact_f;
  73 +
  74 + // Compute the shape merging cost
  75 + const float p1 = static_cast<float>(n1->m_Perimeter);
  76 + const float p2 = static_cast<float>(n2->m_Perimeter);
  77 + const unsigned int boundary = (GraphOperatorType::FindEdge(n1, n2))->m_Boundary;
  78 + const float p3 = p1 + p2 - 2 * static_cast<float>(boundary);
  79 +
  80 + const BoundingBox merged_bbox = ContourOperations::MergeBoundingBoxes(n1->m_Bbox, n2->m_Bbox);
  81 + const float bb1_perimeter = static_cast<float>(2*n1->m_Bbox.m_W + 2*n1->m_Bbox.m_H);
  82 + const float bb2_perimeter = static_cast<float>(2*n2->m_Bbox.m_W + 2*n2->m_Bbox.m_H);
  83 + const float mbb_perimeter = static_cast<float>(2 * merged_bbox.m_W + 2 * merged_bbox.m_H);
  84 +
  85 + smooth_f = a_sum*p3/mbb_perimeter - a1*p1/bb1_perimeter - a2*p2/bb2_perimeter;
  86 + compact_f = a_sum*p3/std::sqrt(a_sum) - a1*p1/std::sqrt(a1) - a2*p2/std::sqrt(a2);
  87 +
  88 + shape_cost = this->m_Param.m_ShapeWeight * compact_f + (1-this->m_Param.m_ShapeWeight) * smooth_f;
  89 +
  90 + return (spect_cost + (1-this->m_Param.m_ShapeWeight)*shape_cost);
  91 + }
  92 + else
  93 + return spect_cost;
  94 + }
  95 +
  96 + template<class TImage>
  97 + void
  98 + BaatzSegmenter<TImage>::UpdateSpecificAttributes(NodePointerType n1, NodePointerType n2)
  99 + {
  100 + const float a1 = static_cast<float>(n1->m_Area);
  101 + const float a2 = static_cast<float>(n2->m_Area);
  102 + const float a_sum = a1 + a2;
  103 +
  104 + for(unsigned int b = 0; b < this->m_NumberOfComponentsPerPixel; ++b)
  105 + {
  106 + n1->m_Means[b] = (a1 * n1->m_Means[b] + a2 * n2->m_Means[b]) / a_sum;
  107 + n1->m_SquareMeans[b] += n2->m_SquareMeans[b];
  108 + n1->m_SpectralSum[b] += n2->m_SpectralSum[b];
  109 + n1->m_Std[b] = std::sqrt((n1->m_SquareMeans[b] - 2 * n1->m_Means[b] * n1->m_SpectralSum[b] +
  110 + a_sum * n1->m_Means[b] * n1->m_Means[b]) / a_sum);
  111 + }
  112 + }
  113 +
  114 + template<class TImage>
  115 + void
  116 + BaatzSegmenter<TImage>::RunSegmentation()
  117 + {
  118 + GraphOperatorType::InitNodes(this->m_Graph, *this, this->m_InputFileName, FOUR);
  119 +
  120 + bool prev_merged =
  121 + GraphOperatorType::PerfomAllIterationsWithLMBFAndConstThreshold(this->m_Graph, *this,
  122 + this->m_Threshold, this->m_NumberOfIterations,
  123 + this->m_ImageWidth, this->m_ImageHeight);
  124 +
  125 + if(prev_merged && this->m_DoBFSegmentation)
  126 + {
  127 + prev_merged =
  128 + GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(this->m_Graph, *this,
  129 + this->m_Threshold, this->m_NumberOfIterations,
  130 + this->m_ImageWidth, this->m_ImageHeight);
  131 + }
  132 +
  133 + if(!this->m_ClusteredImageFileName.empty())
  134 + IOType::WriteOutputRGBImage(this->m_Graph,
  135 + this->m_ImageWidth,
  136 + this->m_ImageHeight,
  137 + this->m_ClusteredImageFileName);
  138 +
  139 + if(!this->m_LabelImageFileName.empty())
  140 + IOType::WriteLabelImage(this->m_Graph,
  141 + this->m_ImageWidth,
  142 + this->m_ImageHeight,
  143 + this->m_LabelImageFileName);
  144 +
  145 + if(!this->m_ContourImageFileName.empty())
  146 + IOType::WriteContourImage(this->m_Graph,
  147 + this->m_InputFileName,
  148 + this->m_LabelImageFileName);
  149 + }
  150 +} // end of namespace lsrm
  151 +
  152 +#endif
  153 +
  154 +
  155 +
  156 +
  157 +
  158 +
  159 +
... ...
Code/lsrmContourOperations.cpp 0 → 100644
... ... @@ -0,0 +1,689 @@
  1 +#include "lsrmContourOperations.h"
  2 +
  3 +namespace lsrm
  4 +{
  5 + BoundingBox
  6 + ContourOperations::MergeBoundingBoxes(const BoundingBox& bb1,
  7 + const BoundingBox& bb2)
  8 + {
  9 + long unsigned int min_ux, min_uy, max_xw, max_yh;
  10 + BoundingBox bb;
  11 +
  12 + min_ux = std::min(bb1.m_UX, bb2.m_UX);
  13 + min_uy = std::min(bb1.m_UY, bb2.m_UY);
  14 + max_xw = std::max(bb1.m_UX + bb1.m_W, bb2.m_UX + bb2.m_W);
  15 + max_yh = std::max(bb1.m_UY + bb1.m_H, bb2.m_UY + bb2.m_H);
  16 +
  17 + bb.m_UX = min_ux;
  18 + bb.m_UY = min_uy;
  19 + bb.m_W = max_xw - min_ux;
  20 + bb.m_H = max_yh - min_uy;
  21 +
  22 + return bb;
  23 + }
  24 +
  25 + long unsigned int
  26 + ContourOperations::ToBoundingBoxCoordinates(const long unsigned int id,
  27 + const BoundingBox& bb,
  28 + const unsigned int width)
  29 + {
  30 + long unsigned int img_x = id % width;
  31 + long unsigned int img_y = id / width;
  32 + long unsigned int bb_x = img_x - bb.m_UX;
  33 + long unsigned int bb_y = img_y - bb.m_UY;
  34 + return (bb_y * bb.m_W + bb_x);
  35 + }
  36 +
  37 + typename ContourOperations::PixelList
  38 + ContourOperations::GenerateBorderPixels(long unsigned int id,
  39 + const Contour& contour,
  40 + const unsigned int width)
  41 + {
  42 + PixelList pixels;
  43 + pixels.push_back(id);
  44 +
  45 + short curr_move, prev_move;
  46 +
  47 + if(contour.size() > 4)
  48 + {
  49 + for(ContourConstIterator m = contour.begin()+1; m != contour.end(); m++)
  50 + {
  51 + curr_move = (*m)[0] + 2*(*m)[1];
  52 + prev_move = (*(m-1))[0] + 2*(*(m-1))[1];
  53 +
  54 + if(curr_move == 0)
  55 + {
  56 + if(prev_move == 0)
  57 + {
  58 + id-=width;
  59 + pixels.push_back(id);
  60 + }
  61 + else if(prev_move == 1) // 1
  62 + {
  63 + id++;
  64 + pixels.push_back(id);
  65 + id-=width;
  66 + pixels.push_back(id);
  67 + }
  68 + }
  69 + else if(curr_move == 1)
  70 + {
  71 + if(prev_move == 1)
  72 + {
  73 + id++;
  74 + pixels.push_back(id);
  75 + }
  76 + else if(prev_move == 2)
  77 + {
  78 + id+=width;
  79 + pixels.push_back(id);
  80 + id++;
  81 + pixels.push_back(id);
  82 + }
  83 + }
  84 + else if(curr_move == 2)
  85 + {
  86 + if(prev_move == 3)
  87 + {
  88 + id--;
  89 + pixels.push_back(id);
  90 + id+=width;
  91 + pixels.push_back(id);
  92 + }
  93 + else if(prev_move == 2)
  94 + {
  95 + id+=width;
  96 + pixels.push_back(id);
  97 + }
  98 + }
  99 + else if(curr_move == 3)
  100 + {
  101 + if(prev_move==0)
  102 + {
  103 + id-=width;
  104 + pixels.push_back(id);
  105 + id--;
  106 + pixels.push_back(id);
  107 + }
  108 + else if(prev_move == 3)
  109 + {
  110 + id--;
  111 + pixels.push_back(id);
  112 + }
  113 + }
  114 + }
  115 + // Remove duplicated pixels
  116 + std::sort(pixels.begin(), pixels.end());
  117 + PixelIterator it = std::unique(pixels.begin(), pixels.end());
  118 + pixels.resize(std::distance(pixels.begin(), it));
  119 + }
  120 + return pixels;
  121 + }
  122 +
  123 + bool
  124 + ContourOperations::GetCollisionAtNorth(const long unsigned int idx,
  125 + const short * grid,
  126 + const unsigned int grid_width)
  127 + {
  128 + long int tmp_idx = idx - grid_width;
  129 +
  130 + if(grid[tmp_idx] > 0)
  131 + return true;
  132 +
  133 + while(tmp_idx > 0)
  134 + {
  135 + if(grid[tmp_idx] > 0)
  136 + return true;
  137 + tmp_idx -= grid_width;
  138 + }
  139 + return false;
  140 + }
  141 +
  142 + bool
  143 + ContourOperations::GetCollisionAtNorthEast(const long unsigned int idx,
  144 + const short * grid,
  145 + const unsigned int grid_width)
  146 + {
  147 + long int tmp_idx = idx - grid_width + 1;
  148 +
  149 + if(grid[tmp_idx] > 0)
  150 + return true;
  151 +
  152 + while(tmp_idx > 0 && (tmp_idx % grid_width) != 0)
  153 + {
  154 + if(grid[tmp_idx] > 0)
  155 + return true;
  156 + tmp_idx = tmp_idx - grid_width + 1;
  157 + }
  158 + return false;
  159 + }
  160 +
  161 + bool
  162 + ContourOperations::GetCollisionAtEast(const long unsigned int idx,
  163 + const short * grid,
  164 + const unsigned int grid_width)
  165 + {
  166 + long int tmp_idx = idx + 1;
  167 +
  168 +
  169 + if(grid[tmp_idx] > 0)
  170 + return true;
  171 +
  172 + while((tmp_idx % grid_width) != 0)
  173 + {
  174 + if(grid[tmp_idx] > 0)
  175 + return true;
  176 + ++tmp_idx;
  177 + }
  178 + return false;
  179 + }
  180 +
  181 + bool
  182 + ContourOperations::GetCollisionAtSouthEast(const long unsigned int idx,
  183 + const short * grid,
  184 + const unsigned int grid_width,
  185 + const unsigned int grid_height)
  186 + {
  187 + long int tmp_idx = idx + 1 + grid_width;
  188 +
  189 + if(grid[tmp_idx] > 0)
  190 + return true;
  191 +
  192 + while(tmp_idx < (grid_width * grid_height) &&
  193 + (tmp_idx % grid_width) != 0)
  194 + {
  195 + if(grid[tmp_idx] > 0)
  196 + return true;
  197 + tmp_idx = tmp_idx + 1 + grid_width;
  198 + }
  199 + return false;
  200 + }
  201 +
  202 + bool
  203 + ContourOperations::GetCollisionAtSouth(const long unsigned int idx,
  204 + const short * grid,
  205 + const unsigned int grid_width,
  206 + const unsigned int grid_height)
  207 + {
  208 + long int tmp_idx = idx + grid_width;
  209 +
  210 + if(grid[tmp_idx] > 0)
  211 + return true;
  212 +
  213 + while(tmp_idx < (grid_width * grid_height))
  214 + {
  215 + if(grid[tmp_idx] > 0)
  216 + return true;
  217 + tmp_idx += grid_width;
  218 + }
  219 + return false;
  220 + }
  221 +
  222 + bool
  223 + ContourOperations::GetCollisionAtSouthWest(const long unsigned int idx,
  224 + const short * grid,
  225 + const unsigned int grid_width,
  226 + const unsigned int grid_height)
  227 + {
  228 + long int tmp_idx = idx + grid_width - 1;
  229 +
  230 + if(grid[tmp_idx] > 0)
  231 + return true;
  232 +
  233 + while(tmp_idx < (grid_width * grid_height) &&
  234 + (tmp_idx % grid_width) != grid_width - 1)
  235 + {
  236 + if(grid[tmp_idx] > 0)
  237 + return true;
  238 + tmp_idx = tmp_idx + grid_width - 1;
  239 + }
  240 + return false;
  241 + }
  242 +
  243 + bool
  244 + ContourOperations::GetCollisionAtWest(const long unsigned int idx,
  245 + const short * grid,
  246 + const unsigned int grid_width)
  247 + {
  248 + long int tmp_idx = idx - 1;
  249 +
  250 + if(grid[tmp_idx] > 0)
  251 + return true;
  252 +
  253 + while((tmp_idx % grid_width) != grid_width - 1)
  254 + {
  255 + if(grid[tmp_idx] > 0)
  256 + return true;
  257 + tmp_idx -= 1;
  258 + }
  259 + return false;
  260 + }
  261 +
  262 + bool
  263 + ContourOperations::GetCollisionAtNorthWest(const long unsigned int idx,
  264 + const short * grid,
  265 + const unsigned int grid_width)
  266 + {
  267 + long int tmp_idx = idx - grid_width - 1;
  268 +
  269 + if(grid[tmp_idx] > 0)
  270 + return true;
  271 +
  272 + while(tmp_idx > 0 && (tmp_idx % grid_width) != grid_width - 1)
  273 + {
  274 + if(grid[tmp_idx] > 0)
  275 + return true;
  276 + tmp_idx = tmp_idx - grid_width - 1;
  277 + }
  278 + return false;
  279 + }
  280 +
  281 + bool
  282 + ContourOperations::IsInternalPixel(const long unsigned int idx,
  283 + const short * grid,
  284 + const unsigned int grid_width,
  285 + const unsigned int grid_height)
  286 + {
  287 + if(!GetCollisionAtNorth(idx, grid, grid_width))
  288 + return false;
  289 + if(!GetCollisionAtNorthEast(idx, grid, grid_width))
  290 + return false;
  291 + if(!GetCollisionAtEast(idx, grid, grid_width))
  292 + return false;
  293 + if(!GetCollisionAtSouthEast(idx, grid, grid_width, grid_height))
  294 + return false;
  295 + if(!GetCollisionAtSouth(idx, grid, grid_width, grid_height))
  296 + return false;
  297 + if(!GetCollisionAtSouthWest(idx, grid, grid_width, grid_height))
  298 + return false;
  299 + if(!GetCollisionAtWest(idx, grid, grid_width))
  300 + return false;
  301 + if(!GetCollisionAtNorthWest(idx, grid, grid_width))
  302 + return false;
  303 +
  304 + return true;
  305 + }
  306 +
  307 + void
  308 + ContourOperations::UpdateSight(short * sight, const short direction)
  309 + {
  310 + sight[0] = (direction + 3) % 4; // look at the left
  311 + sight[1] = direction; // look in front
  312 + sight[2] = (direction + 1) % 4; // look at the right
  313 + sight[3] = (direction + 2) % 4; // look behind
  314 + }
  315 +
  316 + void
  317 + ContourOperations::EncodeContourForTopPixel(long unsigned int& curr_mat_id,
  318 + short& direction,
  319 + short& pos,
  320 + Contour& curr_contour,
  321 + const unsigned int grid_width)
  322 + {
  323 + if(pos == 1)
  324 + {
  325 + if(direction==0)
  326 + pos = 4;
  327 + else if(direction == 1)
  328 + {
  329 + pos = 1;
  330 + curr_contour.push_back(0);
  331 + }
  332 + }
  333 + else if(pos == 2)
  334 + {
  335 + pos = 4;
  336 + curr_contour.push_back(2);
  337 + curr_contour.push_back(3);
  338 + curr_contour.push_back(0);
  339 + }
  340 + else if(pos == 3)
  341 + {
  342 + pos = 4;
  343 + curr_contour.push_back(3);
  344 + curr_contour.push_back(0);
  345 + }
  346 + else if(pos == 4)
  347 + curr_contour.push_back(0);
  348 +
  349 +
  350 + curr_mat_id -= grid_width;
  351 + direction = 0;
  352 + }
  353 +
  354 + void
  355 + ContourOperations::EncodeContourForRightPixel(long unsigned int& curr_mat_id,
  356 + short& direction,
  357 + short& pos,
  358 + Contour& curr_contour)
  359 + {
  360 + if(pos == 1)
  361 + curr_contour.push_back(1);
  362 + else if(pos == 2)
  363 + {
  364 + if(direction == 1)
  365 + pos = 1;
  366 + else if(direction == 2)
  367 + curr_contour.push_back(1);
  368 + }
  369 + else if(pos == 3)
  370 + {
  371 + pos = 1;
  372 + curr_contour.push_back(3);
  373 + curr_contour.push_back(0);
  374 + curr_contour.push_back(1);
  375 + }
  376 + else if(pos == 4)
  377 + {
  378 + pos = 1;
  379 + curr_contour.push_back(0);
  380 + curr_contour.push_back(1);
  381 + }
  382 +
  383 + curr_mat_id++;
  384 + direction = 1;
  385 + }
  386 +
  387 + void
  388 + ContourOperations::EncodeContourForBottomPixel(long unsigned int& curr_mat_id,
  389 + short& direction,
  390 + short& pos,
  391 + Contour& curr_contour,
  392 + const unsigned int grid_width)
  393 + {
  394 + if(pos == 1)
  395 + {
  396 + pos = 2;
  397 + curr_contour.push_back(1);
  398 + curr_contour.push_back(2);
  399 + }
  400 + else if(pos == 2)
  401 + curr_contour.push_back(2);
  402 + else if(pos == 3)
  403 + {
  404 + if(direction == 2)
  405 + pos = 2;
  406 + else if(direction == 3)
  407 + curr_contour.push_back(2);
  408 + }
  409 + else if(pos == 4)
  410 + {
  411 + pos = 2;
  412 + curr_contour.push_back(0);
  413 + curr_contour.push_back(1);
  414 + curr_contour.push_back(2);
  415 + }
  416 +
  417 + curr_mat_id += grid_width;
  418 + direction = 2;
  419 + }
  420 +
  421 + void
  422 + ContourOperations::EncodeContourForLeftPixel(long unsigned int& curr_mat_id,
  423 + short& direction,
  424 + short& pos,
  425 + Contour& curr_contour)
  426 + {
  427 + if(pos == 1)
  428 + {
  429 + pos = 3;
  430 + curr_contour.push_back(1);
  431 + curr_contour.push_back(2);
  432 + curr_contour.push_back(3);
  433 + }
  434 + else if(pos == 2)
  435 + {
  436 + pos = 3;
  437 + curr_contour.push_back(2);
  438 + curr_contour.push_back(3);
  439 + }
  440 + else if(pos == 3)
  441 + curr_contour.push_back(3);
  442 + else if(pos == 4)
  443 + {
  444 + if(direction == 3)
  445 + pos = 3;
  446 + else if(direction == 0)
  447 + {
  448 + pos = 4;
  449 + curr_contour.push_back(3);
  450 + }
  451 + }
  452 +
  453 + curr_mat_id --;
  454 + direction = 3;
  455 + }
  456 +
  457 + Contour
  458 + ContourOperations::CreateContourFromBorderPixels(const long unsigned int id,
  459 + const BoundingBox& bbox,
  460 + const short * grid,
  461 + const unsigned int width)
  462 + {
  463 + /* Location on one of the 4 corners of a pixel */
  464 + short direction, pos;
  465 + /* Field of view */
  466 + short sight[4];
  467 +
  468 + long unsigned int start_mat_id = ToBoundingBoxCoordinates(id, bbox, width);
  469 + long unsigned int curr_mat_id = start_mat_id;
  470 +
  471 + /* New contour */
  472 + Contour curr_contour;
  473 +
  474 + // Build the first move
  475 + long int neighbors[4];
  476 + FOURNeighborhood(neighbors, curr_mat_id, bbox.m_W, bbox.m_H);
  477 +
  478 + // Only 2 cases: neighbor at the right or at the bottom
  479 + if(grid[neighbors[1]] > 0) // right
  480 + {
  481 + pos = 1; // top left corner of the pixel
  482 + curr_mat_id++; // go to the pixel at the right
  483 + direction = 1; // direction is along the right
  484 + curr_contour.push_back(1); // Add the move to the right
  485 + }
  486 + else if(grid[neighbors[2]] > 0) // bottom
  487 + {
  488 + pos = 2; // top right corner of the pixel
  489 + curr_mat_id += bbox.m_W; // go to the pixel at the bottom
  490 + direction = 2; // direction is along the bottom
  491 + curr_contour.push_back(1); // add move to the right
  492 + curr_contour.push_back(2); // add move to the bottom
  493 + }
  494 +
  495 + // Keep going this same reasonning until we reach the start pixel (start_mat_id)
  496 + while(curr_mat_id != start_mat_id)
  497 + {
  498 + UpdateSight(sight, direction);
  499 + FOURNeighborhood(neighbors, curr_mat_id, bbox.m_W, bbox.m_H);
  500 +
  501 + // All the cases are possible: top, right, bottom, left
  502 + for(short d=0; d<4; d++)
  503 + {
  504 + if(neighbors[sight[d]] > -1)
  505 + {
  506 + if(grid[neighbors[sight[d]]] > 0)
  507 + {
  508 + if(sight[d] == 0)
  509 + {
  510 + EncodeContourForTopPixel(curr_mat_id, direction, pos, curr_contour, bbox.m_W);
  511 + break;
  512 + }
  513 + else if(sight[d] == 1)
  514 + {
  515 + EncodeContourForRightPixel(curr_mat_id, direction, pos, curr_contour);
  516 + break;
  517 + }
  518 + else if(sight[d] == 2)
  519 + {
  520 + EncodeContourForBottomPixel(curr_mat_id, direction, pos, curr_contour, bbox.m_W);
  521 + break;
  522 + }
  523 + else if(sight[d] == 3)
  524 + {
  525 + EncodeContourForLeftPixel(curr_mat_id, direction, pos, curr_contour);
  526 + break;
  527 + }
  528 + }
  529 + }
  530 + }
  531 +
  532 + // It is possible to reach the start pixel whithout finishing to encode the whole contour.
  533 + // 00111
  534 + // 11100
  535 + // In this case the direction value is always 3, the position value is always 3
  536 + // We have to check if there is a neighbor at left of the current direction
  537 + if(curr_mat_id == start_mat_id)
  538 + {
  539 + if(pos == 3 && direction == 3)
  540 + {
  541 + UpdateSight(sight, direction);
  542 + FOURNeighborhood(neighbors, curr_mat_id, bbox.m_W, bbox.m_H);
  543 +
  544 + if(neighbors[sight[0]] > -1)
  545 + {
  546 + if(grid[neighbors[sight[0]]] == 1)
  547 + {
  548 + pos = 3;
  549 + direction = 2;
  550 + curr_contour.push_back(2);
  551 + curr_mat_id += bbox.m_W;
  552 + }
  553 + }
  554 + }
  555 + }
  556 + }
  557 +
  558 + // We reach the start pixel but maybe we did not encode the whole contour
  559 + if(pos == 3)
  560 + {
  561 + curr_contour.push_back(3);
  562 + curr_contour.push_back(0);
  563 + }
  564 + else if(pos == 4)
  565 + curr_contour.push_back(0);
  566 +
  567 + return curr_contour;
  568 + }
  569 +
  570 + Contour
  571 + ContourOperations::MergeContours(const long unsigned int start_id1,
  572 + const long unsigned int start_id2,
  573 + const Contour& c1,
  574 + const Contour& c2,
  575 + const BoundingBox& bb,
  576 + const unsigned int width,
  577 + const unsigned int height)
  578 + {
  579 + /* Grid with the dimension of the bounding box initialized to 0 */
  580 + short matrix_bbox[bb.m_W * bb.m_H];
  581 + memset(matrix_bbox, 0, bb.m_W * bb.m_H * sizeof(short));
  582 +
  583 + /* Generate the border pixels wrt to c1 and c2 */
  584 + PixelList pixels1 = GenerateBorderPixels(start_id1, c1, width);
  585 + {
  586 + PixelList pixels2 = GenerateBorderPixels(start_id2, c2, width);
  587 + pixels1.insert(pixels1.end(), pixels2.begin(), pixels2.end());
  588 + }
  589 +
  590 + /* Each case of the grid where a pixel is located is set to 1 */
  591 + for(PixelConstIterator pit = pixels1.begin(); pit != pixels1.end(); ++pit)
  592 + matrix_bbox[ToBoundingBoxCoordinates(*pit, bb, width)] = 1;
  593 +
  594 + /* Remove internal pixels */
  595 + long unsigned int bb_id;
  596 + bool is_internal;
  597 + long int neighbors[8];
  598 + for(long unsigned int y = 1; y < bb.m_H - 1; ++y)
  599 + {
  600 + for(long unsigned int x = 1; x < bb.m_W - 1; ++x)
  601 + {
  602 + bb_id = y*bb.m_W + x;
  603 + is_internal = true;
  604 + if(matrix_bbox[bb_id] > 0)
  605 + {
  606 + EIGHTNeighborhood(neighbors, bb_id, bb.m_W, bb.m_H);
  607 + for(short j = 0; j < 8; ++j)
  608 + {
  609 + if(matrix_bbox[neighbors[j]] < 1)
  610 + {
  611 + is_internal = false;
  612 + break;
  613 + }
  614 + }
  615 +
  616 + if(is_internal)
  617 + matrix_bbox[bb_id] = 2;
  618 + }
  619 + }
  620 + }
  621 +
  622 + for(long unsigned int i = 0; i < bb.m_H * bb.m_W; ++i)
  623 + {
  624 + if(matrix_bbox[i] > 1)
  625 + matrix_bbox[i] = 0;
  626 + }
  627 +