Commit b0b808a6 authored by Yannick TANGUY's avatar Yannick TANGUY
Browse files

Initialize new repo for LSD CMLA algorithm

parents
cmake_minimum_required(VERSION 2.8.9)
project(OTBLSD_CMLA)
set(OTBLSD_CMLA_LIBRARIES OTBLSD_CMLA)
#set(BUILD_MODULE_AS_STANDALONE)
#otb_module_impl()
if(NOT OTB_SOURCE_DIR)
find_package(OTB REQUIRED)
list(APPEND CMAKE_MODULE_PATH ${OTB_CMAKE_DIR})
include(OTBModuleExternal)
else()
otb_module_impl()
endif()
OTB_CREATE_APPLICATION(
NAME LSD_CMLA
SOURCES
otbLSDCMLA.cxx
ParallelSegmentsDetector.cxx
RigtAnglesDetector.cxx
LINK_LIBRARIES ${OTBLSD_CMLA_LIBRARIES})
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.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.
=========================================================================*/
#include "otbVectorData.h"
#include "otbVectorDataFileReader.h"
#include "otbVectorDataFileWriter.h"
#include "otbVectorDataToAlignVectorDataFilter.h"
#include "otbStandardWriterWatcher.h"
#include <time.h>
int main(int argc, char * argv[])
{
if (argc != 6)
{
std::cerr << "Usage: " << argv[0] <<
"infname alignshape_name ";
std::cerr <<"angleThreshold distanceThreshold lengthThreshold" << std::endl;
return EXIT_FAILURE;
}
const char * infname = argv[1];
const char * alignshape_name = argv[2];
double angleThreshold = atof(argv[3]);
double distanceThreshold = atof(argv[4]);
double lengthThreshold = atof(argv[5]);
//double ratioThreshold = atof(argv[7]);
time_t begin,end;
time(&begin);
const unsigned int Dimension = 2;
typedef unsigned char PixelType;
typedef double PrecisionType;
typedef otb::VectorData<PrecisionType> VectorDataType;
typedef otb::VectorDataFileReader<VectorDataType> VectorDataReaderType;
typedef otb::VectorDataFileWriter<VectorDataType> VectorDataWriterType;
VectorDataReaderType::Pointer vreader = VectorDataReaderType::New();
vreader->SetFileName(infname);
vreader->GenerateOutputInformation();
VectorDataWriterType::Pointer shapeWriter = VectorDataWriterType::New();
typedef otb::VectorDataToAlignVectorDataFilter<VectorDataType> AlignFilterType;
AlignFilterType::Pointer alignFilter = AlignFilterType::New();
alignFilter->SetInput(vreader->GetOutput());
alignFilter->SetAngleThreshold(angleThreshold);
alignFilter->SetDistanceThreshold(distanceThreshold);
alignFilter->SetLengthThreshold(lengthThreshold);
//alignFilter->SetRatioThreshold(ratioThreshold);
shapeWriter->SetInput(alignFilter->GetOutput());
shapeWriter->SetFileName(alignshape_name);
shapeWriter->Update();
time(&end);
double temps = difftime(end,begin);
int heures = temps/3600;
int minutes = heures*60;
int secondes = minutes*60;
std::cout <<"\n l'execution du programme "<<argv[0]<<" à duré "<<temps<<" secondes\n"<<std::endl;
return EXIT_SUCCESS;
}
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.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.
=========================================================================*/
#include "otbVectorData.h"
#include "otbVectorDataFileReader.h"
#include "otbVectorDataFileWriter.h"
#include "otbVectorDataToRightAngleVectorDataFilter.h"
int main(int argc, char * argv[])
{
if (argc != 5)
{
std::cerr << "Usage: " << argv[0] <<
" infname outfname ";
std::cerr <<"angleThreshold distanceThreshold" << std::endl;
return EXIT_FAILURE;
}
const char * infname = argv[1];
const char * outfname = argv[2];
double angleThreshold = atof(argv[3]);
double distanceThreshold = atof(argv[4]);
typedef otb::VectorData<> VectorDataType;
typedef otb::VectorDataFileReader<VectorDataType> ReaderType;
typedef otb::VectorDataFileWriter<VectorDataType> WriterType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(infname);
reader->GenerateOutputInformation();
WriterType::Pointer writer = WriterType::New();
typedef otb::VectorDataToRightAngleVectorDataFilter<VectorDataType>
RightAngleFilterType;
RightAngleFilterType::Pointer rightAngleFilter = RightAngleFilterType::New();
rightAngleFilter->SetInput(reader->GetOutput());
rightAngleFilter->SetAngleThreshold(angleThreshold);
rightAngleFilter->SetDistanceThreshold(distanceThreshold);
writer->SetInput(rightAngleFilter->GetOutput());
writer->SetFileName(outfname);
writer->Update();
return EXIT_SUCCESS;
}
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.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.
=========================================================================*/
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbImage.h"
#include "otbMultiToMonoChannelExtractROI.h"
#include "otbOGRDataSourceWrapper.h"
#include "otbOGRFeatureWrapper.h"
#include "itkWin32Header.h"
#ifdef _WIN
#define OTB_FILE_SEPARATOR '\\'
#else
#define OTB_FILE_SEPARATOR '/'
#endif
#include <vcl_algorithm.h>
extern "C"
{
#include "lsd.h"
}
namespace otb
{
namespace Wrapper
{
class LSDCMLA : public Application
{
public:
typedef LSDCMLA Self;
typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef FloatImageType ImageType;
typedef ImageType::InternalPixelType ImagePixelType;
typedef FloatVectorImageType VectorImageType;
typedef VectorImageType::InternalPixelType VectorImagePixelType;
typedef otb::MultiToMonoChannelExtractROI<ImagePixelType,ImagePixelType> ExtractROIFilterType;
typedef itk::ImageRegionConstIterator<ImageType> ImageIteratorType;
itkNewMacro(Self);
itkTypeMacro(Merging, otb::Application);
private:
void DoInit()
{
SetName("LSDCMLA");
SetDescription("TODO");
SetDocName("Line Segment Detector (CMLA implementation)");
SetDocLongDescription("TODO");
SetDocLimitations("TODO");
SetDocAuthors("Julien Michel");
SetDocSeeAlso("");
AddParameter(ParameterType_InputImage, "in", "Input image");
SetParameterDescription( "in", "The input image." );
AddParameter(ParameterType_OutputFilename, "out", "Output gis file");
SetParameterDescription( "out", "The vector file containing the output lines." );
AddParameter(ParameterType_Float,"scale","Gaussian scaling");
SetParameterDescription("scale","Scale image by Gaussian filter before processing.");
SetDefaultParameterFloat("scale",0.8);
SetMinimumParameterFloatValue("scale",0);
AddParameter(ParameterType_Float,"sigmacoef","Ratio for sigma of gaussian filter");
SetParameterDescription("sigmacoef","Sigma for Gaussian filter is computed as sigma_coef/scale.");
SetDefaultParameterFloat("sigmacoef",0.6);
SetMinimumParameterFloatValue("sigmacoef",0);
AddParameter(ParameterType_Float,"quant","Gradient quantization error");
SetParameterDescription("quant","Bound to quantization error on the gradient norm.");
SetDefaultParameterFloat("quant",2);
SetMinimumParameterFloatValue("quant",0);
AddParameter(ParameterType_Float,"angth","Gradient angle tolerance");
SetParameterDescription("angth","Gradient angle tolerance in degrees.");
SetDefaultParameterFloat("angth",22.5);
SetMinimumParameterFloatValue("angth",0);
AddParameter(ParameterType_Float,"logeps","Detection threshold");
SetParameterDescription("logeps","Detection threshold");
SetDefaultParameterFloat("logeps",0);
AddParameter(ParameterType_Float,"densityth","Minimal density of regions");
SetParameterDescription("densityth","Minimal density of region points in a rectangle to be accepted.");
SetDefaultParameterFloat("densityth",0.7);
SetMinimumParameterFloatValue("densityth",0);
SetMaximumParameterFloatValue("densityth",1);
AddParameter(ParameterType_Int, "nbins", "Number of bins in ordering");
SetParameterDescription("nbins", "Number of bins in ordering of gradient modulus.");
SetDefaultParameterInt("nbins", 1024);
SetMinimumParameterIntValue("nbins", 1);
AddParameter(ParameterType_Int, "tilesizex", "Size of tiles in pixel (X-axis)");
SetParameterDescription("tilesizex", "Size of tiles along the X-axis.");
SetDefaultParameterInt("tilesizex", 500);
SetMinimumParameterIntValue("tilesizex", 1);
AddParameter(ParameterType_Int, "tilesizey", "Size of tiles in pixel (Y-axis)");
SetParameterDescription("tilesizey", "Size of tiles along the Y-axis.");
SetDefaultParameterInt("tilesizey", 500);
SetMinimumParameterIntValue("tilesizey", 1);
}
void DoUpdateParameters()
{
}
void DoExecute()
{
unsigned long sizeTilesX = GetParameterInt("tilesizex");
unsigned long sizeTilesY = GetParameterInt("tilesizey");
//Acquisition of the input image dimensions
VectorImageType::Pointer imageIn = GetParameterImage("in");
imageIn->UpdateOutputInformation();
unsigned long sizeImageX = imageIn->GetLargestPossibleRegion().GetSize()[0];
unsigned long sizeImageY = imageIn->GetLargestPossibleRegion().GetSize()[1];
unsigned int nbTilesX = sizeImageX/sizeTilesX + (sizeImageX%sizeTilesX > 0 ? 1 : 0);
unsigned int nbTilesY = sizeImageY/sizeTilesY + (sizeImageY%sizeTilesY > 0 ? 1 : 0);
double scale = GetParameterFloat("scale");
double sigma_coef = GetParameterFloat("sigmacoef");
double quant = GetParameterFloat("quant");
double ang_th = GetParameterFloat("angth");
double log_eps = GetParameterFloat("logeps");
double density_th = GetParameterFloat("densityth");
unsigned int n_bins = GetParameterInt("nbins");
otb::ogr::DataSource::Pointer ogrDS;
otb::ogr::Layer layer(NULL, false);
OGRSpatialReference oSRS(imageIn->GetProjectionRef().c_str());
std::vector<std::string> options;
ogrDS = otb::ogr::DataSource::New(GetParameterString("out"), otb::ogr::DataSource::Modes::Overwrite);
std::string layername = itksys::SystemTools::GetFilenameName(GetParameterString("out"));
std::string extension = itksys::SystemTools::GetFilenameLastExtension(GetParameterString("out"));
layername = layername.substr(0,layername.size()-extension.size());
layer = ogrDS->CreateLayer(layername, &oSRS, wkbMultiLineString, options);
OGRFieldDefn fieldWidth("width",OFTReal);
OGRFieldDefn fieldP("p",OFTReal);
OGRFieldDefn fieldLogNFA("-log(NFA)",OFTReal);
layer.CreateField(fieldWidth,true);
layer.CreateField(fieldP,true);
layer.CreateField(fieldLogNFA,true);
otbAppLogINFO(<<"Number of tiles: "<<nbTilesX<<" x "<<nbTilesY);
//Sums calculation per label
unsigned int currentProgress = 0;
unsigned long totalNbSegments = 0;
for(unsigned int row = 0; row < nbTilesY; row++)
{
for(unsigned int column = 0; column < nbTilesX; column++)
{
double progress = 100 * static_cast<double>(column+nbTilesX*row)/static_cast<double>(nbTilesX*nbTilesY);
while(progress >= currentProgress)
{
otbAppLogINFO(<<"Progress: "<<currentProgress<<"% ("<<totalNbSegments<<" found so far)");
currentProgress+=10;
}
unsigned long startX = column*sizeTilesX;
unsigned long startY = row*sizeTilesY;
unsigned long sizeX = vcl_min(sizeTilesX,sizeImageX-startX);
unsigned long sizeY = vcl_min(sizeTilesY,sizeImageY-startY);
//Tiles extraction of the segmented image
ExtractROIFilterType::Pointer imageROI = ExtractROIFilterType::New();
imageROI->SetInput(imageIn);
imageROI->SetStartX(startX);
imageROI->SetStartY(startY);
imageROI->SetSizeX(sizeX);
imageROI->SetSizeY(sizeY);
imageROI->SetChannel(1);
imageROI->Update();
double * image = new double[(size_t)(sizeX*sizeY)];
ImageIteratorType it(imageROI->GetOutput(),imageROI->GetOutput()->GetLargestPossibleRegion());
unsigned long idx = 0;
for(it.GoToBegin();!it.IsAtEnd() && idx < sizeX*sizeY;++it,++idx)
{
image[idx] = it.Get();
}
int n = 0;
int regXNotUsed = 0;
int regYNotUsed = 0;
double * segs = LineSegmentDetection(&n, image, sizeX, sizeY,
scale, sigma_coef, quant,
ang_th, log_eps, density_th,
n_bins,NULL,&regXNotUsed,&regYNotUsed);
totalNbSegments+=n;
for(unsigned int idx = 0; idx < n;++idx)
{
OGRMultiLineString newMultiLineString;
OGRLineString newLineString;
itk::ContinuousIndex<double,2> a,b;
a[0] = segs[idx * 7];
a[1] = segs[idx * 7 + 1];
b[0] = segs[idx * 7 + 2];
b[1] = segs[idx * 7 + 3];
double width = segs[idx * 7 + 4];
double p = segs[idx * 7 + 5];
double lognfa = segs[idx * 7 + 6];
ImageType::PointType pa, pb;
imageROI->GetOutput()->TransformContinuousIndexToPhysicalPoint(a,pa);
imageROI->GetOutput()->TransformContinuousIndexToPhysicalPoint(b,pb);
// We assume isotropic spacing here
width *= vcl_abs(imageROI->GetOutput()->GetSpacing()[0]);
newLineString.addPoint(pa[0],pa[1],0);
newLineString.addPoint(pb[0],pb[1],0);
newMultiLineString.addGeometry(&newLineString);
otb::ogr::Feature newFeature(layer.GetLayerDefn());
newFeature.SetGeometry(&newMultiLineString);
newFeature.ogr().SetField("Width",width);
newFeature.ogr().SetField("-log(NFA)",lognfa);
newFeature.ogr().SetField("p",p);
layer.CreateFeature(newFeature);
}
layer.ogr().CommitTransaction();
// Release memory
delete [] image;
free( (void *) segs );
}
}
ogrDS->SyncToDisk();
otbAppLogINFO(<<"Progress: 100% ("<<totalNbSegments<<" found)");
}
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::LSDCMLA)
/*----------------------------------------------------------------------------
LSD - Line Segment Detector on digital images
This code is part of the following publication and was subject
to peer review:
"LSD: a Line Segment Detector" by Rafael Grompone von Gioi,
Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall,
Image Processing On Line, 2012. DOI:10.5201/ipol.2012.gjmr-lsd
http://dx.doi.org/10.5201/ipol.2012.gjmr-lsd
Copyright (c) 2007-2011 rafael grompone von gioi <grompone@gmail.com>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Affero General Public License as
published by the Free Software Foundation, either version 3 of the
License, or (at your option) any later version.
This program 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 Affero General Public License for more details.
You should have received a copy of the GNU Affero General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
/** @file lsd.h
LSD module header
@author rafael grompone von gioi <grompone@gmail.com>
*/
/*----------------------------------------------------------------------------*/
#ifndef LSD_HEADER
#define LSD_HEADER
/*----------------------------------------------------------------------------*/
/** LSD Full Interface
@param n_out Pointer to an int where LSD will store the number of
line segments detected.
@param img Pointer to input image data. It must be an array of
doubles of size X x Y, and the pixel at coordinates
(x,y) is obtained by img[x+y*X].
@param X X size of the image: the number of columns.
@param Y Y size of the image: the number of rows.
@param scale When different from 1.0, LSD will scale the input image
by 'scale' factor by Gaussian filtering, before detecting
line segments.
Example: if scale=0.8, the input image will be subsampled
to 80% of its size, before the line segment detector
is applied.
Suggested value: 0.8
@param sigma_scale When scale!=1.0, the sigma of the Gaussian filter is:
sigma = sigma_scale / scale, if scale < 1.0
sigma = sigma_scale, if scale >= 1.0
Suggested value: 0.6
@param quant Bound to the quantization error on the gradient norm.
Example: if gray levels are quantized to integer steps,
the gradient (computed by finite differences) error
due to quantization will be bounded by 2.0, as the
worst case is when the error are 1 and -1, that
gives an error of 2.0.
Suggested value: 2.0
@param ang_th Gradient angle tolerance in the region growing
algorithm, in degrees.
Suggested value: 22.5
@param log_eps Detection threshold, accept if -log10(NFA) > log_eps.
The larger the value, the more strict the detector is,
and will result in less detections.
(Note that the 'minus sign' makes that this
behavior is opposite to the one of NFA.)
The value -log10(NFA) is equivalent but more
intuitive than NFA:
- -1.0 gives an average of 10 false detections on noise
- 0.0 gives an average of 1 false detections on noise
- 1.0 gives an average of 0.1 false detections on nose
- 2.0 gives an average of 0.01 false detections on noise
.
Suggested value: 0.0
@param density_th Minimal proportion of 'supporting' points in a rectangle.
Suggested value: 0.7
@param n_bins Number of bins used in the pseudo-ordering of gradient
modulus.
Suggested value: 1024
@param reg_img Optional output: if desired, LSD will return an
int image where each pixel indicates the line segment
to which it belongs. Unused pixels have the value '0',
while the used ones have the number of the line segment,
numbered 1,2,3,..., in the same order as in the
output list. If desired, a non NULL int** pointer must
be assigned, and LSD will make that the pointer point
to an int array of size reg_x x reg_y, where the pixel
value at (x,y) is obtained with (*reg_img)[x+y*reg_x].
Note that the resulting image has the size of the image
used for the processing, that is, the size of the input
image scaled by the given factor 'scale'. If scale!=1
this size differs from XxY and that is the reason why
its value is given by reg_x and reg_y.
Suggested value: NULL
@param reg_x Pointer to an int where LSD will put the X size
'reg_img' image, when asked for.
Suggested value: NULL
@param reg_y Pointer to an int where LSD will put the Y size
'reg_img' image, when asked for.
Suggested value: NULL
@return A double array of size 7 x n_out, containing the list
of line segments detected. The array contains first
7 values of line segment number 1, then the 7 values
of line segment number 2, and so on, and it finish
by the 7 values of line segment number n_out.
The seven values are:
- x1,y1,x2,y2,width,p,-log10(NFA)
.
for a line segment from coordinates (x1,y1) to (x2,y2),
a width 'width', an angle precision of p in (0,1) given
by angle_tolerance/180 degree, and NFA value 'NFA'.
If 'out' is the returned pointer, the 7 values of
line segment number 'n+1' are obtained with
'out[7*n+0]' to 'out[7*n+6]'.
*/
double * LineSegmentDetection( int * n_out,
double * img, int X, int Y,
double scale, double sigma_scale, double quant,
double ang_th, double log_eps, double density_th,
int n_bins,
int ** reg_img, int * reg_x, int * reg_y );
/*----------------------------------------------------------------------------*/
/** LSD Simple Interface with Scale and Region output.
@param n_out Pointer to an int where LSD will store the number of
line segments detected.
@param img Pointer to input image data. It must be an array of
doubles of size X x Y, and the pixel at coordinates
(x,y) is obtained by img[x+y*X].
@param X X size of the image: the number of columns.
@param Y Y size of the image: the number of rows.
@param scale When different from 1.0, LSD will scale the input image
by 'scale' factor by Gaussian filtering, before detecting
line segments.
Example: if scale=0.8, the input image will be subsampled
to 80% of its size, before the line segment detector
is applied.
Suggested value: 0.8
@param reg_img Optional output: if desired, LSD will return an
int image where each pixel indicates the line segment
to which it belongs. Unused pixels have the value '0',