polygons.cpp

/* 
 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Dmitry Logashenko
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */


#include "common/util/file_util.h"

/*
 * Implementation of the Point-in-Polygon algorithms
 */
#include "polygons.h"

namespace ug {
namespace d3f {

/**
 * Adds a point at the end of the list of corners of the polygon. Updates
 * the bounding box.
 */
void PolygonalRegions::t_polygon_item::append (number x, number y)
{
	MathVector<2> pnt (x, y);
	corner.push_back (pnt);
	
//	update the bounding box
	for (size_t i = 0; i < 2; i++)
		if (pnt[i] < bb[0][i]) bb[0][i] = pnt[i];
		else if (pnt[i] > bb[1][i]) bb[1][i] = pnt[i];
}

/**
 * Tests if a point is left or right of an infinite straight line.
 * Return true if pnt is left of the line or on the line, false otherwise.
 */
inline bool PolygonalRegions::t_polygon_item::pnt_is_on_the_left
(
	const MathVector<2> & pnt_0, ///< one point on the line
	const MathVector<2> & pnt_1, ///< another point on the line
	const MathVector<2> & pnt ///< point to check
) const
{
	return (pnt_1[0] - pnt_0[0]) * (pnt[1] - pnt_0[1])
			- (pnt[0] -  pnt_0[0]) * (pnt_1[1] - pnt_0[1]) >= 0;
}

/**
 * Checks whether a given point is in 'this' polygon.
 * The function returns false if no, true if yes.
 * This is an implementation of the winding-number algorithm taken from
 * http://geomalgorithms.com/a03-_inclusion.html
 * (Copyright 2000 softSurfer, 2012 Dan Sunday, free for public use)
 */
bool PolygonalRegions::t_polygon_item::pnt_inside_WN
(
	const MathVector<2> pnt ///< the point
) const
{
	/* Check whether we are in the bounding box: */
	if (pnt[0] < bb[0][0] || pnt[0] > bb[1][0]
		|| pnt[1] < bb[0][1] || pnt[1] > bb[1][1])
		return false; /* the point is definitely not in the polygon */
	
	int wn = 0; /* the winding number counter */
	
	/* Loop over all edges of the polygon: */
	for (size_t i = 0; i < corner.size (); i++)
	{
		size_t j = (i + 1) % corner.size ();
		/* Edge from corner[i] to corner[j] */
		if (corner[i][1] <= pnt[1])
		{
			if (corner[j][1] > pnt[1]) /* an upward crossing */
				if (pnt_is_on_the_left (corner[i], corner[j], pnt))
					++wn;
		}
		else
		{
			if (corner[j][1] <= pnt[1]) /* a downward crossing */
				if (! pnt_is_on_the_left (corner[i], corner[j], pnt))
					--wn;
		}
	}
	return wn != 0;
}

/**
 * Checks whether a given point is in a given element.
 * The function returns false if no, true if yes.
 * This is an implementation of the crossing-number algorithm taken from
 * http://geomalgorithms.com/a03-_inclusion.html
 * (Copyright 2000 softSurfer, 2012 Dan Sunday, free for public use)
 */
bool PolygonalRegions::t_polygon_item::pnt_inside_CN
(
	const MathVector<2> pnt ///< the point
) const
{
	/* Check whether we are in the bounding box: */
	if (pnt[0] < bb[0][0] || pnt[0] > bb[1][0]
		|| pnt[1] < bb[0][1] || pnt[1] > bb[1][1])
		return false; /* the point is definitely not in the polygon */
	
	int cn = 0; /* the crossing number counter */

	/* Loop over all edges of the polygon: */
	for (size_t i = 0; i < corner.size (); i++)
	{
		size_t j = (i + 1) % corner.size ();
		if ((corner[i][1] <= pnt[1] && corner[j][1] > pnt[1]) /* an upward crossing */
			|| (corner[i][1] > pnt[1] && corner[j][1] <= pnt[1])) /* a downward crossing */
		{
			/* Get the edge-ray intersection x-coordinate: */
			number t = (pnt[1]  - corner[i][1]) / (corner[j][1] - corner[i][1]);
			if (pnt[0] <  corner[i][0] + t * (corner[j][0] - corner[i][0]))
				++cn; /* a valid crossing by the ray */
		}
	}
	return (cn & 1) != 0; /* false if cn is even (out of the polygon), and true if cn is odd (in the polygon) */
}

/**
 * Loads polygon corners and id's from a given stream.
 */
void PolygonalRegions::load_points_from (std::istream & input)
{
	std::string input_line;
	std::istringstream line_stream;
	
	line_stream.exceptions (std::istream::failbit | std::istream::badbit);
	
	while (! input.eof ())
	{
		number x, y;
		int id;
		
		if (input.fail ())
			UG_THROW ("PolygonalRegions: Could not load the points from the file!");
		std::getline (input, input_line);
		if (input_line.length () == 0)
			continue;
		try
		{
			line_stream.str (input_line);
			line_stream.clear ();
			line_stream >> x >> y >> id;
		}
		catch (std::istream::failure & e)
		{
			UG_THROW ("PolygonalRegions: Failed to parse line '" << input_line << "' in the input file.");
		};
		
		append_point (x, y, id);
	}
	
	if (! good ())
		UG_THROW ("PolygonalRegions: Illegal specification of the polygons.");
}

/**
 * Loads polygon corners and id's from a given file.
 */
void PolygonalRegions::load_points_from (const char * file_name)
{
	std::string full_file_name = FindFileInStandardPaths (file_name);
	UG_COND_THROW (full_file_name.empty (), "PolygonalRegions: Coulnd't locate file '" << file_name << "'.");
	
	std::ifstream input (full_file_name.c_str (), std::ifstream::in);
	UG_COND_THROW (input.fail (), "PolygonalRegions: Cannot open data file '" << full_file_name << "' for input!");
	load_points_from (input);
}

} // namespace d3f
} // end namespace ug

/* End of File */