PathsBuilder.cpp 45.4 KB
Newer Older
1
2
3
4
5
/**
 * calculate obstacles' voronoi edge paths
 * @author Tobias Weber <tweber@ill.fr>
 * @date jun-2021
 * @license GPLv3, see 'LICENSE' file
6
7
8
 *
 * ----------------------------------------------------------------------------
 * TAS-Paths (part of the Takin software suite)
Tobias WEBER's avatar
Tobias WEBER committed
9
 * Copyright (C) 2021  Tobias WEBER (Institut Laue-Langevin (ILL),
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 *                     Grenoble, France).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3 of the License.
 *
 * 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 <http://www.gnu.org/licenses/>.
 * ----------------------------------------------------------------------------
24
25
26
27
 */

#include "PathsBuilder.h"

28
29
30
31
32
33
34
35
36
37
38
39
40
#include <iostream>
#include <thread>
#include <future>
#include <cmath>
#include <cstdint>

#include "mingw_hacks.h"
#include <boost/asio.hpp>
namespace asio = boost::asio;

using t_task = std::packaged_task<void()>;
using t_taskptr = std::shared_ptr<t_task>;

41

Tobias WEBER's avatar
Tobias WEBER committed
42
43
44
45
46
47
48
49
50
// pixel values for various configuration space conditions
#define PIXEL_VALUE_FORBIDDEN_ANGLE  0xf0
#define PIXEL_VALUE_COLLISION        0xff
#define PIXEL_VALUE_NOCOLLISION      0x00


/**
 * constructor
 */
51
PathsBuilder::PathsBuilder()
52
	: m_sigProgress{std::make_shared<t_sig_progress>()}
Tobias WEBER's avatar
Tobias WEBER committed
53
{ }
54
55


Tobias WEBER's avatar
Tobias WEBER committed
56
57
58
/**
 * destructor
 */
59
PathsBuilder::~PathsBuilder()
Tobias WEBER's avatar
Tobias WEBER committed
60
{ }
61
62
63
64
65


void PathsBuilder::Clear()
{
	//m_img.Clear();
66
67
	m_wallsindextree.Clear();

68
69
70
71
72
	m_wallcontours.clear();
	m_fullwallcontours.clear();

	m_lines.clear();
	m_linegroups.clear();
73

74
	m_voro_results.Clear();
75
76
77
}


78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
/**
 * get path length, taking into account the motor speeds
 */
t_real PathsBuilder::GetPathLength(const t_vec2& _vec) const
{
	// directly calculate length if motor speeds are not used
	if(!m_use_motor_speeds)
		return tl2::norm<t_vec2>(_vec);


	// move analysator instead of monochromator?
	bool kf_fixed = true;
	if(m_tascalc)
	{
		if(!std::get<1>(m_tascalc->GetKfix()))
			kf_fixed = false;
	}

Tobias WEBER's avatar
Tobias WEBER committed
96
97
	const Instrument& instr = m_instrspace->GetInstrument();

98
99
	// monochromator 2theta angular speed (alternatively analyser speed if kf is not fixed)
	t_real a2_speed = kf_fixed
Tobias WEBER's avatar
Tobias WEBER committed
100
101
		? instr.GetMonochromator().GetAxisAngleOutSpeed()
		: instr.GetAnalyser().GetAxisAngleOutSpeed();
102
103

	// sample 2theta angular speed
Tobias WEBER's avatar
Tobias WEBER committed
104
	t_real a4_speed = instr.GetSample().GetAxisAngleOutSpeed();
105
106
107
108
109
110
111
112
113

	t_vec2 vec = _vec;
	vec[0] /= a4_speed;
	vec[1] /= a2_speed;

	return tl2::norm<t_vec2>(vec);
}


114
115
116
117
118
119
/**
 * show progress messages on the console
 */
void PathsBuilder::AddConsoleProgressHandler()
{
	auto handler = []
120
		(CalculationState state, t_real progress,
121
122
		const std::string& msg) -> bool
	{
123
124
125
126
127
		std::string statedescr{};
		if(state == CalculationState::STARTED)
			statedescr = " -- START";
		else if(state == CalculationState::FAILED)
			statedescr = " -- FAIL";
128
		else if(state == CalculationState::SUCCEEDED)
129
130
			statedescr = " -- SUCCESS";

Tobias WEBER's avatar
Tobias WEBER committed
131
132
		std::cout << std::fixed << "["
			<< std::setw(3) << (int)(progress * 100.)
133
			<< "%" << statedescr << "] " << msg << std::endl;
134
135
136
137
138
139
140
		return true;
	};

	AddProgressSlot(handler);
}


Tobias WEBER's avatar
Tobias WEBER committed
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
/**
 * convert a pixel of the plot image into the angular range of the plot
 */
t_vec2 PathsBuilder::PixelToAngle(const t_vec2& pix, bool deg, bool inc_sense) const
{
	return PixelToAngle(pix[0], pix[1], deg, inc_sense);
}


/**
 * convert angular coordinates to a pixel in the plot image
 */
t_vec2 PathsBuilder::AngleToPixel(const t_vec2& angle, bool deg, bool inc_sense) const
{
	return AngleToPixel(angle[0], angle[1], deg, inc_sense);
}


Tobias WEBER's avatar
Tobias WEBER committed
159
/**
Tobias WEBER's avatar
Tobias WEBER committed
160
 * convert a pixel of the plot image into the angular range of the plot
Tobias WEBER's avatar
Tobias WEBER committed
161
 */
