79#include <unordered_map>
81#include <DGtal/io/viewers/PolyscopeViewer.h>
82#include <DGtal/base/Common.h>
83#include <DGtal/helpers/StdDefs.h>
84#include <DGtal/io/Color.h>
85#include <DGtal/io/readers/GenericReader.h>
86#include <DGtal/io/writers/GenericWriter.h>
87#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
88#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
89#include "DGtal/images/ImageSelector.h"
90#include "DGtal/images/imagesSetsUtils/ImageFromSet.h"
92#include <DGtal/topology/SurfelAdjacency.h>
93#include <DGtal/topology/CubicalComplex.h>
94#include <DGtal/topology/CubicalComplexFunctions.h>
96#include <DGtal/topology/VoxelComplex.h>
97#include <DGtal/topology/VoxelComplexFunctions.h>
98#include "DGtal/topology/NeighborhoodConfigurations.h"
99#include "DGtal/topology/tables/NeighborhoodTables.h"
102#include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
103#include "DGtal/geometry/volumes/distance/VoronoiMap.h"
104#include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
109using namespace DGtal;
111using namespace DGtal::Z3i;
113int main(
int argc,
char*
const argv[]){
117 std::string inputFileName;
118 std::string sk_string {
"1isthmus"};
119 string select_string {
"dmax"};
120 string foreground {
"black"};
121 string outputFilenameImg;
122 string outputFilenameSDP;
124 int thresholdMin {0};
125 int thresholdMax {255};
127 bool profile {
false};
128 bool verbose {
false};
129 bool visualize {
true};
130 bool useInputImgToExp {
false};
132 app.description(
"Compute the thinning of a volume using the CriticalKernels framework\nBasic usage: criticalKernelsThinning3D --input <volFileName> --skel <ulti,end, 1isthmus, isthmus> --select [ -f <white,black> -m <minlevel> -M <maxlevel> -v ] [--persistence <value> ] --persistence <value> ] \n options for --skel {ulti end 1isthmus isthmus} \n options for --select = {dmax random first} \n Example: \n criticalKernelsThinning3D --input ${DGtal}/examples/samples/Al.100.vol --select dmax --skel 1isthmus --persistence 1 --visualize --verbose --exportImage ./Al100_dmax_1isthmus_p1.vol \n");
133 app.add_option(
"-i,--input,1", inputFileName,
"Input vol file." )
135 ->check(CLI::ExistingFile);
137 app.add_option(
"--skel,-s", sk_string,
"Type of skeletonization. Options: 1isthmus, isthmus, end, ulti." )
138 -> check(CLI::IsMember({
"ulti",
"end",
"isthmus",
"1isthmus"}));
139 app.add_option(
"--select,-c", select_string,
"Select the ordering for skeletonization. Options: dmax, random, first")
140 -> check(CLI::IsMember({
"random",
"dmax",
"first"}));
141 app.add_option(
"--foreground,-f",foreground,
"Foreground color in binary image" )
142 -> check(CLI::IsMember({
"white",
"black"}));
143 app.add_option(
"--thresholdMin,-m", thresholdMin,
"Threshold min (excluded) to define binary shape" );
144 app.add_option(
"--thresholdMax,-M", thresholdMax,
"Threshold max (included) to define binary shape" );
145 app.add_option(
"--persistence,-p",persistence,
"Persistence value, implies use of persistence algorithm if p>=1" )
146 ->check(CLI::PositiveNumber);
147 app.add_flag(
"--profile", profile,
"Profile algorithm");
148 app.add_flag(
"--verbose,-v",verbose,
"Verbose output");
149 app.add_option(
"--exportImage,-o",outputFilenameImg,
"Export the resulting set of points to a image compatible with GenericWriter.");
150 app.add_option(
"--exportSDP,-e",outputFilenameSDP,
"Export the resulting set of points in a simple (sequence of discrete point (sdp))." );
151 app.add_flag(
"--visualize,-t", visualize,
"Visualize result in viewer");
152 app.add_flag(
"--useInputImgToExp,-k", useInputImgToExp,
"Use input image type to export result (allowing to keep same domain (and same image spacing when using ITK)).");
154 app.get_formatter()->column_width(40);
155 CLI11_PARSE(app, argc, argv);
161 std::cout <<
"Skel: " << sk_string << std::endl;
162 std::cout <<
"Select: " << select_string << std::endl;
163 std::cout <<
"Persistence: " << persistence << std::endl;
164 std::cout <<
"Input: " << inputFileName << std::endl;
166 trace.beginBlock(
"Reading input");
167 using Domain = Z3i::Domain ;
170 using Image = ImageSelector < Z3i::Domain, unsigned char, ITKIMAGEDATA_CONTAINER_I>::Type ;
172 using Image = ImageSelector < Z3i::Domain, unsigned char>::Type ;
176 Image image = GenericReader<Image>::import(inputFileName);
179 DigitalSet image_set (image.domain());
180 SetFromImage<Z3i::DigitalSet>::append<Image>(
182 thresholdMin, thresholdMax);
187 using DigitalTopology = DT26_6;
189 DGtal::DigitalSetByAssociativeContainer<Domain ,
190 std::unordered_set< typename Domain::Point> >;
191 using Complex = DGtal::VoxelComplex<KSpace>;
193 auto & sk = sk_string;
195 KSpace::Point d1( KSpace::Point::diagonal( 1 ) );
196 ks.init(image.domain().lowerBound() - d1 ,
197 image.domain().upperBound() + d1 ,
true);
199 trace.beginBlock(
"construct with table");
201 vc.construct(image_set, functions::loadTable(simplicity::tableSimple26_6 ));
203 trace.beginBlock(
"load isthmus table");
204 boost::dynamic_bitset<> isthmus_table;
206 isthmus_table = *functions::loadTable(isthmusicity::tableIsthmus);
207 else if (sk ==
"1isthmus")
208 isthmus_table = *functions::loadTable(isthmusicity::tableOneIsthmus);
209 auto pointMap = *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
214 std::function< bool(
const Complex&,
const Cell&) > Skel ;
215 if (sk ==
"ulti") Skel = skelUltimate<Complex>;
216 else if (sk ==
"end") Skel = skelEnd<Complex>;
217 else if (sk ==
"isthmus" || sk ==
"1isthmus")
218 Skel = [&isthmus_table, &pointMap](
const Complex & fc,
219 const Complex::Cell & c){
220 return skelWithTable(isthmus_table, pointMap, fc, c);
222 else throw std::runtime_error(
"Invalid skel string");
223 auto start = std::chrono::system_clock::now();
229 trace.beginBlock(
"Create Distance Map");
230 using L3Metric = ExactPredicateLpSeparableMetric<Z3i::Space, 3>;
231 using DT = DistanceTransformation<Z3i::Space, DigitalSet, L3Metric>;
233 DT dt(image.domain(), image_set, l3);
236 std::function< std::pair<typename Complex::Cell, typename Complex::Data>(
const Complex::Clique&) > Select ;
237 auto & sel = select_string;
238 if (sel ==
"random") Select = selectRandom<Complex>;
239 else if (sel ==
"first") Select = selectFirst<Complex>;
240 else if (sel ==
"dmax"){
242 [&dt](
const Complex::Clique & clique){
243 return selectMaxValue<DT, Complex>(dt,clique);
245 }
else throw std::runtime_error(
"Invalid skel string");
247 trace.beginBlock(
"Thinning");
249 if (persistence == 0)
250 vc_new = asymetricThinningScheme< Complex >(
251 vc, Select, Skel, verbose);
253 vc_new = persistenceAsymetricThinningScheme< Complex >(
254 vc, Select, Skel, persistence, verbose);
257 auto end = std::chrono::system_clock::now();
258 auto elapsed = std::chrono::duration_cast<std::chrono::seconds> (end - start) ;
259 if (profile) std::cout <<
"Time elapsed: " << elapsed.count() << std::endl;
262 DigitalSet thin_set(image.domain());
263 vc_new.dumpVoxels(thin_set);
264 const auto & all_set = image_set;
266 if (outputFilenameSDP !=
"")
269 out.open(outputFilenameSDP.c_str());
270 for (
auto &p : thin_set)
272 out << p[0] <<
" " << p[1] <<
" " << p[2] << std::endl;
276 if (outputFilenameImg !=
"")
279 std::cout <<
"outputFilename" << outputFilenameImg << std::endl;
281 unsigned int foreground_value = 255;
282 if (useInputImgToExp){
283 for(
auto p: image.domain()){image.setValue(p, 0);}
284 ImageFromSet<Image>::append(image, thin_set, foreground_value);
285 image >> outputFilenameImg;
287 auto thin_image = ImageFromSet<Image>::create(thin_set, foreground_value,
false, useInputImgToExp);
288 thin_image >> outputFilenameImg;
294 polyscope::options::programName =
"DGtalTools: criticalKernelsThinning3D";
296 PolyscopeViewer<> viewer( ks );
298 viewer.drawColor(Color(255, 255, 255, 255));
302 viewer.drawColor(Color(40, 200, 55, 10));