diff --git a/ChangeLog.md b/ChangeLog.md index 529d824a..8bd2c570 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -5,7 +5,9 @@ domain equal to the input image (instead using the resulting bouding box set). (Bertrand Kerautret [#431](https://github.com/DGtal-team/DGtalTools/pull/431)) - + - criticalKernelsThinning3d: it can now export OBJ files for input surface and + output skeletons + (Jacques-Olivier Lachaud [#414](https://github.com/DGtal-team/DGtalTools/pull/414)) - *converters* - heightfield2shading: new option to add a matcap rendering (from normal @@ -36,6 +38,7 @@ longvol exporting). (Bertrand Kerautret [#420](https://github.com/DGtal-team/DGtalTools/pull/420)) - volInfo: get information from a volumetric file. (David Coeurjolly, [#430](https://github.com/DGtal-team/DGtalTools/pull/430)) + # DGtalTools 1.2 - *global* - Fix itk2vol and fix ITK cmake configuaration that was making issues with the ITK image read. diff --git a/doc/images/resPolygonDigitizer.png b/doc/images/resPolygonDigitizer.png new file mode 100644 index 00000000..5982d6da Binary files /dev/null and b/doc/images/resPolygonDigitizer.png differ diff --git a/doc/images/resPolygonDigitizer2.png b/doc/images/resPolygonDigitizer2.png new file mode 100644 index 00000000..7437e1b6 Binary files /dev/null and b/doc/images/resPolygonDigitizer2.png differ diff --git a/generators/2dSimplePolygonDigitizer.cpp b/generators/2dSimplePolygonDigitizer.cpp new file mode 100644 index 00000000..43c2af08 --- /dev/null +++ b/generators/2dSimplePolygonDigitizer.cpp @@ -0,0 +1,346 @@ +/** + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + **/ +/** + * @file 2dSimplePolygonDigitizer.cpp + * @ingroup Tools + * @author Phuc Ngo (\c hoai-diem-phuc.ngo@loria.fr) + * LORIA - Lorraine Univeristy , France + * + * @date 2021/04/06 + * + * Gauss Digitization of a simple polyline. + * + * This file is part of the DGtal library. + */ + +/////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include +#include + +#include "CLI11.hpp" + +#include "DGtal/base/Common.h" +#include "DGtal/helpers/StdDefs.h" + +//image +#include "DGtal/io/writers/GenericWriter.h" +#include "DGtal/images/ImageSelector.h" +#include "DGtal/io/boards/Board2D.h" + +//contour +#include "DGtal/io/readers/PointListReader.h" + +using namespace DGtal; + + +/** + @page 2dSimplePolygonDigitizer 2dSimplePolygonDigitizer + @brief Compute the Gauss Digitization of a simple closed polyline (no hole or self-intersection). + + The digitizer compute the set of integer points inside the input polyline. + + @b Usage: 2dSimplePolygonDigitizer [input] [output] + + @b Allowed @b options @b are: + + @code + + Positionals: + 1 TEXT:FILE REQUIRED Input polyline filename (sdp). + + Options: + + Positionals: + 1 TEXT:FILE REQUIRED Input polyline file name. + + Options: + -h,--help Print this help message and exit + -i,--input TEXT:FILE REQUIRED Input sdp filename. + -o,--output TEXT=result.pgm the output image filename (pgm or svg) + @endcode + + @b Example: + @code + $2dSimplePolygonDigitizer -i ${DGtal}/tests/samples/flower-30-8-3.sdp -o sample.pgm + @endcode + You will obtain such image: + @image html resPolygonDigitizer.png "Resulting image" + + @b Example @b with @b more @b complex @b contours: + + The file located in $DGtal/examples/samples/contourS.sdp + + @code + $ 2dSimplePolygonDigitizer -i $DGtal/examples/samples/contourS.sdp -o sample2.pgm + @endcode + + You will obtain such image: + @image html resPolygonDigitizer2.png "Resulting image" + + @see @ref 2dSimplePolygonDigitizer + @ref 2dSimplePolygonDigitizer.cpp + + */ + +std::pair getBoundingBox(const std::vector& vecPts); +std::vector GaussDigization(const std::vector& polygon); + +int main( int argc, char** argv ) +{ + + // parse command line using CLI ---------------------------------------------- + CLI::App app; + std::string inputFileName; + std::string outputFileName="result.pgm"; + + app.description("compute the set of integer points inside the input polyline.\n Basic example:\n \t 2dSimplePolygonDigitizer -i inputPolyline.sdp -o contourDisplay.pgm"); + app.add_option("-i,--input,1", inputFileName, "Input filename (freeman chain of sdp)." ) + ->required() + ->check(CLI::ExistingFile); + app.add_option("-o,--output,2", outputFileName, "Output filename", true); + + app.get_formatter()->column_width(40); + CLI11_PARSE(app, argc, argv); + // END parse command line using CLI ---------------------------------------------- + + std::string inputExt = inputFileName.substr(inputFileName.find_last_of(".")+1); + if(inputExt != "sdp"){ + trace.error() << "input file should be sdp file" << std::endl; + return EXIT_FAILURE; + } + + std::vector poly = PointListReader::getPointsFromFile(inputFileName); + std::vector gd = GaussDigization(poly); + + std::string outputExt = outputFileName.substr(outputFileName.find_last_of(".")+1); + std::string outputBaseName = outputFileName.substr(0, outputFileName.find_last_of(".")); + + if(outputExt == "pgm") { + //Write results to pgm image + std::pair bb = getBoundingBox(gd); + typedef ImageSelector ::Type Image; + Z2i::Point p1 = bb.first-Z2i::Point(1,1); + Z2i::Point p2 = bb.second+Z2i::Point(1,1); + Image image(Z2i::Domain(p1,p2)); + for(std::vector::const_iterator it=gd.begin(); it!=gd.end(); it++) { + Z2i::Point p ((*it)[0],(*it)[1]); + image.setValue(p,255); + } + PGMWriter::exportPGM(outputFileName,image); + } + else { + //Write results to svg image + outputFileName = outputBaseName + ".svg"; + DGtal::Board2D board; + board << DGtal::SetMode("PointVector", "Both"); + for(std::vector::const_iterator it=gd.begin(); it!=gd.end(); it++) { + board << *it; + } + board.setPenColor(DGtal::Color::Red); + board.setLineWidth(2); + for(size_t it=0; it getBoundingBox(const std::vector& vecPts) { + int minX = vecPts.at(0)[0]; + int maxX = vecPts.at(0)[0]; + int minY = vecPts.at(0)[1]; + int maxY = vecPts.at(0)[1]; + + for(size_t it=1; it x) minX = x; + if(maxX < x) maxX = x; + if(minY > y) minY = y; + if(maxY < y) maxY = y; + } + Z2i::Point tl (minX,minY); + Z2i::Point rb (maxX,maxY); + return std::make_pair(tl,rb); +} + +/* Source adapted from : https://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/ */ +// Given three colinear points p, q, r, the function checks if +// point q lies on line segment 'pr' +template +T max(T v1, T v2) { + return v1 > v2 ? v1 : v2; +} +template +T min(T v1, T v2) { + return v1 > v2 ? v2 : v1; +} +template +bool onSegment(TPoint p, TPoint q, TPoint r) +{ + if (q[0] <= max(p[0], r[0]) && q[0] >= min(p[0], r[0]) && q[1] <= max(p[1], r[1]) && q[1] >= min(p[1], r[1])) + return true; + return false; +} + +// To find orientation of ordered triplet (p, q, r). +// The function returns following values +// 0 --> p, q and r are colinear +// 1 --> Clockwise +// 2 --> Counterclockwise +template +int orientation(TPoint p, TPoint q, TPoint r) +{ + int val = (q[1] - p[1]) * (r[0] - q[0]) - (q[0] - p[0]) * (r[1] - q[1]); + if (val == 0) return 0; // colinear + return (val > 0)? 1: 2; // clock or counterclock wise +} + +// The function that returns true if line segment 'p1q1' +// and 'p2q2' intersect. +template +bool doIntersect(TPoint p1, TPoint q1, TPoint p2, TPoint q2) +{ + // Find the four orientations needed for general and + // special cases + int o1 = orientation(p1, q1, p2); + int o2 = orientation(p1, q1, q2); + int o3 = orientation(p2, q2, p1); + int o4 = orientation(p2, q2, q1); + + // General case + if (o1 != o2 && o3 != o4) + return true; + + // Special Cases + // p1, q1 and p2 are colinear and p2 lies on segment p1q1 + if (o1 == 0 && onSegment(p1, p2, q1)) return true; + + // p1, q1 and p2 are colinear and q2 lies on segment p1q1 + if (o2 == 0 && onSegment(p1, q2, q1)) return true; + + // p2, q2 and p1 are colinear and p1 lies on segment p2q2 + if (o3 == 0 && onSegment(p2, p1, q2)) return true; + + // p2, q2 and q1 are colinear and q1 lies on segment p2q2 + if (o4 == 0 && onSegment(p2, q1, q2)) return true; + + return false; // Doesn't fall in any of the above cases +} + +//Find direction that does not pass through any vertex of the polygon +template +TPoint findDirection(const std::vector& polygon, TPoint p, TPoint d1, TPoint d2 ) +{ + TPoint d = d1; + bool belong = true; + //Find the ray that doest pass by any vertices of polygon + do { + belong = false; + for(size_t it=0; it +bool isInsidePolygon(const std::vector& polygon, TPoint p) +{ + int n = polygon.size(); + // There must be at least 3 vertices in polygon[] + if (n < 3) return false; + + // Create a point for line segment from p to infinite + // Find the direction without containing any polygon vertices + TPoint extreme1 (p[0], INT_MAX/1e6); + TPoint extreme2 (INT_MAX/1e6, p[1]); + TPoint extreme = findDirection(polygon, p, extreme1, extreme2); + + // Check p is a vertex of polygon, then it is inside + if(extreme==p) + return true; + // Count intersections of the above line with sides of polygon + int count = 0, i = 0; + do + { + int next = (i+1)%n; + + // Check if the line segment from 'p' to 'extreme' intersects + // with the line segment from 'polygon[i]' to 'polygon[next]' + if (doIntersect(polygon[i], polygon[next], p, extreme)) + { + // If the point 'p' is colinear with line segment 'i-next', + // then check if it lies on segment. If it lies, return true, + // otherwise false + if (orientation(polygon[i], p, polygon[next]) == 0) + return onSegment(polygon[i], p, polygon[next]); + + count++; + } + i = next; + } while (i != 0); + + // Return true if count is odd, false otherwise + return count&1; // Same as (count%2 == 1) +} + +// Gauss digitizer +std::vector GaussDigization(const std::vector& polygon) +{ + typedef Z2i::Point TPoint; + std::vector gd; + std::pair bb = getBoundingBox(polygon); + int minx = int(bb.first[0])-1; + int miny = int(bb.first[1])-1; + int maxx = int(bb.second[0])+1; + int maxy = int(bb.second[1])+1; + for(int x = minx; x(polygon,p); + if(ok) { + TPoint pi(p[0],p[1]); + gd.push_back(pi); + } + } + } + return gd; +} diff --git a/generators/CMakeLists.txt b/generators/CMakeLists.txt index 2ed65667..790c10a9 100644 --- a/generators/CMakeLists.txt +++ b/generators/CMakeLists.txt @@ -2,6 +2,7 @@ SET(DGTAL_TOOLS_SRC shapeGenerator contourGenerator 3dParametricCurveDigitizer + 2dSimplePolygonDigitizer ) FOREACH(FILE ${DGTAL_TOOLS_SRC}) diff --git a/volumetric/criticalKernelsThinning3D.cpp b/volumetric/criticalKernelsThinning3D.cpp index 0e00d820..628f26fb 100644 --- a/volumetric/criticalKernelsThinning3D.cpp +++ b/volumetric/criticalKernelsThinning3D.cpp @@ -55,7 +55,8 @@ -e,--exportSDP TEXT Export the resulting set of points in a simple (sequence of discrete point (sdp)). -t,--visualize Visualize result in viewer -k,--keepInputDomain Keep the resulting image domain equal to the input image (instead using the resulting bouding box set). - + -O,--exportOBJ TEXT Export the resulting set of points in an OBJ file. + -I,--exportInputOBJ TEXT Export the input set of points in an OBJ file. @endcode @@ -80,6 +81,7 @@ #endif #include #include +#include #include #include #include "DGtal/images/imagesSetsUtils/SetFromImage.h" @@ -88,9 +90,9 @@ #include "DGtal/images/imagesSetsUtils/ImageFromSet.h" #include -#include #include #include +#include #include #include @@ -119,7 +121,9 @@ int main(int argc, char* const argv[]){ string foreground {"black"}; string outputFilenameImg; string outputFilenameSDP; - + string outputFilenameOBJ; + string outputFilenameInputOBJ; + int thresholdMin {0}; int thresholdMax {255}; int persistence {0}; @@ -147,6 +151,8 @@ int main(int argc, char* const argv[]){ app.add_flag("--verbose,-v",verbose, "Verbose output"); app.add_option("--exportImage,-o",outputFilenameImg, "Export the resulting set of points to a image compatible with GenericWriter."); app.add_option("--exportSDP,-e",outputFilenameSDP, "Export the resulting set of points in a simple (sequence of discrete point (sdp))." ); + app.add_option("--exportOBJ,-O",outputFilenameOBJ, "Export the resulting set of points in an OBJ file." ); + app.add_option("--exportInputOBJ,-I",outputFilenameInputOBJ, "Export the input set of points in an OBJ file." ); #ifdef WITH_QGLVIEWER app.add_flag("--visualize,-t", visualize, "Visualize result in viewer"); #endif @@ -274,7 +280,58 @@ int main(int argc, char* const argv[]){ } } - if (outputFilenameImg != "") + + //-------------- export OBJ ------------------------------------------- + if ( outputFilenameInputOBJ != "" ) + { + Board3D< Space,KSpace > board( ks ); + // Display lines that are not in the mesh. + SCell surfel; + board << SetMode3D( surfel.className(), "Basic"); + board.setLineColor( Color::Black ); + board.setFillColor( Color( 200, 200, 255, 255 ) ); + for ( auto p : all_set ) + { + SCell voxel = ks.sSpel( p ); + for ( Dimension k = 0; k < 3; k++ ) { + Point q = p; + q[ k ] += 1; + if ( all_set.find( q ) == all_set.end() ) + board << ks.sIncident( voxel, k, true ); + q[ k ] -= 2; + if ( all_set.find( q ) == all_set.end() ) + board << ks.sIncident( voxel, k, false ); + } + } + board.saveOBJ( outputFilenameInputOBJ ); + } + + //-------------- export OBJ ------------------------------------------- + if ( outputFilenameOBJ != "" ) + { + Board3D< Space,KSpace > board( ks ); + // Display lines that are not in the mesh. + SCell surfel; + board << SetMode3D( surfel.className(), "Basic"); + board.setLineColor( Color::Black ); + board.setFillColor( Color::Red ); + for ( auto p : thin_set ) + { + SCell voxel = ks.sSpel( p ); + for ( Dimension k = 0; k < 3; k++ ) { + Point q = p; + q[ k ] += 1; + if ( thin_set.find( q ) == thin_set.end() ) + board << ks.sIncident( voxel, k, true ); + q[ k ] -= 2; + if ( thin_set.find( q ) == thin_set.end() ) + board << ks.sIncident( voxel, k, false ); + } + } + board.saveOBJ( outputFilenameOBJ ); + } + +if (outputFilenameImg != "") { if(verbose) std::cout << "outputFilename" << outputFilenameImg << std::endl; @@ -290,12 +347,12 @@ int main(int argc, char* const argv[]){ } } #ifdef WITH_QGLVIEWER - if(visualize) +if(visualize) { int argc(1); char** argv(nullptr); QApplication app(argc, argv); - Viewer3D<> viewer; + Viewer3D<> viewer( ks ); viewer.setWindowTitle("criticalKernelsThinning3D"); viewer.show();