<commit_before>#include <limits>
#include "Skeletonize.h"
#include <vigra/multi_distance.hxx>
#include <vigra/multi_gridgraph.hxx>
#include <vigra/multi_labeling.hxx>
#include <vigra/multi_impex.hxx>
#include <util/timing.h>
#include <util/ProgramOptions.h>
#include <util/Logger.h>

util::ProgramOption optionSkeletonBoundaryWeight(
		util::_long_name        = "skeletonBoundaryWeight",
		util::_description_text = "The weight of the boundary term to find the tube's skeletons.",
		util::_default_value    = 1);

util::ProgramOption optionSkeletonMaxNumSegments(
		util::_long_name        = "skeletonMaxNumSegments",
		util::_description_text = "The maximal number of segments to extract for a skeleton.",
		util::_default_value    = 10);

util::ProgramOption optionSkeletonMinSegmentLength(
		util::_long_name        = "skeletonMinSegmentLength",
		util::_description_text = "The mininal length of a segment (including the boundary penalty) to extract for a skeleton.",
		util::_default_value    = 0);

util::ProgramOption optionSkeletonMinSegmentLengthRatio(
		util::_long_name        = "skeletonMinSegmentLengthRatio",
		util::_description_text = "The mininal length of a segment (including the boundary penalty) as a ration of the largest segment to extract for a skeleton.",
		util::_default_value    = 1);

util::ProgramOption optionSkeletonSkipExplainedNodes(
		util::_long_name        = "skeletonSkipExplainedNodes",
		util::_description_text = "Don't add segments to nodes that are already explained by the current skeleton. "
		                          "Nodes are explained, if they fall within a sphere around any current skeleton node. "
		                          "The size of the sphere is determined by boundary distance * skeletonExplanationWeight.");

util::ProgramOption optionSkeletonExplanationWeight(
		util::_long_name        = "skeletonExplanationWeight",
		util::_description_text = "A factor to multiply with the boundary distance to create 'explanation spheres'. "
		                          "See skeletonSkipExplainedNodes.",
		util::_default_value    = 1);

util::ProgramOption optionSkeletonDownsampleVolume(
		util::_long_name        = "skeletonDownsampleVolume",
		util::_description_text = "downsample the volume dimensions by the largest power of two that does not change connectivity.");

logger::LogChannel skeletonizelog("skeletonizelog", "[Skeletonize] ");

Skeletonize::Skeletonize(ExplicitVolume<unsigned char>& volume) :
	_downSampleVolume(optionSkeletonDownsampleVolume),
	_positionMap(_graph),
	_distanceMap(_graph),
	_boundaryWeight(optionSkeletonBoundaryWeight),
	_dijkstra(_graph, _distanceMap),
	_minSegmentLength(optionSkeletonMinSegmentLength),
	_minSegmentLengthRatio(optionSkeletonMinSegmentLengthRatio),
	_skipExplainedNodes(optionSkeletonSkipExplainedNodes),
	_explanationWeight(optionSkeletonExplanationWeight) {

	Timer t(__FUNCTION__);

	setVolume(volume);
	createVoxelGraph();
}

void
Skeletonize::setVolume(const ExplicitVolume<unsigned char>& volume) {

	if (_downSampleVolume)
		downsampleVolume(volume);
	else
		_volume = volume;

	_boundaryDistance = vigra::MultiArray<3, float>(_volume.data().shape(), 0.0);
}