Tobias WEBER's avatar
Tobias WEBER committed
162
t_vec2 PathsBuilder::PixelToAngle(t_real img_x, t_real img_y, bool deg, bool inc_sense) const
Tobias WEBER's avatar
Tobias WEBER committed
163
{
Tobias WEBER's avatar
Tobias WEBER committed
164
	t_real x = std::lerp(m_sampleScatteringRange[0], m_sampleScatteringRange[1],
Tobias WEBER's avatar
Tobias WEBER committed
165
		img_x / t_real(m_img.GetWidth()));
Tobias WEBER's avatar
Tobias WEBER committed
166
	t_real y = std::lerp(m_monoScatteringRange[0], m_monoScatteringRange[1],
Tobias WEBER's avatar
Tobias WEBER committed
167
		img_y / t_real(m_img.GetHeight()));
Tobias WEBER's avatar
Tobias WEBER committed
168

Tobias WEBER's avatar
Tobias WEBER committed
169
170
171
172
173
174
	if(deg)
	{
		x *= t_real(180) / tl2::pi<t_real>;
		y *= t_real(180) / tl2::pi<t_real>;
	}

Tobias WEBER's avatar
Tobias WEBER committed
175
176
177
178
179
	const t_real *sensesCCW = nullptr;
	if(m_tascalc)
		sensesCCW = m_tascalc->GetScatteringSenses();

	if(inc_sense && sensesCCW)
Tobias WEBER's avatar
Tobias WEBER committed
180
	{
181
182
183
184
185
		// move analysator instead of monochromator?
		std::size_t mono_idx = 0;
		if(!std::get<1>(m_tascalc->GetKfix()))
			mono_idx = 2;

Tobias WEBER's avatar
Tobias WEBER committed
186
		x *= sensesCCW[1];
187
		y *= sensesCCW[mono_idx];
Tobias WEBER's avatar
Tobias WEBER committed
188
189
	}

Tobias WEBER's avatar
Tobias WEBER committed
190
	return tl2::create<t_vec2>({x, y});
Tobias WEBER's avatar
Tobias WEBER committed
191
192
193
194
}


/**
Tobias WEBER's avatar
Tobias WEBER committed
195
 * convert angular coordinates to a pixel in the plot image
Tobias WEBER's avatar
Tobias WEBER committed
196
 */
Tobias WEBER's avatar
Tobias WEBER committed
197
t_vec2 PathsBuilder::AngleToPixel(t_real angle_x, t_real angle_y, bool deg, bool inc_sense) const
Tobias WEBER's avatar
Tobias WEBER committed
198
199
200
201
202
203
204
{
	if(deg)
	{
		angle_x *= tl2::pi<t_real> / t_real(180);
		angle_y *= tl2::pi<t_real> / t_real(180);
	}

Tobias WEBER's avatar
Tobias WEBER committed
205
206
207
208
209
	const t_real *sensesCCW = nullptr;
	if(m_tascalc)
		sensesCCW = m_tascalc->GetScatteringSenses();

	if(inc_sense && sensesCCW)
Tobias WEBER's avatar
Tobias WEBER committed
210
	{
211
212
213
214
215
		// move analysator instead of monochromator?
		std::size_t mono_idx = 0;
		if(!std::get<1>(m_tascalc->GetKfix()))
			mono_idx = 2;

Tobias WEBER's avatar
Tobias WEBER committed
216
		angle_x *= sensesCCW[1];
217
		angle_y *= sensesCCW[mono_idx];
Tobias WEBER's avatar
Tobias WEBER committed
218
219
220
	}


Tobias WEBER's avatar
Tobias WEBER committed
221
	t_real x = std::lerp(t_real(0.), t_real(m_img.GetWidth()),
Tobias WEBER's avatar
Tobias WEBER committed
222
		(angle_x - m_sampleScatteringRange[0]) / (m_sampleScatteringRange[1] - m_sampleScatteringRange[0]));
Tobias WEBER's avatar
Tobias WEBER committed
223
	t_real y = std::lerp(t_real(0.), t_real(m_img.GetHeight()),
Tobias WEBER's avatar
Tobias WEBER committed
224
		(angle_y - m_monoScatteringRange[0]) / (m_monoScatteringRange[1] - m_monoScatteringRange[0]));
Tobias WEBER's avatar
Tobias WEBER committed
225

Tobias WEBER's avatar
Tobias WEBER committed
226
	return tl2::create<t_vec2>({x, y});
Tobias WEBER's avatar
Tobias WEBER committed
227
228
229
}


Tobias WEBER's avatar
Tobias WEBER committed
230
231
232
/**
 * returns the full or the simplified wall contours
 */
Tobias WEBER's avatar
Tobias WEBER committed
233
const std::vector<std::vector<PathsBuilder::t_contourvec>>&
Tobias WEBER's avatar
Tobias WEBER committed
234
235
236
237
238
239
240
241
242
PathsBuilder::GetWallContours(bool full) const
{
	if(full)
		return m_fullwallcontours;

	return m_wallcontours;
}


243
244
245
246
247
248
249
250
251
252
253
254
255
256
/**
 * indicate that a new workflow starts
 */
void PathsBuilder::StartPathMeshWorkflow()
{
	(*m_sigProgress)(CalculationState::STARTED, 1, "Workflow starting.");
}


/**
 * indicate that the current workflow has ended
 */
void PathsBuilder::FinishPathMeshWorkflow(bool success)
{
257
	CalculationState state = success ? CalculationState::SUCCEEDED : CalculationState::FAILED;
258
259
260
261
	(*m_sigProgress)(state, 1, "Workflow has finished.");
}


262
263
/**
 * calculate the obstacle regions in the angular configuration space
264
 * the monochromator a2/a3 variables can alternatively refer to the analyser a5/a6 in case kf is not fixed
265
 */
Tobias WEBER's avatar
Tobias WEBER committed
266
bool PathsBuilder::CalculateConfigSpace(
Tobias WEBER's avatar
Tobias WEBER committed
267
	t_real da2, t_real da4,
Tobias WEBER's avatar
Tobias WEBER committed
268
269
	t_real starta2, t_real enda2,
	t_real starta4, t_real enda4)
