Commit 91ee89f7 authored by Pierre Lassalle's avatar Pierre Lassalle

GRM library

parent db6aa268
GRM Library: Quick Start Tutorial
Requirements:
* Orfeo Toolbox library available here: http://orfeo-toolbox.org/otb/download.html
* Boost library
* CMake (version 2.8 minimum)
On Linux:
Step 1: copy the source code into your favorite path. In the following we assume that
your path is PATH/grm/
Step 2: Go to PATH/grm/
Step 3: mkdir grm_build
Step 4: cd grm_build/
Step 5: Generate the Makefile:
cmake \
-DCMAKE_CXX_FLAGS:STRING=-std=c++11 \
-DOTB_DIR:PATH=/usr/lib/otb \
-DITK_DIR:PATH=/usr/lib/cmake/ITK-4.6 \
-DGDAL_CONFIG:FILEPATH=/usr/bin/gdal-config \
-DGDAL_INCLUDE_DIR:PATH=/usr/include/gdal \
-DGDAL_LIBRARY:FILEPATH=/usr/lib/libgdal.so \
../code/
Step 6: make
Step 7: ./BaatzSegmentation --help
3 region merging criteria have been already inmplemented, up to you to add a new one.
If you are a generous man, you can send me a message to propose your new criterion
at the following email address: lassallepierre34@gmail.com
Enjoy !
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)
#Check compiler version
if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
# require at least gcc 4.8
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 4.8)
message(FATAL_ERROR "GCC version must be at least 4.8!")
endif()
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
# require at least clang 3.2
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.2)
message(FATAL_ERROR "Clang version must be at least 3.2!")
endif()
else()
message(WARNING "You are using an unsupported compiler! Compilation has only been tested with Clang and GCC.")
endif()
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")
else()
message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.")
endif()
file(
GLOB_RECURSE
HEADERS
"src/*.h"
)
file(
GLOB_RECURSE
TEMPLATES
"src/*.txx"
)
file(
GLOB_RECURSE
SOURCES
"src/*.cxx"
)
set(PROJ_INCLUDE_DIRS "")
foreach(_headerFile ${HEADERS})
get_filename_component(_dir ${_headerFile} PATH)
list(APPEND PROJ_INCLUDE_DIRS ${_dir})
endforeach()
foreach(_templateFile ${TEMPLATES})
get_filename_component(_dir ${_templateFile} PATH)
list(APPEND PROJ_INCLUDE_DIRS ${_dir})
endforeach()
list(REMOVE_DUPLICATES PROJ_INCLUDE_DIRS)
include_directories(${PROJ_INCLUDE_DIRS})
ADD_EXECUTABLE(
EuclideanDistanceSegmentation
apps/EuclideanDistanceSegmentation.cxx
${SOURCES}
)
TARGET_LINK_LIBRARIES(EuclideanDistanceSegmentation OTBCommon OTBIO boost_program_options)
ADD_EXECUTABLE(
FLSASegmentation
apps/FLSASegmentation.cxx
${SOURCES}
)
TARGET_LINK_LIBRARIES(FLSASegmentation OTBCommon OTBIO boost_program_options)
ADD_EXECUTABLE(
BaatzSegmentation
apps/BaatzSegmentation.cxx
${SOURCES}
)
TARGET_LINK_LIBRARIES(BaatzSegmentation OTBCommon OTBIO boost_program_options)
/*=========================================================================
Program: Example of application of the Generic Region Merging Library
Language: C++
author: Lassalle Pierre
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
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.
=========================================================================*/
// Boost headers (mandatory)
#include <boost/program_options/cmdline.hpp>
#include <boost/program_options/parsers.hpp>
#include <boost/program_options/options_description.hpp>
#include <boost/program_options/variables_map.hpp>
// Your algorithm header file
#include "baatz-algorithm.h"
namespace po = boost::program_options;
bool init_args(int argc, char ** argv, po::options_description& desc, po::variables_map& vm);
int main(int argc, char **argv)
{
po::options_description desc("Configuration of the Baatz & Schäpe region merging segmentation");
po::variables_map vm;
if(init_args(argc, argv, desc, vm))
{
typedef otb::VectorImage<float, 2> ImageType;
typedef BaatzAlgorithmRM<ImageType> SegmenterType;
BaatzParams params = {vm["cw"].as<float>(), vm["sw"].as<float>(), vm["sp"].as<float>()};
SegmenterType seg_;
seg_.SetInput(vm["input"].as<std::string>());
seg_.SetOutputRGB(vm["output_rgb"].as<std::string>());
seg_.SetParameters(params);
seg_.InitFromImage();
seg_.Segmentation();
seg_.WriteRGBOutputImage();
return 0;
}
else
return 1;
}
bool init_args(int argc, char ** argv, po::options_description& desc, po::variables_map& vm)
{
desc.add_options()
("help", "print mandatory arguments")
("input", po::value<std::string>(), "set input image file")
("output_rgb", po::value<std::string>(), "set output rgb image file")
("cw", po::value<float>(), "set the spectral weight")
("sw", po::value<float>(), "set the shape weight")
("sp", po::value<float>(), "set the scale parameter");
po::store(po::parse_command_line(argc, argv, desc), vm);
po::notify(vm);
if (vm.count("help")) {
std::cout << desc << "\n";
return false;
}
if (!vm.count("input")) {
std::cout << "The input image file was not set (--input).\n";
std::cout << desc << "\n";
return false;
}
if (!vm.count("output_rgb")) {
std::cout << "The output rgb file was not set (--output_rgb).\n";
std::cout << desc << "\n";
return false;
}
if(!vm.count("cw")){
std::cout << "The spectral weight was not set (--cw).\n";
std::cout << desc << "\n";
return false;
}
if(!vm.count("sw")){
std::cout << "The shape weight was not set (--sw).\n";
std::cout << desc << "\n";
return false;
}
if(!vm.count("sp")){
std::cout << "The scale parameter was not set (--sp).\n";
std::cout << desc << "\n";
return false;
}
return true;
}
/*=========================================================================
Program: Example of application of the Generic Region Merging Library
Language: C++
author: Lassalle Pierre
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
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.
=========================================================================*/
// Boost headers (mandatory)
#include <boost/program_options/cmdline.hpp>
#include <boost/program_options/parsers.hpp>
#include <boost/program_options/options_description.hpp>
#include <boost/program_options/variables_map.hpp>
// Your algorithm header file
#include "euc-dist-algorithm.h"
namespace po = boost::program_options;
bool init_args(int argc, char ** argv, po::options_description& desc, po::variables_map& vm);
int main(int argc, char **argv)
{
po::options_description desc("Configuration of the Euclidean distance region merging segmentation");
po::variables_map vm;
if(init_args(argc, argv, desc, vm))
{
typedef otb::VectorImage<float, 2> ImageType;
typedef EuclideanDistanceRM<ImageType> SegmenterType;
EucDistParams params = {vm["euc_thresh"].as<float>(), 0};
SegmenterType seg_;
seg_.SetInput(vm["input"].as<std::string>());
seg_.SetOutputRGB(vm["output_rgb"].as<std::string>());
seg_.SetParameters(params);
seg_.InitFromImage();
seg_.Segmentation();
seg_.WriteRGBOutputImage();
return 0;
}
else
return 1;
}
bool init_args(int argc, char ** argv, po::options_description& desc, po::variables_map& vm)
{
desc.add_options()
("help", "print mandatory arguments")
("input", po::value<std::string>(), "set input image file")
("output_rgb", po::value<std::string>(), "set output rgb image file")
("euc_thresh", po::value<float>(), "set the euclidean threshold");
po::store(po::parse_command_line(argc, argv, desc), vm);
po::notify(vm);
if (vm.count("help")) {
std::cout << desc << "\n";
return false;
}
if (!vm.count("input")) {
std::cout << "The input image file was not set (--input).\n";
std::cout << desc << "\n";
return false;
}
if (!vm.count("output_rgb")) {
std::cout << "The output rgb file was not set (--output_rgb).\n";
std::cout << desc << "\n";
return false;
}
if(!vm.count("euc_thresh")){
std::cout << "The Euclidean threshold was not set (--euc_thresh).\n";
std::cout << desc << "\n";
return false;
}
return true;
}
/*=========================================================================
Program: Example of application of the Generic Region Merging Library
Language: C++
author: Lassalle Pierre
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
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.
=========================================================================*/
// Boost headers (mandatory)
#include <boost/program_options/cmdline.hpp>
#include <boost/program_options/parsers.hpp>
#include <boost/program_options/options_description.hpp>
#include <boost/program_options/variables_map.hpp>
// Your algorithm header file
#include "fls-algorithm.h"
namespace po = boost::program_options;
bool init_args(int argc, char ** argv, po::options_description& desc, po::variables_map& vm);
int main(int argc, char **argv)
{
po::options_description desc("Configuration of the Full Lambda Schedule region merging segmentation");
po::variables_map vm;
if(init_args(argc, argv, desc, vm))
{
typedef otb::VectorImage<float, 2> ImageType;
typedef FLSAlgorithmRM<ImageType> SegmenterType;
FLSParams params = vm["thresh"].as<float>();
SegmenterType seg_;
seg_.SetInput(vm["input"].as<std::string>());
seg_.SetOutputRGB(vm["output_rgb"].as<std::string>());
seg_.SetParameters(params);
seg_.InitFromImage();
seg_.Segmentation();
seg_.WriteRGBOutputImage();
return 0;
}
else
return 1;
}
bool init_args(int argc, char ** argv, po::options_description& desc, po::variables_map& vm)
{
desc.add_options()
("help", "print mandatory arguments")
("input", po::value<std::string>(), "set input image file")
("output_rgb", po::value<std::string>(), "set output rgb image file")
("thresh", po::value<float>(), "set the threshold");
po::store(po::parse_command_line(argc, argv, desc), vm);
po::notify(vm);
if (vm.count("help")) {
std::cout << desc << "\n";
return false;
}
if (!vm.count("input")) {
std::cout << "The input image file was not set (--input).\n";
std::cout << desc << "\n";
return false;
}
if (!vm.count("output_rgb")) {
std::cout << "The output rgb file was not set (--output_rgb).\n";
std::cout << desc << "\n";
return false;
}
if(!vm.count("thresh")){
std::cout << "The threshold was not set (--thresh).\n";
std::cout << desc << "\n";
return false;
}
return true;
}
/*=========================================================================
Program: Implementation of the Baatz & Schäpe criterion using the GRM library
Language: C++
author: Lassalle Pierre
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
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.
=========================================================================*/
#ifndef __BAATZ_ALGORITHM_H
#define __BAATZ_ALGORITHM_H
#include "grm-region.h"
#include "grm-algorithm.h"
class BaatzRegion : public grm::RegionBase<BaatzRegion>
{
public:
float m_Area;
float m_Perimeter;
std::vector<float> m_Means;
std::vector<float> m_SquareMeans;
std::vector<float> m_SpectralSum;
std::vector<float> m_Std;
};
struct BaatzParams
{
float m_SpectralWeight;
float m_ShapeWeight;
float m_Scale;
};
template<class TInputImage>
class BaatzAlgorithmRM : public grm::RegionMergingAlgorithm<TInputImage,
BaatzRegion,
BaatzParams>
{
public:
// Mandatory typedef
typedef grm::RegionMergingAlgorithm<TInputImage, BaatzRegion, BaatzParams> Superclass;
typedef BaatzRegion RegionType;
typedef std::shared_ptr<RegionType> RegionPointerType;
typedef std::vector<RegionPointerType> RegionListType;
typedef BaatzParams ParameterType;
typedef grm::Interface<RegionType> InterfaceType;
typedef typename Superclass::MergingCostFunction MergingCostFunction;
BaatzAlgorithmRM(){}
virtual ~BaatzAlgorithmRM(){}
void InitFromImage();
void UpdateAttribute(RegionPointerType, RegionPointerType);
virtual void Segmentation();
};
#include "baatz-algorithm.txx"
#endif
/*=========================================================================
Program: Implementation of the Baatz & Schäpe criterion using the GRM library
Language: C++
author: Lassalle Pierre
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
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.
=========================================================================*/
#ifndef __BAATZ_ALGORITHM_TXX
#define __BAATZ_ALGORITHM_TXX
#include "baatz-algorithm.h"
#include <otbVectorImage.h>
#include <otbImageFileReader.h>
#include <itkImageRegionIterator.h>
template<class TInputImage>
void
BaatzAlgorithmRM<TInputImage>::InitFromImage()
{
typedef otb::ImageFileReader<TInputImage> InputImageFileReaderType;
typedef itk::ImageRegionIterator<TInputImage> InputImageIteratorType;
// Read the input image
auto input_reader = InputImageFileReaderType::New();
input_reader->SetFileName(this->m_Input);
input_reader->Update();
// Retrieve useful information from the input image
auto regionToProcess = input_reader->GetOutput()->GetLargestPossibleRegion();
unsigned int rows = regionToProcess.GetSize()[1];
unsigned int cols = regionToProcess.GetSize()[0];
unsigned int bands = input_reader->GetOutput()->GetNumberOfComponentsPerPixel();
// Step 1: For efficiency the user is advised to allocate the initial number of regions
this->m_RegionList.reserve(rows*cols);
// Step 2: the user has to intialize the specific attributes of the initial regions
InputImageIteratorType it(input_reader->GetOutput(), regionToProcess);
for(it.GoToBegin(); !it.IsAtEnd(); ++it)
{
auto new_region = std::make_shared<BaatzRegion>();
new_region->m_Area = 1;
new_region->m_Perimeter = 4;
new_region->m_Means.reserve(bands);
new_region->m_SquareMeans.reserve(bands);
new_region->m_SpectralSum.reserve(bands);
new_region->m_Std.reserve(bands);
for(unsigned int b = 0; b<bands; ++b)
{
new_region->m_Means.push_back(it.Get()[b]);
new_region->m_SquareMeans.push_back(pow(it.Get()[b], 2));
new_region->m_SpectralSum.push_back(it.Get()[b]);
new_region->m_Std.push_back(0);
}
// The user has to push the initial regions
this->m_RegionList.push_back(new_region);
}
// Step 3: The neighborhhod and some internal attributes have to be initialized
// the user has to call the method InitRegions
this->InternalInit(rows, cols);
}
// Let assume that regions r1 and r2 have to be merged and r2 is merged into r1
// The user has to describe how the specific attributes of r1 have to be updated
// by defining the mandatory method UpdateAttribute(RegionPointerType, RegionPointerType)
template<class TInputImage>
void
BaatzAlgorithmRM<TInputImage>::UpdateAttribute(RegionPointerType r1, RegionPointerType r2)
{
float a1 = r1->m_Area, a2 = r2->m_Area, a_tot = a1 + a2;
for(unsigned int b=0; b<r1->m_Means.size(); ++b)
{
r1->m_Means[b] = ((a1*r1->m_Means[b]) + (a2*r2->m_Means[b])) / a_tot;
r1->m_SquareMeans[b] += r2->m_SquareMeans[b];
r1->m_SpectralSum[b] += r2->m_SpectralSum[b];
r1->m_Std[b] = std::sqrt((r1->m_SquareMeans[b] -
2*r1->m_Means[b]*r1->m_SpectralSum[b] +
a_tot*r1->m_Means[b]*r1->m_Means[b]) / a_tot);
}
r1->m_Area += r2->m_Area;
r1->m_Perimeter += r2->m_Perimeter - 2 * r1->GetConnections(r2);
}
// This method is obviously mandatory but its implementation depends on some
// subtilities of your algorithm
// What is common is the the iterative merging process of the region, what can be different
// is the stopping criterion and the behavior of the threshold
template<class TInputImage>
void
BaatzAlgorithmRM<TInputImage>::Segmentation()
{
bool prev_merged = true;
float threshold = this->m_Parameters.m_Scale*this->m_Parameters.m_Scale;
float spec_w = this->m_Parameters.m_SpectralWeight;
float shap_w = this->m_Parameters.m_ShapeWeight;
MergingCostFunction cost_func = [&](RegionPointerType r1, RegionPointerType r2)->double
{
unsigned int bands = r1->m_Means.size();
float a1 = r1->m_Area, a2 = r2->m_Area, a_tot = a1 + a2;
// Compute the spectral merging cost
double spec_cost = 0.0;
float mean, square_mean, sum, std, tmp_cost;
for(unsigned int b = 0; b < bands; b++)
{
mean = ((a1*r1->m_Means[b]) + (a2*r2->m_Means[b])) / a_tot;
square_mean = r1->m_SquareMeans[b] + r2->m_SquareMeans[b];
sum = r1->m_SpectralSum[b] + r2->m_SpectralSum[b];
std = std::sqrt((square_mean - 2*mean*sum + a_tot*mean*mean)/a_tot);
tmp_cost = a_tot * std - a1*r1->m_Std[b] - a2*r2->m_Std[b];
spec_cost += tmp_cost;
}
spec_cost *= spec_w;
if(spec_cost < threshold)
{
double shape_cost, smooth_f, compact_f;
// Compute the shape merging cost
float p1 = r1->m_Perimeter, p2 = r2->m_Perimeter,
p3 = p1 + p2 - 2*r1->GetConnections(r2);
auto mmbox = r1->GetResultingBbox(r2);
float bb1 = 2*r1->GetBboxSize(0) + 2*r1->GetBboxSize(1),
bb2 = 2*r2->GetBboxSize(0) + 2*r2->GetBboxSize(1),
bb3 = 2*mmbox[2] + 3*mmbox[3];
smooth_f = a_tot*p3/bb3 - a1*p1/bb1 - a2*p2/bb2;
compact_f = a_tot*p3/std::sqrt(a_tot) - a1*p1/std::sqrt(a1) - a2*p2/std::sqrt(a2);
shape_cost = shap_w * compact_f + (1-shap_w) * smooth_f;
return (spec_cost + (1-spec_w)*shape_cost);
}
else
return spec_cost;
};
while(prev_merged)
{
std::cout << "." << std::flush;
prev_merged = false;
this->m_RMHandler.UpdateMergingCosts(this->m_RegionList, cost_func);
// Do a loop on the regions of the graph
for(auto& r : this->m_RegionList)
{
// For each explored region, we determine if the LMBF is met
auto ref_region = this->m_RMHandler.CheckLMBF(r, threshold);
// If it holds then the pointer to ref_region is not null
if(ref_region != nullptr)
{
// Process to merge the regions
// Step 1: Merge the specific attributes
UpdateAttribute(ref_region, ref_region->GetClosestNeighbor());
// Step 2: Internal update (mandatory)
this->m_RMHandler.Update(ref_region, ref_region->GetClosestNeighbor());
prev_merged = true;
}
}
// After one iteration, you have to remove the expired regions, i.e. regions
// who have merged into a larger one
this->m_RMHandler.RemoveExpiredVertices(this->m_RegionList);
}
std::cout << "\n";
}
#endif
/*=========================================================================
Program: Implementation of the euclidean distance criterion using the GRM library
Language: C++
author: Lassalle Pierre
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
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.
=========================================================================*/
#ifndef __EUC_DIST_ALGORITHM_H
#define __EUC_DIST_ALGORITHM_H
#include "grm-region.h"
#include "grm-algorithm.h"
class EucDistRegion : public grm::RegionBase<EucDistRegion>
{
public:
EucDistRegion(){}
~EucDistRegion(){}
float m_Area;
std::vector<float> m_Means;
};
struct EucDistParams
{
float m_FinalThresh;
unsigned int m_CurrThresh;
};
template<class TInputImage>
class EuclideanDistanceRM : public grm::RegionMergingAlgorithm<TInputImage,
EucDistRegion,
EucDistParams>
{
public:
// Mandatory typedef.
typedef grm::RegionMergingAlgorithm<TInputImage, EucDistRegion, EucDistParams> Superclass;
typedef EucDistRegion RegionType;
typedef std::shared_ptr<RegionType> RegionPointerType;
typedef std::vector<RegionPointerType> RegionListType;
typedef EucDistParams ParameterType;
typedef grm::Interface<RegionType> InterfaceType;
typedef typename Superclass::MergingCostFunction MergingCostFunction;
EuclideanDistanceRM(){}
virtual ~EuclideanDistanceRM(){}
void InitFromImage();
void UpdateAttribute(RegionPointerType, RegionPointerType);