void
Skeletonize::downsampleVolume(const ExplicitVolume<unsigned char>& volume) {

	vigra::TinyVector<float, 3> origRes = {
			volume.getResolutionX(),
			volume.getResolutionY(),
			volume.getResolutionZ()};


	vigra::TinyVector<int, 3> origSize = {
			(int)volume.width(),
			(int)volume.height(),
			(int)volume.depth()};

	float finestRes = -1;
	int   finestDimension;
	for (int d = 0; d < 3; d++)
		if (finestRes < 0 || finestRes > origRes[d]) {

			finestRes = origRes[d];
			finestDimension = d;
		}

	// the largest downsample factor to consider
	int downsampleFactor = 8;

	bool goodLevelFound = false;
	while (!goodLevelFound) {

		LOG_DEBUG(skeletonizelog)
				<< "trying to downsample finest dimension by factor "
				<< downsampleFactor << std::endl;

		vigra::TinyVector<int, 3>   factors;
		vigra::TinyVector<float, 3> targetRes;
		vigra::TinyVector<int, 3>   targetSize;

		factors[finestDimension]    = downsampleFactor;
		targetRes[finestDimension]  = origRes[finestDimension]*downsampleFactor;
		targetSize[finestDimension] = origSize[finestDimension]/downsampleFactor;

		// the target resolution of the finest dimension, when downsampled with 
		// current factor
		float targetFinestRes = finestRes*downsampleFactor;

		// for each other dimension, find best downsample factor
		for (int d = 0; d < 3; d++) {

			if (d == finestDimension)
				continue;

			int bestFactor = 0;
			float minResDiff = 0;

			for (int f = downsampleFactor; f != 0; f /= 2) {

				float targetRes = origRes[d]*f;
				float resDiff = std::abs(targetFinestRes - targetRes);

				if (bestFactor == 0 || resDiff < minResDiff) {

					bestFactor = f;
					minResDiff = resDiff;
				}
			}

			factors[d]    = bestFactor;
			targetRes[d]  = origRes[d]*bestFactor;
			targetSize[d] = origSize[d]/bestFactor;
		}

		LOG_DEBUG(skeletonizelog)
				<< "best downsampling factors for each dimension are "
				<< factors << std::endl;

		_volume = ExplicitVolume<unsigned char>(targetSize[0], targetSize[1], targetSize[2]);
		_volume.setResolution(targetRes[0], targetRes[1], targetRes[2]);
		_volume.setBoundingBox(volume.getBoundingBox());

		// copy volume
		for (int z = 0; z < targetSize[2]; z++)
		for (int y = 0; y < targetSize[1]; y++)
		for (int x = 0; x < targetSize[0]; x++)
			_volume(x, y, z) = volume(x*factors[0], y*factors[1], z*factors[2]);

		vigra::exportVolume(_volume.data(), vigra::VolumeExportInfo("downsampled/downsampled", ".tif"));
		//UTIL_THROW_EXCEPTION(Exception, "breakpoint");

		int numRegions;
		try {

			// check for downsampling errors
			vigra::MultiArray<3, unsigned int> labels(_volume.data().shape());
			numRegions = vigra::labelMultiArrayWithBackground(
					_volume.data(),
					labels);
<commit_msg>removed debugging file output in Skeletonize<commit_after>#include <limits>
#include "Skeletonize.h"
#include <vigra/multi_distance.hxx>
#include <vigra/multi_gridgraph.hxx>
#include <vigra/multi_labeling.hxx>
#include <vigra/multi_impex.hxx>
#include <util/timing.h>
#include <util/ProgramOptions.h>
#include <util/Logger.h>

util::ProgramOption optionSkeletonBoundaryWeight(
		util::_long_name        = "skeletonBoundaryWeight",
		util::_description_text = "The weight of the boundary term to find the tube's skeletons.",
		util::_default_value    = 1);

util::ProgramOption optionSkeletonMaxNumSegments(
		util::_long_name        = "skeletonMaxNumSegments",
		util::_description_text = "The maximal number of segments to extract for a skeleton.",
		util::_default_value    = 10);

util::ProgramOption optionSkeletonMinSegmentLength(
		util::_long_name        = "skeletonMinSegmentLength",
		util::_description_text = "The mininal length of a segment (including the boundary penalty) to extract for a skeleton.",
		util::_default_value    = 0);

util::ProgramOption optionSkeletonMinSegmentLengthRatio(
		util::_long_name        = "skeletonMinSegmentLengthRatio",
		util::_description_text = "The mininal length of a segment (including the boundary penalty) as a ration of the largest segment to extract for a skeleton.",
		util::_default_value    = 1);

util::ProgramOption optionSkeletonSkipExplainedNodes(
		util::_long_name        = "skeletonSkipExplainedNodes",
		util::_description_text = "Don't add segments to nodes that are already explained by the current skeleton. "
		                          "Nodes are explained, if they fall within a sphere around any current skeleton node. "
		                          "The size of the sphere is determined by boundary distance * skeletonExplanationWeight.");

util::ProgramOption optionSkeletonExplanationWeight(
		util::_long_name        = "skeletonExplanationWeight",
		util::_description_text = "A factor to multiply with the boundary distance to create 'explanation spheres'. "
		                          "See skeletonSkipExplainedNodes.",
		util::_default_value    = 1);

util::ProgramOption optionSkeletonDownsampleVolume(
		util::_long_name        = "skeletonDownsampleVolume",
		util::_description_text = "downsample the volume dimensions by the largest power of two that does not change connectivity.");

logger::LogChannel skeletonizelog("skeletonizelog", "[Skeletonize] ");

Skeletonize::Skeletonize(ExplicitVolume<unsigned char>& volume) :
	_downSampleVolume(optionSkeletonDownsampleVolume),
	_positionMap(_graph),
	_distanceMap(_graph),
	_boundaryWeight(optionSkeletonBoundaryWeight),
	_dijkstra(_graph, _distanceMap),
	_minSegmentLength(optionSkeletonMinSegmentLength),
	_minSegmentLengthRatio(optionSkeletonMinSegmentLengthRatio),
	_skipExplainedNodes(optionSkeletonSkipExplainedNodes),
	_explanationWeight(optionSkeletonExplanationWeight) {

	Timer t(__FUNCTION__);

	setVolume(volume);
	createVoxelGraph();
}

void
Skeletonize::setVolume(const ExplicitVolume<unsigned char>& volume) {

	if (_downSampleVolume)
		downsampleVolume(volume);
	else
		_volume = volume;

	_boundaryDistance = vigra::MultiArray<3, float>(_volume.data().shape(), 0.0);
}

void
Skeletonize::downsampleVolume(const ExplicitVolume<unsigned char>& volume) {

	vigra::TinyVector<float, 3> origRes = {
			volume.getResolutionX(),
			volume.getResolutionY(),
			volume.getResolutionZ()};


	vigra::TinyVector<int, 3> origSize = {
			(int)volume.width(),
			(int)volume.height(),
			(int)volume.depth()};

	float finestRes = -1;
	int   finestDimension;
	for (int d = 0; d < 3; d++)
		if (finestRes < 0 || finestRes > origRes[d]) {

			finestRes = origRes[d];
			finestDimension = d;
		}

	// the largest downsample factor to consider
	int downsampleFactor = 8;

	bool goodLevelFound = false;
	while (!goodLevelFound) {

		LOG_DEBUG(skeletonizelog)
				<< "trying to downsample finest dimension by factor "
				<< downsampleFactor << std::endl;

		vigra::TinyVector<int, 3>   factors;
		vigra::TinyVector<float, 3> targetRes;
		vigra::TinyVector<int, 3>   targetSize;

		factors[finestDimension]    = downsampleFactor;
		targetRes[finestDimension]  = origRes[finestDimension]*downsampleFactor;
		targetSize[finestDimension] = origSize[finestDimension]/downsampleFactor;

		// the target resolution of the finest dimension, when downsampled with 
		// current factor
		float targetFinestRes = finestRes*downsampleFactor;

		// for each other dimension, find best downsample factor
		for (int d = 0; d < 3; d++) {

			if (d == finestDimension)
				continue;

			int bestFactor = 0;
			float minResDiff = 0;

			for (int f = downsampleFactor; f != 0; f /= 2) {

				float targetRes = origRes[d]*f;
				float resDiff = std::abs(targetFinestRes - targetRes);

				if (bestFactor == 0 || resDiff < minResDiff) {

					bestFactor = f;
					minResDiff = resDiff;
				}
			}

			factors[d]    = bestFactor;
			targetRes[d]  = origRes[d]*bestFactor;
			targetSize[d] = origSize[d]/bestFactor;
		}

		LOG_DEBUG(skeletonizelog)
				<< "best downsampling factors for each dimension are "
				<< factors << std::endl;

		_volume = ExplicitVolume<unsigned char>(targetSize[0], targetSize[1], targetSize[2]);
		_volume.setResolution(targetRes[0], targetRes[1], targetRes[2]);
		_volume.setBoundingBox(volume.getBoundingBox());

		// copy volume
		for (int z = 0; z < targetSize[2]; z++)
		for (int y = 0; y < targetSize[1]; y++)
		for (int x = 0; x < targetSize[0]; x++)
			_volume(x, y, z) = volume(x*factors[0], y*factors[1], z*factors[2]);

		int numRegions;
		try {

			// check for downsampling errors
			vigra::MultiArray<3, unsigned int> labels(_volume.data().shape());
			numRegions = vigra::labelMultiArrayWithBackground(
					_volume.data(),
					labels);

{'fraction_non_alphanumeric': 0.11249786288254403, 'fraction_numerical': 0.007180714652077278, 'mean_word_length': 4.643511818620357, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2796320', 'n_tokens_mistral': 4036, 'n_tokens_neox': 3597, 'n_words': 876}