270
{
271
	if(!m_instrspace)
Tobias WEBER's avatar
Tobias WEBER committed
272
		return false;
273

Tobias WEBER's avatar
Tobias WEBER committed
274
275
	m_sampleScatteringRange[0] = starta4;
	m_sampleScatteringRange[1] = enda4;
Tobias WEBER's avatar
Tobias WEBER committed
276
277
	m_monoScatteringRange[0] = starta2;
	m_monoScatteringRange[1] = enda2;
Tobias WEBER's avatar
Tobias WEBER committed
278

Tobias WEBER's avatar
Tobias WEBER committed
279
280
	std::ostringstream ostrmsg;
	ostrmsg << "Calculating configuration space in " << m_maxnum_threads << " threads...";
281
	(*m_sigProgress)(CalculationState::STEP_STARTED, 0, ostrmsg.str());
282

Tobias WEBER's avatar
Tobias WEBER committed
283
	const t_real *sensesCCW = nullptr;
284
285
	std::size_t mono_idx = 0;
	bool kf_fixed = true;
Tobias WEBER's avatar
Tobias WEBER committed
286
	if(m_tascalc)
287
	{
Tobias WEBER's avatar
Tobias WEBER committed
288
289
		sensesCCW = m_tascalc->GetScatteringSenses();

290
291
292
293
294
295
296
297
		// move analysator instead of monochromator?
		if(!std::get<1>(m_tascalc->GetKfix()))
		{
			kf_fixed = false;
			mono_idx = 2;
		}
	}

Tobias WEBER's avatar
Tobias WEBER committed
298
299
	const Instrument& instr = m_instrspace->GetInstrument();

300
301
	// analyser angle (alternatively monochromator angle if kf is not fixed)
	t_real a6 = kf_fixed
Tobias WEBER's avatar
Tobias WEBER committed
302
303
		? instr.GetAnalyser().GetAxisAngleOut()	      // a6 or
		: instr.GetMonochromator().GetAxisAngleOut(); // a2
304

305
	// include scattering senses
Tobias WEBER's avatar
Tobias WEBER committed
306
	if(sensesCCW)
307
	{
Tobias WEBER's avatar
Tobias WEBER committed
308
309
310
		da4 *= sensesCCW[1];
		starta4 *= sensesCCW[1];
		enda4 *= sensesCCW[1];
311

312
313
314
		da2 *= sensesCCW[mono_idx];
		starta2 *= sensesCCW[mono_idx];
		enda2 *= sensesCCW[mono_idx];
315
316
317
318
319
	}

	// create colour map and image
	std::size_t img_w = (enda4-starta4) / da4;
	std::size_t img_h = (enda2-starta2) / da2;
Tobias WEBER's avatar
Tobias WEBER committed
320
	//std::cout << "Image size: " << img_w << " x " << img_h << "." << std::endl;
321
322
323
	m_img.Init(img_w, img_h);

	// create thread pool
Tobias WEBER's avatar
Tobias WEBER committed
324
	asio::thread_pool pool(m_maxnum_threads);
325
326
327
328

	std::vector<t_taskptr> tasks;
	tasks.reserve(img_h);

Tobias WEBER's avatar
Tobias WEBER committed
329

330
	// set image pixels
Tobias WEBER's avatar
Tobias WEBER committed
331
	std::atomic<std::size_t> num_pixels = 0;
332
333
	for(std::size_t img_row=0; img_row<img_h; ++img_row)
	{
334
		auto task = [this, img_w, img_row, a6, kf_fixed, &num_pixels]()
335
336
337
338
339
		{
			InstrumentSpace instrspace_cpy = *this->m_instrspace;

			for(std::size_t img_col=0; img_col<img_w; ++img_col)
			{
Tobias WEBER's avatar
Tobias WEBER committed
340
				t_vec2 angle = PixelToAngle(img_col, img_row, false, true);
Tobias WEBER's avatar
Tobias WEBER committed
341
342
				t_real a4 = angle[0];
				t_real a2 = angle[1];
343
344
				t_real a3 = a4 * 0.5;

Tobias WEBER's avatar
Tobias WEBER committed
345
346
				Instrument& instr = instrspace_cpy.GetInstrument();

347
				// set scattering angles (a2 and a6 are flipped in case kf is not fixed)
Tobias WEBER's avatar
Tobias WEBER committed
348
349
350
				instr.GetMonochromator().SetAxisAngleOut(kf_fixed ? a2 : a6);
				instr.GetSample().SetAxisAngleOut(a4);
				instr.GetAnalyser().SetAxisAngleOut(kf_fixed ? a6 : a2);
351

352
				// set crystal angles (a1 and a5 are flipped in case kf is not fixed)
Tobias WEBER's avatar
Tobias WEBER committed
353
354
355
				instr.GetMonochromator().SetAxisAngleInternal(kf_fixed ? 0.5*a2 : 0.5*a6);
				instr.GetSample().SetAxisAngleInternal(a3);
				instr.GetAnalyser().SetAxisAngleInternal(kf_fixed ? 0.5*a6 : 0.5*a2);
356
357

				// set image value
Tobias WEBER's avatar
Tobias WEBER committed
358
359
360
361
				bool angle_ok = instrspace_cpy.CheckAngularLimits();

				if(!angle_ok)
				{
Tobias WEBER's avatar
Tobias WEBER committed
362
					m_img.SetPixel(img_col, img_row, PIXEL_VALUE_FORBIDDEN_ANGLE);
Tobias WEBER's avatar
Tobias WEBER committed
363
364
365
366
				}
				else
				{
					bool colliding = instrspace_cpy.CheckCollision2D();
Tobias WEBER's avatar
Tobias WEBER committed
367
					m_img.SetPixel(img_col, img_row, colliding ? PIXEL_VALUE_COLLISION : PIXEL_VALUE_NOCOLLISION);
Tobias WEBER's avatar
Tobias WEBER committed
368
				}
Tobias WEBER's avatar
Tobias WEBER committed
369
370

				++num_pixels;
371
			}
Tobias WEBER's avatar
Tobias WEBER committed
372
373

			//std::cout << a2/tl2::pi<t_real>*180. << " finished" << std::endl;
374
375
376
377
378
379
380
		};

		t_taskptr taskptr = std::make_shared<t_task>(task);
		tasks.push_back(taskptr);
		asio::post(pool, [taskptr]() { (*taskptr)(); });
	}

Tobias WEBER's avatar
Tobias WEBER committed
381

382
	// get results
Tobias WEBER's avatar
Tobias WEBER committed
383
	std::size_t num_tasks = tasks.size();
384
385
	// send no more than four-percent update signals
	std::size_t signal_skip = num_tasks / 25;
Tobias WEBER's avatar
Tobias WEBER committed
386
387

	for(std::size_t taskidx=0; taskidx<num_tasks; ++taskidx)
388
	{
Tobias WEBER's avatar
Tobias WEBER committed
389
		// prevent sending too many progress signals
390
		if(signal_skip && (taskidx % signal_skip == 0))
391
		{
392
			if(!(*m_sigProgress)(CalculationState::RUNNING, t_real(taskidx) / t_real(num_tasks), ostrmsg.str()))
Tobias WEBER's avatar
Tobias WEBER committed
393
394
395
396
			{
				pool.stop();
				break;
			}
397
		}
398
399

		tasks[taskidx]->get_future().get();
Tobias WEBER's avatar
Tobias WEBER committed
400
		//std::cout << taskidx << " / " << num_tasks << " finished" << std::endl;
401
402
403
	}

	pool.join();
404
	(*m_sigProgress)(CalculationState::STEP_SUCCEEDED, 1, ostrmsg.str());
Tobias WEBER's avatar
Tobias WEBER committed
405

Tobias WEBER's avatar
Tobias WEBER committed
406
	//std::cout << "pixels total: " << img_h*img_w << ", calculated: " << num_pixels << std::endl;
Tobias WEBER's avatar
Tobias WEBER committed
407
	return num_pixels == img_h*img_w;
408
409
410
}


