Commit ff8510c6 authored by Pierre Lassalle's avatar Pierre Lassalle

Remove older grm version

parent 3b980eec
#=========================================================================
# Program: Generic Region Merging Library (GRM)
# Language: C++
# author: Lassalle Pierre
# Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
# See grmlib-copyright.txt for details.
# This software is distributed WITHOUT ANY WARRANTY; without even
# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE. See the above copyright notices for more information.
#=========================================================================
add_executable(RegionMergingSegmentation RegionMergingSegmentation.cxx)
target_link_libraries(RegionMergingSegmentation OTBGRM)
#include <iostream>
#include <otbImage.h>
#include <otbVectorImage.h>
#include <otbImageFileReader.h>
#include <otbImageFileWriter.h>
#include "lsrmBaatzSegmenter.h"
int main(int argc, char *argv[])
{
if(argc != 7)
{
std::cerr << "Usage: ./" << argv[0] << "\n"
<< "\t[input image path] : (.jpg, .png, .tif)\n"
<< "\t[output clustered image] : (.jpg, .png, .tif)\n"
<< "\t[output label image] : (.tif)\n"
<< "\t[spectral weight] : range from 0 to 1\n"
<< "\t[shape weight] : range from 0 to 1\n"
<< "\t[scale threshold] : unlimited positive value"
<< std::endl;
return 1;
}
lsrm::BaatzParam params;
const char * inFileName = argv[1];
const char * clusterFileName = argv[2];
const char * labelFileName = argv[3];
params.m_SpectralWeight = atof(argv[4]);
params.m_ShapeWeight = atof(argv[5]);
float sqrt_scale = atof(argv[6]);
const float scale = sqrt_scale * sqrt_scale;
typedef float PixelType;
typedef unsigned long int LabelPixelType;
typedef unsigned char ClusterPixelType;
typedef otb::VectorImage<PixelType, 2> InputImageType;
typedef otb::Image<LabelPixelType, 2> LabelImageType;
typedef otb::VectorImage<ClusterPixelType, 2> ClusterImageType;
typedef otb::ImageFileReader<InputImageType> InputImageReaderType;
typedef otb::ImageFileWriter<LabelImageType> LabelImageWriterType;
typedef otb::ImageFileWriter<ClusterImageType> ClusterImageWriterType;
typedef lsrm::BaatzSegmenter<InputImageType> SegmenterType;
auto imgReader = InputImageReaderType::New();
imgReader->SetFileName(inFileName);
imgReader->Update();
SegmenterType segmenter;
segmenter.SetParam(params);
segmenter.SetThreshold(scale);
segmenter.SetInput(imgReader->GetOutput());
segmenter.Update();
auto labelWriter = LabelImageWriterType::New();
labelWriter->SetFileName(labelFileName);
labelWriter->SetInput(segmenter.GetLabeledClusteredOutput());
labelWriter->Update();
auto clusterWriter = ClusterImageWriterType::New();
clusterWriter->SetFileName(clusterFileName);
clusterWriter->SetInput(segmenter.GetClusteredImageOutput());
clusterWriter->Update();
return 0;
}
#=========================================================================
# Program: Generic Region Merging Library (GRM)
# Language: C++
# author: Lassalle Pierre
# Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
# See grmlib-copyright.txt for details.
# This software is distributed WITHOUT ANY WARRANTY; without even
# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE. See the above copyright notices for more information.
#=========================================================================
project(GRM)
cmake_minimum_required(VERSION 2.8)
find_package(OTB)
IF(OTB_FOUND)
include(${OTB_USE_FILE})
ELSE(OTB_FOUND)
message(FATAL_ERROR
"Cannot build OTB project without OTB. Please set OTB_DIR.")
ENDIF(OTB_FOUND)
#Activate c++11
include(CheckCXXCompilerFlag)
CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11)
if(COMPILER_SUPPORTS_CXX11)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -fpermissive -Wall -Wmaybe-uninitialized")
else()
message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.")
endif()
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Code)
add_subdirectory(Code)
add_subdirectory(Applications)
#=========================================================================
# Program: Generic Region Merging Library (GRM)
# Language: C++
# author: Lassalle Pierre
# Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
# See grmlib-copyright.txt for details.
# This software is distributed WITHOUT ANY WARRANTY; without even
# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE. See the above copyright notices for more information.
#=========================================================================
add_library(OTBGRM
lsrmNeighborhood.cpp)
target_link_libraries(OTBGRM OTBCommon OTBIO)
#ifndef __LSRM_BAATZ_SEGMENTER_H
#define __LSRM_BAATZ_SEGMENTER_H
#include "lsrmSegmenter.h"
/*
* Tutorial : Implementation of the Baatz & Schape criterion
*
* Details about the criterion can be found in the following publication:
*
* Martin Baatz and Arno Schape. Multiresolution segmentation: an optimization approach for high quality multi-scale image segmentation.
* Angewandte Geographische Informationsverarbeitung XII, pages 12–23, 2000.
*
* The steps are ordered has to be followed in a chronological way for a full understanding.
* This tutorial does not aim at explaining all the details of the large scale region merging
* but gives all the elements to use it properly.
*/
namespace lsrm
{
/*
* Step 1 :
* We define the specific attributes required for the Baatz & Schape criterion.
* Regions are represented by nodes in a graph and the lsrm library has an internal
* representation of the nodes with the class lsrm::Node.
* --> A node contains a unique number to be indentified, it corresponds to the vectorized
* coordinates of its first pixel and its name is m_Id.
* To retrieve the 2D coordinates from the value of m_Id, it is necessary to know the dimension
* of the image (the number of rows and columns), then :
* x = m_Id % columns and y = m_Id / columns.
*
* --> A node contains the perimeter (m_Perimeter) of its corresponding region, it is the length of the external
* boundary of the region. For example a region of one pixel has a perimeter of 4.
*
* --> A node contains the area (m_Area) of its corresponding region, it is the number of pixels within the
* region.
*
* --> A node contains the minimum rectangular bounding box (parrallel to the image axis) of the region (m_Bbox).
* the boundix box is determined by the structure lsrm::BoundingBox with the coordinates of the upper left
* corner pixel and the width and the height (m_UX, m_UY, m_W, m_H).
*
* After Reading the paper about the Baatz & Schape criterion, we can conclude that we have to define additional
* spectral attributes for a node: its mean vector, the square mean vector, the sum of all the pixels
* and the standard deviation vector.
*
* We define a new class BaatzNode which inherits from the template structure Node. Notice that the template
* is the derived node type, then the template type is BaatzNode.
*/
struct BaatzNode : Node<BaatzNode>
{
std::vector<float> m_Means;
std::vector<float> m_SquareMeans;
std::vector<float> m_SpectralSum;
std::vector<float> m_Std;
};
/*
* Step 2
*
* The Baatz & Schape criterion has user-defined parameters, the spectral weight
* which is the relative importance given to the spectral components and the shape
* weight for the shape components.
*
* We define then a structure wrapping these two parameters we called BaatzParam.
*/
struct BaatzParam
{
float m_SpectralWeight;
float m_ShapeWeight;
};
/*
* Step 3 :
*
*/
template<class TImage>
class BaatzSegmenter : public Segmenter< TImage, BaatzNode, BaatzParam>
{
public:
/* Some convenient typedefs */
typedef Segmenter<TImage, BaatzNode, BaatzParam> Superclass;
typedef TImage ImageType;
typedef typename Superclass::GraphType GraphType;
typedef typename Superclass::NodePointerType NodePointerType;
typedef typename Superclass::GraphOperatorType GraphOperatorType;
typedef GraphToOtbImage<GraphType> IOType;
BaatzSegmenter();
void Update();
float ComputeMergingCost(NodePointerType n1, NodePointerType n2);
void UpdateSpecificAttributes(NodePointerType n1, NodePointerType n2);
void InitFromImage();
};
} // end of namespace lsrm
#include "lsrmBaatzSegmenter.txx"
#endif
#ifndef __LSRM_BAATZ_SEGMENTER_TXX
#define __LSRM_BAATZ_SEGMENTER_TXX
#include <otbImageFileReader.h>
#include <itkImageRegionIterator.h>
namespace lsrm
{
template<class TImage>
BaatzSegmenter<TImage>::BaatzSegmenter() : Superclass()
{
this->m_DoBFSegmentation = true;
this->m_NumberOfIterations = 75;
}
template<class TImage>
void
BaatzSegmenter<TImage>::InitFromImage()
{
typedef itk::ImageRegionIterator<TImage> ImageIterator;
this->m_ImageWidth = this->m_InputImage->GetLargestPossibleRegion().GetSize()[0];
this->m_ImageHeight =this->m_InputImage->GetLargestPossibleRegion().GetSize()[1];
this->m_NumberOfComponentsPerPixel = this->m_InputImage->GetNumberOfComponentsPerPixel();
std::size_t idx = 0;
ImageIterator it(this->m_InputImage, this->m_InputImage->GetLargestPossibleRegion());
for(it.GoToBegin(); !it.IsAtEnd(); ++it)
{
this->m_Graph.m_Nodes[idx]->m_Means.reserve(this->m_NumberOfComponentsPerPixel);
this->m_Graph.m_Nodes[idx]->m_SquareMeans.reserve(this->m_NumberOfComponentsPerPixel);
this->m_Graph.m_Nodes[idx]->m_SpectralSum.reserve(this->m_NumberOfComponentsPerPixel);
this->m_Graph.m_Nodes[idx]->m_Std.assign(this->m_NumberOfComponentsPerPixel, 0.0f);
for(std::size_t b = 0; b < this->m_NumberOfComponentsPerPixel; ++b)
{
this->m_Graph.m_Nodes[idx]->m_Means.push_back(it.Get()[b]);
this->m_Graph.m_Nodes[idx]->m_SquareMeans.push_back((it.Get()[b])*(it.Get()[b]));
this->m_Graph.m_Nodes[idx]->m_SpectralSum.push_back(it.Get()[b]);
}
++idx;
}
}
template<class TImage>
float
BaatzSegmenter<TImage>::ComputeMergingCost(NodePointerType n1, NodePointerType n2)
{
const unsigned int a1 = n1->m_Area, a2 = n2->m_Area, a_sum = a1 + a2;
float spect_cost = 0.0f;
float mean, square_mean, sum, std;
for (unsigned int b = 0; b < this->m_NumberOfComponentsPerPixel; b++)
{
mean = ((a1 * n1->m_Means[b]) + (a2 * n2->m_Means[b])) / a_sum;
square_mean = n1->m_SquareMeans[b] + n2->m_SquareMeans[b];
sum = n1->m_SpectralSum[b] + n2->m_SpectralSum[b];
std = std::sqrt((square_mean - 2*mean*sum + a_sum * mean* mean) / a_sum);
spect_cost += (a_sum * std - a1 * n1->m_Std[b] - a2 * n2->m_Std[b]);
}
spect_cost *= this->m_Param.m_ShapeWeight;
if(spect_cost < this->m_Threshold)
{
float shape_cost, smooth_f, compact_f;
// Compute the shape merging cost
const float p1 = static_cast<float>(n1->m_Perimeter);
const float p2 = static_cast<float>(n2->m_Perimeter);
const unsigned int boundary = (GraphOperatorType::FindEdge(n1, n2))->m_Boundary;
const float p3 = p1 + p2 - 2 * static_cast<float>(boundary);
const BoundingBox merged_bbox = MergeBoundingBoxes(n1->m_Bbox, n2->m_Bbox);
const float bb1_perimeter = static_cast<float>(2*n1->m_Bbox.m_W + 2*n1->m_Bbox.m_H);
const float bb2_perimeter = static_cast<float>(2*n2->m_Bbox.m_W + 2*n2->m_Bbox.m_H);
const float mbb_perimeter = static_cast<float>(2 * merged_bbox.m_W + 2 * merged_bbox.m_H);
smooth_f = a_sum*p3/mbb_perimeter - a1*p1/bb1_perimeter - a2*p2/bb2_perimeter;
compact_f = a_sum*p3/std::sqrt(a_sum) - a1*p1/std::sqrt(a1) - a2*p2/std::sqrt(a2);
shape_cost = this->m_Param.m_ShapeWeight * compact_f + (1-this->m_Param.m_ShapeWeight) * smooth_f;
return (spect_cost + (1-this->m_Param.m_ShapeWeight)*shape_cost);
}
else
return spect_cost;
}
template<class TImage>
void
BaatzSegmenter<TImage>::UpdateSpecificAttributes(NodePointerType n1, NodePointerType n2)
{
const float a1 = static_cast<float>(n1->m_Area);
const float a2 = static_cast<float>(n2->m_Area);
const float a_sum = a1 + a2;
for(unsigned int b = 0; b < this->m_NumberOfComponentsPerPixel; ++b)
{
n1->m_Means[b] = (a1 * n1->m_Means[b] + a2 * n2->m_Means[b]) / a_sum;
n1->m_SquareMeans[b] += n2->m_SquareMeans[b];
n1->m_SpectralSum[b] += n2->m_SpectralSum[b];
n1->m_Std[b] = std::sqrt((n1->m_SquareMeans[b] - 2 * n1->m_Means[b] * n1->m_SpectralSum[b] +
a_sum * n1->m_Means[b] * n1->m_Means[b]) / a_sum);
}
}
template<class TImage>
void
BaatzSegmenter<TImage>::Update()
{
GraphOperatorType::InitNodes(this->m_InputImage, this->m_Graph, *this, FOUR);
bool prev_merged = false;
if(this->m_NumberOfIterations > 0)
{
prev_merged =
GraphOperatorType::PerfomAllIterationsWithLMBFAndConstThreshold(this->m_Graph, *this,
this->m_Threshold, this->m_NumberOfIterations,
this->m_ImageWidth, this->m_ImageHeight);
if(prev_merged && this->m_DoBFSegmentation)
{
prev_merged =
GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(this->m_Graph, *this,
this->m_Threshold, this->m_NumberOfIterations,
this->m_ImageWidth, this->m_ImageHeight);
}
}
else
{
assert(this->m_DoBFSegmentation == true);
prev_merged =
GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(this->m_Graph, *this,
this->m_Threshold, this->m_NumberOfIterations,
this->m_ImageWidth, this->m_ImageHeight);
}
}
} // end of namespace lsrm
#endif
#ifndef __LSRM_DATA_STRUCTURES_H
#define __LSRM_DATA_STRUCTURES_H
#include <stdexcept>
#include <bitset>
#include <vector>
namespace lsrm
{
struct BoundingBox
{
/* X coordinate of the upper left. */
long unsigned int m_UX;
/* Y coordinate of the upper left. */
long unsigned int m_UY;
/* Width */
unsigned int m_W;
/* Height */
unsigned int m_H;
};
BoundingBox MergeBoundingBoxes(const BoundingBox& bb1,
const BoundingBox& bb2)
{
long unsigned int min_ux, min_uy, max_xw, max_yh;
BoundingBox bb;
min_ux = std::min(bb1.m_UX, bb2.m_UX);
min_uy = std::min(bb1.m_UY, bb2.m_UY);
max_xw = std::max(bb1.m_UX + bb1.m_W, bb2.m_UX + bb2.m_W);
max_yh = std::max(bb1.m_UY + bb1.m_H, bb2.m_UY + bb2.m_H);
bb.m_UX = min_ux;
bb.m_UY = min_uy;
bb.m_W = max_xw - min_ux;
bb.m_H = max_yh - min_uy;
return bb;
}
} // end of namespace lsrm
#endif
#ifndef __LSRM_GRAPH_H
#define __LSRM_GRAPH_H
#include "lsrmDataStructures.h"
namespace lsrm
{
struct BaseNode
{
/* Node already merged. */
bool m_Valid;
/* Node has to be removed from the graph. */
bool m_Expired;
/* Does the node merge at the previous iteration */
bool m_IsMerged;
/* Perimeter of the region */
unsigned int m_Perimeter;
/* Area (number of inner pixels) of the region */
unsigned int m_Area;
/*
Node is identified by the location
of the first pixel of the region.
*/
long unsigned int m_Id;
/*
Bounding box of the region
in the image.
*/
BoundingBox m_Bbox;
/* List of pixels contained in the regions */
std::vector<unsigned long int> m_Pixels;
};
template<class DerivedNode>
struct NeighborType
{
typedef std::weak_ptr<DerivedNode> WeakDerived;
typedef std::shared_ptr<DerivedNode> SharedDerived;
WeakDerived m_Target;
float m_Cost;
unsigned int m_Boundary;
bool m_CostUpdated;
NeighborType(WeakDerived ptr, double w, unsigned int c) :
m_Target(ptr), m_Cost(w), m_Boundary(c), m_CostUpdated(false) {}
inline SharedDerived GetRegion()
{
SharedDerived ptr(m_Target.lock());
if(!ptr)
throw std::runtime_error("lss_GenericLMBFRegionMergingHandler.h - NeighborType::GetRegion - Region pointer is not valid");
return ptr;
}
};
template<class DerivedNode>
struct Node : BaseNode
{
typedef NeighborType<DerivedNode> CRPTNeighborType;
std::vector<CRPTNeighborType> m_Edges;
};
template<class TNode>
struct Graph
{
typedef TNode NodeType;
typedef std::shared_ptr<NodeType> NodePointerType;
typedef typename NodeType::CRPTNeighborType EdgeType;
typedef std::vector<NodePointerType> NodeListType;
typedef typename NodeListType::iterator NodeIteratorType;
typedef typename NodeListType::const_iterator NodeConstIteratorType;
typedef std::vector<EdgeType> EdgeListType;
typedef typename EdgeListType::iterator EdgeIteratorType;
typedef typename EdgeListType::const_iterator EdgeConstIteratorType;
std::vector< NodePointerType > m_Nodes;
};
} // end of namespace lsrm
#endif
#ifndef __LSRM_GRAPH_OPERATIONS_H
#define __LSRM_GRAPH_OPERATIONS_H
#include "lsrmGraph.h"
#include "lsrmNeighborhood.h"
#include <iostream>
#include <cassert>
#include <limits>
#include <map>
#include <utility>
#include <set>
namespace lsrm
{
template<class TSegmenter>
class GraphOperations
{
public:
/* Some convenient typedefs */
typedef TSegmenter SegmenterType;
typedef typename SegmenterType::ImageType ImageType;
typedef typename SegmenterType::GraphType GraphType;
typedef typename GraphType::NodeType NodeType;
typedef typename GraphType::EdgeType EdgeType;
typedef typename GraphType::NodePointerType NodePointerType;
typedef typename GraphType::NodeListType NodeList;
typedef typename GraphType::NodeIteratorType NodeIterator;
typedef typename GraphType::NodeConstIteratorType NodeConstIterator;
typedef typename GraphType::EdgeListType EdgeList;
typedef typename GraphType::EdgeIteratorType EdgeIterator;
typedef typename GraphType::EdgeConstIteratorType EdgeConstIterator;
/*
* Given the size of the input image and the mask of the
* neighborhood, we initialize a new graph of nodes
*
* @params:
* GraphType& graph: reference to a graph of nodes
* const unsigned int width: width of the input image
* const unsigned int height: height of the input image
* CONNECTIVITY mask : mask of the neighborhood (4X4 or 8X8)
*/
static void InitNodes(ImageType * inputImg,
GraphType& graph,
SegmenterType& seg,
CONNECTIVITY mask);
/*
* Given a graph of nodes, we explore all the nodes
* and for each node we compute his merging costs
* with all its neighboring nodes given a function
* to compute the merging cost between two nodes.
*
* @params:
* GraphType& graph: reference to the graph of nodes
* float(*fptr)(NodeType*, NodeType*): pointer to the function
* to compute the merging cost between two adjacent nodes.
*/
static void UpdateMergingCosts(GraphType& graph, SegmenterType& seg);
/*
* Given a node A, we analyse its best node B.
* If the node A is also node B's best node
* then it returns a pointer to node B if node A 's id
* is smaller or a pointer to node A if node B's id is
* smaller
* else it returns a null pointer.
* (Local Mutual Best Fitting Heuristic)
*
* @params:
* NodeType * a : Pointer to the node A
* float t : threshold of the criterion
*/
static NodePointerType CheckLMBF(NodePointerType, float t);
/*
* Given a node A, we analyse its best node B.
* If the criterion is checked and the node B
* is valid then it returns a pointer to the
* node A if node's A id is smaller or a pointer
* to node B if node B's id is smaller
* else it returns a null pointer.
*
* @params:
* NodeType * a : pointer to node A
* float t : threshold of the criterion
*/
static NodePointerType CheckBF(NodePointerType a, float t);
/*
* Given the current node and the target node, it returns
* the edge from the current node targeting to the target
* node.
*
* @params
* const NodeType * n : pointer to the current node.
* const NodeType * target : pointer to the target node.
* @return an iterator pointing to the candidate edge.
*/
static EdgeIterator FindEdge(NodePointerType n, NodePointerType target);
/*
* Given a node a and the node b to be merged into node a,
* it updates the neighbors of node a with respect to the
* neighbors of node b.
*
* @params
* NodeType * a : pointer to node a.
* NodeType * b : pointer to node b.
*/
static void UpdateNeighbors(NodePointerType a, NodePointerType b);
/*
* Given 2 nodes A and B (node B being merged into node A)
* we update the internal attributes of node A with respect
* to node B.
*
* @params:
* NodeType * a: pointer to node A.
* NodeType * b: pointer to node B.
*/
static void UpdateInternalAttributes(NodePointerType a,
NodePointerType b,
const unsigned int width,
const unsigned int height);
/*
* Given a graph, it removes all the expired nodes.
*
* @params
* GraphType& graph : reference to the graph.
*/
static void RemoveExpiredNodes(GraphType& graph);
/*
* Given a graph, a region merging algorithm, a threshold
* and the dimension of the image, it performs one iteration
* of the merging process using the local mutual best fitting
* heuristic.
*
* @params
* GraphType& graph : reference to the graph
* SegmenterType& seg : reference to the region merging algorithm.
* const float threshold : threshold for this iteration.
* const unsigned int width : width of the image.
* const unsigned int height : height of the image.
*
* @return a boolean pointing out if there was at least a fusion
* of nodes.
*/
static bool PerfomOneIterationWithLMBF(GraphType& graph,
SegmenterType& seg,
const float threshold,
const unsigned int width,
const unsigned int height);
/*
* Given a graph, a region merging algorithm, a threshold
* and the dimension of the image, it performs one iteration
* of the merging process using the best fitting
* heuristic.
*
* @params
* GraphType& graph : reference to the graph
* SegmenterType& seg : reference to the region merging algorithm.
* const float threshold : threshold for this iteration.
* const unsigned int width : width of the image.
* const unsigned int height : height of the image.
*
* @return a boolean pointing out if there was at least a fusion
* of nodes.
*/
static bool PerfomOneIterationWithBF(GraphType& graph,
SegmenterType& seg,
const float threshold,
const unsigned int width,
const unsigned int height);
/*
* Given a graph, a region merging algorithm, a threshold,
* the number of iterations to apply and the dimension of the image,
* it performs all the iterations of the merging process using the
* local mutual best fitting heuristic.
* This method can be used when the threshold is constant during
* the region merging process.
*