411
412
413
414
415
416
417
418
419
420
/**
 * save all wall position in an index tree for more efficient position lookup
 */
bool PathsBuilder::CalculateWallsIndexTree()
{
	m_wallsindextree = geo::build_closest_pixel_tree<t_contourvec, decltype(m_img)>(m_img);
	return true;
}


421
/**
422
 * calculate the contour lines of the obstacle regions
423
 */
424
bool PathsBuilder::CalculateWallContours(bool simplify, bool convex_split)
425
{
426
	std::string message{"Calculating obstacle contours..."};
427
	(*m_sigProgress)(CalculationState::STEP_STARTED, 0, message);
428

429
	m_wallcontours = geo::trace_boundary<t_contourvec, decltype(m_img)>(m_img);
Tobias WEBER's avatar
Tobias WEBER committed
430
	m_fullwallcontours = m_wallcontours;
431

432
	(*m_sigProgress)(CalculationState::RUNNING, 0.33, message);
433

Tobias WEBER's avatar
Tobias WEBER committed
434
	if(simplify)
435
	{
Tobias WEBER's avatar
Tobias WEBER committed
436
		// iterate and simplify contour groups
Tobias WEBER's avatar
Tobias WEBER committed
437
438
439
		for(auto& contour : m_wallcontours)
		{
			// replace contour with its convex hull
Tobias WEBER's avatar
Tobias WEBER committed
440
			//std::vector<t_vec2> contour_real = tl2::convert<t_vec2, t_contourvec, std::vector>(contour);
Tobias WEBER's avatar
Tobias WEBER committed
441
			//auto [hull_verts, hull_lines, hull_indices]
Tobias WEBER's avatar
Tobias WEBER committed
442
443
			//	= geo::calc_delaunay<t_vec2>(2, contour_real, true);
			//contour = tl2::convert<t_contourvec, t_vec2, std::vector>(hull_verts);
Tobias WEBER's avatar
Tobias WEBER committed
444
445

			// simplify hull contour
446
			geo::simplify_contour<t_contourvec, t_real>(contour, m_simplify_mindist, m_eps_angular, m_eps);
Tobias WEBER's avatar
Tobias WEBER committed
447
		}
448
	}
Tobias WEBER's avatar
Tobias WEBER committed
449

450
	(*m_sigProgress)(CalculationState::RUNNING, 0.66, message);
Tobias WEBER's avatar
Tobias WEBER committed
451

452
453
	if(convex_split)
	{
Tobias WEBER's avatar
Tobias WEBER committed
454
		// convex split
455
		std::vector<std::vector<t_contourvec>> splitcontours;
Tobias WEBER's avatar
Tobias WEBER committed
456
457
458
459
		splitcontours.reserve(m_wallcontours.size()*2);

		for(auto& contour : m_wallcontours)
		{
460
			//std::reverse(contour.begin(), contour.end());
Tobias WEBER's avatar
Tobias WEBER committed
461
462
			auto splitcontour = geo::convex_split<t_contourvec, t_real>(contour, m_eps);
			if(splitcontour.size())
463
			{
Tobias WEBER's avatar
Tobias WEBER committed
464
				for(auto&& poly : splitcontour)
465
466
				{
					//std::reverse(poly.begin(), poly.end());
467
					splitcontours.emplace_back(std::move(poly));
468
				}
469
470
471
472
473
474
			}
			else
			{
				// no split, use original contour
				splitcontours.push_back(contour);
			}
Tobias WEBER's avatar
Tobias WEBER committed
475
476
		}

477
		m_wallcontours = std::move(splitcontours);
478
	}
Tobias WEBER's avatar
Tobias WEBER committed
479

480
	(*m_sigProgress)(CalculationState::STEP_SUCCEEDED, 1, message);
Tobias WEBER's avatar
Tobias WEBER committed
481
	return true;
482
483
484
}


485
486
487
/**
 * calculate lines segments and groups
 */
Tobias WEBER's avatar
Tobias WEBER committed
488
bool PathsBuilder::CalculateLineSegments(bool use_region_function)
489
{
490
	std::string message{"Calculating obstacle line segments..."};
491
	(*m_sigProgress)(CalculationState::STEP_STARTED, 0, message);
492

493
494
	m_lines.clear();
	m_linegroups.clear();
Tobias WEBER's avatar
Tobias WEBER committed
495
496
	m_points_outside_regions.clear();
	m_inverted_regions.clear();
497

Tobias WEBER's avatar
Tobias WEBER committed
498
	// find an arbitrary point outside all obstacles
Tobias WEBER's avatar
Tobias WEBER committed
499
500
	auto find_point_outside_regions = [this]
	(std::size_t x_start = 0, std::size_t y_start = 0,
Tobias WEBER's avatar
Tobias WEBER committed
501
		bool skip_search = false) -> t_vec2
Tobias WEBER's avatar
Tobias WEBER committed
502
	{
Tobias WEBER's avatar
Tobias WEBER committed
503
		t_vec2 point_outside_regions = tl2::create<t_vec2>({-50, -40});
Tobias WEBER's avatar
Tobias WEBER committed
504
505
506
		bool found_point = false;

		if(!skip_search)
Tobias WEBER's avatar
Tobias WEBER committed
507
		{
Tobias WEBER's avatar
Tobias WEBER committed
508
			for(std::size_t y=y_start; y<m_img.GetHeight(); ++y)
Tobias WEBER's avatar
Tobias WEBER committed
509
			{
Tobias WEBER's avatar
Tobias WEBER committed
510
511
				for(std::size_t x=x_start; x<m_img.GetWidth(); ++x)
				{
Tobias WEBER's avatar
Tobias WEBER committed
512
					if(m_img.GetPixel(x, y) == PIXEL_VALUE_NOCOLLISION)
Tobias WEBER's avatar
Tobias WEBER committed
513
					{
Tobias WEBER's avatar
Tobias WEBER committed
514
515
516
						point_outside_regions =
							tl2::create<t_vec2>({
								static_cast<double>(x),
Tobias WEBER's avatar
Tobias WEBER committed
517
518
								static_cast<double>(y)
							});
Tobias WEBER's avatar
Tobias WEBER committed
519
520
521
522
523
524
						found_point = true;
						break;
					}
				}
				if(found_point)
					break;
Tobias WEBER's avatar
Tobias WEBER committed
525
526
			}
		}
Tobias WEBER's avatar
Tobias WEBER committed
527
528
529

		return point_outside_regions;
	};
Tobias WEBER's avatar
Tobias WEBER committed
530

531
532
533
534
	std::size_t totalverts = 0;
	for(const auto& contour : m_wallcontours)
		totalverts += contour.size();

535
	m_lines.reserve(totalverts/2 + 1);
536
	m_linegroups.reserve(m_wallcontours.size());
Tobias WEBER's avatar
Tobias WEBER committed
537
538
539
540
541
542

	if(!use_region_function)
	{
		m_points_outside_regions.reserve(m_wallcontours.size());
		m_inverted_regions.reserve(m_wallcontours.size());
	}
543
544
545
546
547
548
549

	// contour vertices
	std::size_t linectr = 0;
	for(std::size_t contouridx = 0; contouridx < m_wallcontours.size(); ++contouridx)
	{
		const auto& contour = m_wallcontours[contouridx];
		std::size_t groupstart = linectr;
Tobias WEBER's avatar
Tobias WEBER committed
550
		t_contourvec contour_mean = tl2::zero<t_vec2>(2);
551

552
553
554
555
556
557
		for(std::size_t vert1 = 0; vert1 < contour.size(); ++vert1)
		{
			std::size_t vert2 = (vert1 + 1) % contour.size();

			const t_contourvec& vec1 = contour[vert1];
			const t_contourvec& vec2 = contour[vert2];
558
			contour_mean += vec1;
559

Tobias WEBER's avatar
Tobias WEBER committed
560
561
			t_vec2 linevec1 = vec1;
			t_vec2 linevec2 = vec2;
562

Tobias WEBER's avatar
Tobias WEBER committed
563
564
			m_lines.emplace_back(
				std::make_pair(std::move(linevec1), std::move(linevec2)));
565
566
567

			++linectr;
		}
568

569
570
571
572
		contour_mean /= contour.size();

		// move a point on the contour in the direction of the contour mean
		// to get a point inside the contour
Tobias WEBER's avatar
Tobias WEBER committed
573
		//t_contourvec inside_contour = contour[0] + (contour_mean-contour[0]) / 8;
574
575
576
577

		// find a point outside the contour by moving a pixel away from the minimum vertex
		auto [contour_min, contour_max] = tl2::minmax(contour);
		t_contourvec outside_contour = contour_min;
578
		for(int i = 0; i < 2; ++i)
Tobias WEBER's avatar
Tobias WEBER committed
579
			outside_contour[i] -= 1;
580
581
582

		// mark line group start and end index
		std::size_t groupend = linectr;
583

Tobias WEBER's avatar
Tobias WEBER committed
584
585
		// don't include outer bounding region
		// TODO: test if such a region is there
Tobias WEBER's avatar
Tobias WEBER committed
586
587
		if(contouridx > 0)
		{
Tobias WEBER's avatar
Tobias WEBER committed
588
			m_linegroups.emplace_back(std::make_pair(groupstart, groupend));
Tobias WEBER's avatar
Tobias WEBER committed
589

Tobias WEBER's avatar
Tobias WEBER committed
590
591
592
593
594
595
			if(!use_region_function)
			{
				t_vec2 point_outside_regions =
					find_point_outside_regions(contour[0][0], contour[0][1], true);
				m_points_outside_regions.emplace_back(
					std::move(point_outside_regions));
596

Tobias WEBER's avatar
Tobias WEBER committed
597
598
				//auto pix_incontour = m_img.GetPixel(inside_contour[0], inside_contour[1]);
				auto pix_outcontour = m_img.GetPixel(outside_contour[0], outside_contour[1]);
599
600

#ifdef DEBUG
Tobias WEBER's avatar
Tobias WEBER committed
601
602
603
604
605
				std::cout << "contour " << std::dec << contouridx
					<< ", pixel inside " << inside_contour[0] << ", " << inside_contour[1]
					<< ": " << std::hex << int(pix_incontour) << std::dec
					<< "; pixel outside " << outside_contour[0] << ", " << outside_contour[1]
					<< ": " << std::hex << int(pix_outcontour) << std::endl;
606
607
#endif

Tobias WEBER's avatar
Tobias WEBER committed
608
609
				// normal regions encircle forbidden coordinate points
				// inverted regions encircle allowed coordinate points
Tobias WEBER's avatar
Tobias WEBER committed
610
611
				//m_inverted_regions.push_back(pix_incontour == PIXEL_VALUE_NOCOLLISION);
				m_inverted_regions.push_back(pix_outcontour != PIXEL_VALUE_NOCOLLISION);
Tobias WEBER's avatar
Tobias WEBER committed
612
			}
Tobias WEBER's avatar
Tobias WEBER committed
613
		}
614
	}
615

616
	(*m_sigProgress)(CalculationState::STEP_SUCCEEDED, 1, message);
Tobias WEBER's avatar
Tobias WEBER committed
617
	return true;
618
619
620
}


621
622
623
/**
 * calculate the voronoi diagram
 */
624
625
bool PathsBuilder::CalculateVoronoi(bool group_lines, VoronoiBackend backend,
	bool use_region_function)
626
{
627
	std::string message{"Calculating Voronoi diagram..."};
628
	(*m_sigProgress)(CalculationState::STEP_STARTED, 0, message);
629

630
	// is the vertex in a forbidden region?
631
632
633
634
635
636
637
638
639
640
641
642
	std::function<bool(const t_vec2&)> region_func = [this](const t_vec2& vec) -> bool
	{
		if(vec[0] < 0 || vec[1] < 0)
			return true;

		std::size_t x = vec[0];
		std::size_t y = vec[1];

		if(x >= m_img.GetWidth() || y >= m_img.GetHeight())
			return true;

		// an occupied pixel signifies a forbidden region
Tobias WEBER's avatar
Tobias WEBER committed
643
		if(m_img.GetPixel(x, y) != PIXEL_VALUE_NOCOLLISION)
644
645
646
647
648
			return true;

		return false;
	};

649
650
651
652
653
654
655
	// validation function that checks if (voronoi) vertices are far enough from any wall
	std::function<bool(const t_vec2&)> validation_func = [this](const t_vec2& vec) -> bool
	{
		t_real dist_to_walls = GetDistToNearestWall(vec);
		return dist_to_walls >= m_min_angular_dist_to_walls;
	};

656
657
658
659
660
661
662
	geo::VoronoiLinesRegions<t_vec2, t_line> regions{};
	regions.SetGroupLines(group_lines);
	regions.SetRemoveVoronoiVertices(true);
	regions.SetLineGroups(&m_linegroups);
	regions.SetPointsOutsideRegions(&m_points_outside_regions);
	regions.SetInvertedRegions(&m_inverted_regions);
	regions.SetRegionFunc(use_region_function ? &region_func : nullptr);
Tobias WEBER's avatar
Tobias WEBER committed
663
	regions.SetValidateFunc(m_remove_bisectors_below_min_wall_dist ? &validation_func : nullptr);
664

665
666
667
668
	if(backend == VoronoiBackend::BOOST)
	{
		m_voro_results
			= geo::calc_voro<t_vec2, t_line, t_graph>(
669
				m_lines, m_eps, m_voroedge_eps, &regions);
670
	}
Tobias WEBER's avatar
Tobias WEBER committed
671
#ifdef USE_CGAL
672
673
674
675
	else if(backend == VoronoiBackend::CGAL)
	{
		m_voro_results
			= geo::calc_voro_cgal<t_vec2, t_line, t_graph>(
676
				m_lines, m_eps, m_voroedge_eps, &regions);
677
	}
Tobias WEBER's avatar
Tobias WEBER committed
678
#endif
679
680
681
682
683
684
	else
	{
		// invalid backend selected
		(*m_sigProgress)(CalculationState::FAILED, 1, message);
		return false;
	}
685

686
	(*m_sigProgress)(CalculationState::STEP_SUCCEEDED, 1, message);
Tobias WEBER's avatar
Tobias WEBER committed
687
	return true;
688
689
690
691
}


/**
692
 * save the contour line segments to the lines tool
693
 */
694
bool PathsBuilder::SaveToLinesTool(std::ostream& ostr)
695
{
696
	ostr << "<lines2d>\n";
697
698
	std::vector<std::pair<std::size_t, std::size_t>> group_indices;
	group_indices.reserve(m_linegroups.size());
699

700
	// contour vertices
701
	std::size_t vertctr = 0;
702
	ostr << "<vertices>\n";
703
	for(std::size_t contouridx = 0; contouridx < m_linegroups.size(); ++contouridx)
704
	{
705
		const auto& contour = m_linegroups[contouridx];
706
707
		ostr << "\t<!-- contour " << contouridx << " -->\n";

708
709
		std::size_t group_begin = vertctr;

710
		for(std::size_t lineidx=std::get<0>(contour); lineidx<std::get<1>(contour); ++lineidx)
711
		{
712
			const t_line& line = m_lines[lineidx];
713
714

			ostr << "\t<" << vertctr;
715
716
			ostr << " x=\"" << std::get<0>(line)[0] << "\"";
			ostr << " y=\"" << std::get<0>(line)[1] << "\"";
717
718
719
720
			ostr << "/>\n";
			++vertctr;

			ostr << "\t<" << vertctr;
721
			ostr << " x=\"" << std::get<1>(line)[0] << "\"";
Tobias WEBER's avatar
Tobias WEBER committed
722
			ostr << " y=\"" << std::get<1>(line)[1] << "\"";
723
724
725
			ostr << "/>\n\n";
			++vertctr;
		}
726
727
728

		std::size_t group_end = vertctr;
		group_indices.emplace_back(std::make_pair(group_begin, group_end));
729
	}
730
731
	ostr << "</vertices>\n";

Tobias WEBER's avatar
Tobias WEBER committed
732
733
	// contour groups
	ostr << "\n<groups>\n";
734
	for(std::size_t groupidx = 0; groupidx < group_indices.size(); ++groupidx)
Tobias WEBER's avatar
Tobias WEBER committed
735
	{
736
		const auto& group = group_indices[groupidx];
Tobias WEBER's avatar
Tobias WEBER committed
737
738
739
		ostr << "\t<!-- contour " << groupidx << " -->\n";
		ostr << "\t<" << groupidx << ">\n";

740
741
		ostr << "\t\t<begin>" << std::get<0>(group) << "</begin>\n";
		ostr << "\t\t<end>" << std::get<1>(group) << "</end>\n";
Tobias WEBER's avatar
Tobias WEBER committed
742
743
744
745
746

		ostr << "\t</" << groupidx << ">\n\n";
	}
	ostr << "</groups>\n";

747
	// alternatively: contour regions (obsolete)
Tobias WEBER's avatar
Tobias WEBER committed
748
	/*ostr << "\n<regions>\n";
749
750
751
752
	for(std::size_t contouridx = 0; contouridx<m_wallcontours.size(); ++contouridx)
	{
		const auto& contour = m_wallcontours[contouridx];
		ostr << "\t<!-- contour " << contouridx << " -->\n";
Tobias WEBER's avatar
Tobias WEBER committed
753
		ostr << "\t<" << contouridx << ">\n";
754
755
756
757
758

		for(std::size_t vertidx=0; vertidx<contour.size(); ++vertidx)
		{
			const t_contourvec& vec = contour[vertidx];

Tobias WEBER's avatar
Tobias WEBER committed
759
			ostr << "\t\t<" << vertidx;
760
761
762
763
764
			ostr << " x=\"" << vec[0] << "\"";
			ostr << " y=\"" << vec[1] << "\"";
			ostr << "/>\n";
		}

Tobias WEBER's avatar
Tobias WEBER committed
765
		ostr << "\t</" << contouridx << ">\n\n";
766
	}
Tobias WEBER's avatar
Tobias WEBER committed
767
	ostr << "</regions>\n";*/
768
769

	ostr << "</lines2d>" << std::endl;
770
	return true;
771
}
Tobias WEBER's avatar
Tobias WEBER committed
772
773


774
775
776
777
778
779
780
781
782
783
784
785
786
/**
 * save the contour line segments to the lines tool
 */
bool PathsBuilder::SaveToLinesTool(const std::string& filename)
{
	std::ofstream ofstr(filename);
	if(!ofstr)
		return false;

	return SaveToLinesTool(ofstr);
}


Tobias WEBER's avatar
Tobias WEBER committed
787
788
/**
 * find a path from an initial (a2, a4) to a final (a2, a4)
789
 * the monochromator a2/a3 variables can alternatively refer to the analyser a5/a6 in case kf is not fixed
Tobias WEBER's avatar
Tobias WEBER committed
790
 */
Tobias WEBER's avatar
Tobias WEBER committed
791
InstrumentPath PathsBuilder::FindPath(
Tobias WEBER's avatar
Tobias WEBER committed
792
	t_real a2_i, t_real a4_i,
Tobias WEBER's avatar
Tobias WEBER committed
793
794
	t_real a2_f, t_real a4_f,
	PathStrategy pathstragy)
Tobias WEBER's avatar
Tobias WEBER committed
795
{
796
797
798
799
800
801
802
803
804
	InstrumentPath path{};
	path.ok = false;

	// check if start or target point are within obstacles
	{
		if(!m_instrspace)
			return path;

		const t_real *sensesCCW = nullptr;
805
806
		std::size_t mono_idx = 0;
		bool kf_fixed = true;
807
		if(m_tascalc)
808
		{
809
810
			sensesCCW = m_tascalc->GetScatteringSenses();

811
812
813
814
815
816
817
818
			// move analysator instead of monochromator?
			if(!std::get<1>(m_tascalc->GetKfix()))
			{
				kf_fixed = false;
				mono_idx = 2;
			}
		}

819
		InstrumentSpace instrspace_cpy = *this->m_instrspace;
Tobias WEBER's avatar
Tobias WEBER committed
820
		Instrument& instr = instrspace_cpy.GetInstrument();
821
822
823
824
825
826

		// set instrument angles to start point
		t_real a2 = a2_i;
		t_real a4 = a4_i;
		if(sensesCCW)
		{
827
			a2 *= sensesCCW[mono_idx];
828
829
			a4 *= sensesCCW[1];
		}
830
831
832

		if(kf_fixed)
		{
Tobias WEBER's avatar
Tobias WEBER committed
833
834
			instr.GetMonochromator().SetAxisAngleOut(a2);
			instr.GetMonochromator().SetAxisAngleInternal(0.5 * a2);
835
836
837
		}
		else
		{
Tobias WEBER's avatar
Tobias WEBER committed
838
839
			instr.GetAnalyser().SetAxisAngleOut(a2);
			instr.GetAnalyser().SetAxisAngleInternal(0.5 * a2);
840
		}
Tobias WEBER's avatar
Tobias WEBER committed
841
		instr.GetSample().SetAxisAngleOut(a4);
842
843
844
845
846
847
848
849
850
851
852

		bool in_angular_limits = instrspace_cpy.CheckAngularLimits();
		bool colliding = instrspace_cpy.CheckCollision2D();
		if(!in_angular_limits || colliding)
			return path;

		// set instrument angles to target point
		a2 = a2_f;
		a4 = a4_f;
		if(sensesCCW)
		{
853
			a2 *= sensesCCW[mono_idx];
854
855
			a4 *= sensesCCW[1];
		}
856
857
858

		if(kf_fixed)
		{
Tobias WEBER's avatar
Tobias WEBER committed
859
860
			instr.GetMonochromator().SetAxisAngleOut(a2);
			instr.GetMonochromator().SetAxisAngleInternal(0.5 * a2);
861
862
863
		}
		else
		{
Tobias WEBER's avatar
Tobias WEBER committed
864
865
			instr.GetAnalyser().SetAxisAngleOut(a2);
			instr.GetAnalyser().SetAxisAngleInternal(0.5 * a2);
866
		}
Tobias WEBER's avatar
Tobias WEBER committed
867
		instr.GetSample().SetAxisAngleOut(a4);
868
869
870
871
872
873
874

		in_angular_limits = instrspace_cpy.CheckAngularLimits();
		colliding = instrspace_cpy.CheckCollision2D();
		if(!in_angular_limits || colliding)
			return path;
	}

875

876
	// convert angles to degrees
877
878
879
880
	a2_i *= 180. / tl2::pi<t_real>;
	a4_i *= 180. / tl2::pi<t_real>;
	a2_f *= 180. / tl2::pi<t_real>;
	a4_f *= 180. / tl2::pi<t_real>;
Tobias WEBER's avatar
Tobias WEBER committed
881

Tobias WEBER's avatar
Tobias WEBER committed
882
#ifdef DEBUG
Tobias WEBER's avatar
Tobias WEBER committed
883
	std::cout << "a4_i = " << a4_i << ", a2_i = " << a2_i
Tobias WEBER's avatar
Tobias WEBER committed
884
		<< "; a4_f = " << a4_f << ", a2_f = " << a2_f
885
		<< "." << std::endl;
Tobias WEBER's avatar
Tobias WEBER committed
886
#endif
887

Tobias WEBER's avatar
Tobias WEBER committed
888
	// vertices in configuration space
Tobias WEBER's avatar
Tobias WEBER committed
889
890
	path.vec_i = AngleToPixel(a4_i, a2_i, true);
	path.vec_f = AngleToPixel(a4_f, a2_f, true);
Tobias WEBER's avatar
Tobias WEBER committed
891

Tobias WEBER's avatar
Tobias WEBER committed
892
#ifdef DEBUG
Tobias WEBER's avatar
Tobias WEBER committed
893
	std::cout << "start pixel: (" << path.vec_i[0] << ", " << path.vec_i[1] << std::endl;
Tobias WEBER's avatar
Tobias WEBER committed
894
	std::cout << "target pixel: (" << path.vec_f[0] << ", " << path.vec_f[1] << std::endl;
Tobias WEBER's avatar
Tobias WEBER committed
895
#endif
Tobias WEBER's avatar
Tobias WEBER committed
896
897

	// find closest voronoi vertices
Tobias WEBER's avatar
Tobias WEBER committed
898
	const auto& voro_vertices = m_voro_results.GetVoronoiVertices();
Tobias WEBER's avatar
Tobias WEBER committed
899
900
901
902
903

	// no voronoi vertices available
	if(voro_vertices.size() == 0)
		return path;

Tobias WEBER's avatar
Tobias WEBER committed
904
905
906

	std::size_t idx_i = 0;
	std::size_t idx_f = 0;
907

Tobias WEBER's avatar
Tobias WEBER committed
908
909
910
	// calculation of closest voronoi vertices using the index tree
	if(m_voro_results.GetIndexTreeSize())
	{
Tobias WEBER's avatar
Tobias WEBER committed
911
		std::vector<std::size_t> indices_i =
Tobias WEBER's avatar
Tobias WEBER committed
912
913
914
915
			m_voro_results.GetClosestVoronoiVertices(path.vec_i, 1);
		if(indices_i.size())
			idx_i = indices_i[0];

Tobias WEBER's avatar
Tobias WEBER committed
916
		std::vector<std::size_t> indices_f =
Tobias WEBER's avatar
Tobias WEBER committed
917
918
919
920
			m_voro_results.GetClosestVoronoiVertices(path.vec_f, 1);
		if(indices_f.size())
			idx_f = indices_f[0];
	}
921
922

	// alternate calculation without index tree
Tobias WEBER's avatar
Tobias WEBER committed
923
	else
Tobias WEBER's avatar
Tobias WEBER committed
924
	{
Tobias WEBER's avatar
Tobias WEBER committed
925
926
		t_real mindist_i = std::numeric_limits<t_real>::max();
		t_real mindist_f = std::numeric_limits<t_real>::max();
Tobias WEBER's avatar
Tobias WEBER committed
927

Tobias WEBER's avatar
Tobias WEBER committed
928
929
		for(std::size_t idx_vert = 0; idx_vert < voro_vertices.size(); ++idx_vert)
		{
Tobias WEBER's avatar
Tobias WEBER committed
930
			const t_vec2& cur_vert = voro_vertices[idx_vert];
Tobias WEBER's avatar
Tobias WEBER committed
931
			//std::cout << "cur_vert: " << cur_vert[0] << " " << cur_vert[1] << std::endl;
Tobias WEBER's avatar
Tobias WEBER committed
932

Tobias WEBER's avatar
Tobias WEBER committed
933
934
			t_vec2 diff_i = path.vec_i - cur_vert;
			t_vec2 diff_f = path.vec_f - cur_vert;
Tobias WEBER's avatar
Tobias WEBER committed
935

Tobias WEBER's avatar
Tobias WEBER committed
936
937
			t_real dist_i_sq = tl2::inner<t_vec2>(diff_i, diff_i);
			t_real dist_f_sq = tl2::inner<t_vec2>(diff_f, diff_f);
Tobias WEBER's avatar
Tobias WEBER committed
938

Tobias WEBER's avatar
Tobias WEBER committed
939
940
941
942
943
944
945
946
947
948
949
			if(dist_i_sq < mindist_i)
			{
				mindist_i = dist_i_sq;
				idx_i = idx_vert;
			}

			if(dist_f_sq < mindist_f)
			{
				mindist_f = dist_f_sq;
				idx_f = idx_vert;
			}
Tobias WEBER's avatar
Tobias WEBER committed
950
951
952
953
954
955
956
957
958
		}
	}

#ifdef DEBUG
	std::cout << "Nearest voronoi vertices: " << idx_i << ", " << idx_f << "." << std::endl;
#endif


	// find the shortest path between the voronoi vertices
Tobias WEBER's avatar
Tobias WEBER committed
959
	const auto& voro_graph = m_voro_results.GetVoronoiGraph();
Tobias WEBER's avatar
Tobias WEBER committed
960

Tobias WEBER's avatar
Tobias WEBER committed
961
962
963
964
	// are the graph vertex indices valid?
	if(idx_i >= voro_graph.GetNumVertices() || idx_f >= voro_graph.GetNumVertices())
		return path;

Tobias WEBER's avatar
Tobias WEBER committed
965
	const std::string& ident_i = voro_graph.GetVertexIdent(idx_i);
966
967
968
969
	using t_weight = typename t_graph::t_weight;


	// callback function with which the graph's edge weights can be modified
Tobias WEBER's avatar
Tobias WEBER committed
970
	auto weight_func = [this, &voro_graph, &voro_vertices, pathstragy](
971
972
973
974
975
976
977
		std::size_t idx1, std::size_t idx2) -> std::optional<t_weight>
	{
		// get original graph edge weight
		auto _weight = voro_graph.GetWeight(idx1, idx2);
		if(!_weight)
			return std::nullopt;

Tobias WEBER's avatar
Tobias WEBER committed
978
979
980
		// shortest path -> just use original edge weights
		if(pathstragy == PathStrategy::SHORTEST)
			return _weight;
981
982


Tobias WEBER's avatar
Tobias WEBER committed
983
984
985
		t_weight weight = *_weight;

		// get voronoi vertices of the current edge
986
987
988
		const t_vec2& vertex1 = voro_vertices[idx1];
		const t_vec2& vertex2 = voro_vertices[idx2];

989
990
991
		// get the distances to the wall vertices that are closest to the current voronoi vertices
		t_real dist1 = GetDistToNearestWall(vertex1);
		t_real dist2 = GetDistToNearestWall(vertex2);
Tobias WEBER's avatar
Tobias WEBER committed
992
		t_real min_dist = std::min(dist1, dist2);
Tobias WEBER's avatar
Tobias WEBER committed
993

Tobias WEBER's avatar
Tobias WEBER committed
994
		// modify edge weights using the minimum distance to the next wall
Tobias WEBER's avatar
Tobias WEBER committed
995
		if(pathstragy == PathStrategy::PENALISE_WALLS)
Tobias WEBER's avatar
Tobias WEBER committed
996
			return weight / min_dist;
997
998
999
1000

		return weight;
	};

1001

1002
	// find shortest path given the above weight function
Tobias WEBER's avatar
Tobias WEBER committed
1003
#if TASPATHS_SSSP_IMPL==1
1004
	const auto predecessors = geo::dijk(voro_graph, ident_i, &weight_func);
Tobias WEBER's avatar
Tobias WEBER committed
1005
#elif TASPATHS_SSSP_IMPL==2
1006
	const auto predecessors = geo::dijk_mod(voro_graph, ident_i, &weight_func);
Tobias WEBER's avatar
Tobias WEBER committed
1007
#elif TASPATHS_SSSP_IMPL==3
1008
1009
	const auto [distvecs, predecessors] = geo::bellman(
		voro_graph, ident_i, &weight_func);