From bb03fa41d587efb913bd0959002007883af2d140 Mon Sep 17 00:00:00 2001
From: rinkk <karsten.rink@ufz.de>
Date: Thu, 12 Jan 2023 14:19:37 +0100
Subject: [PATCH] updated makebuildings to handle more complex geometries

---
 makeBuildings/makeBuildings.cpp | 433 +++++++++++++++++++++-----------
 1 file changed, 293 insertions(+), 140 deletions(-)

diff --git a/makeBuildings/makeBuildings.cpp b/makeBuildings/makeBuildings.cpp
index ef924ef..569d40d 100644
--- a/makeBuildings/makeBuildings.cpp
+++ b/makeBuildings/makeBuildings.cpp
@@ -1,15 +1,36 @@
-/**
- * @file   makeBildings.cpp
- * @author Karsten Rink
- * @date   2016/03/03
- * @brief  Takes polygons of building plans and creates simple 3d objects
- *
- * @copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/LICENSE.txt
- */
+//For this to also work with shared points in the original geometry, we need to remove the uniqe poins check in PointVec.cpp
+//for (std::size_t k(0); k < _data_vec->size(); ++k)
+//{
+//    GeoLib::Point* const pnt((*_data_vec)[k]);
+    /*
+    if (!_oct_tree->addPoint(pnt, ret_pnt))
+    {
+        assert(ret_pnt != nullptr);
+        _pnt_id_map[pnt->getID()] = ret_pnt->getID();
+        rm_pos.push_back(k);
+        delete (*_data_vec)[k];
+        (*_data_vec)[k] = nullptr;
+    }
+    else
+    {
+        _pnt_id_map[k] = pnt->getID();
+    }
+    */
+    //  _pnt_id_map[k] = pnt->getID();
+    //}
+
+    /**
+     * @file   makeBildings.cpp
+     * @author Karsten Rink
+     * @date   2016/03/03
+     * @brief  Takes polygons of building plans and creates simple 3d objects
+     *
+     * @copyright
+     * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+     *            Distributed under a Modified BSD License.
+     *              See accompanying file LICENSE.txt or
+     *              http://www.opengeosys.org/LICENSE.txt
+     */
 
 #include <algorithm>
 #include <memory>
@@ -17,148 +38,280 @@
 
 #include <tclap/CmdLine.h>
 
+#include "Applications/FileIO/Legacy/createSurface.h"
+
+#include "GeoLib/EarClippingTriangulation.h"
+
+#include "BaseLib/IO/readStringListFromFile.h"
 #include "GeoLib/GEOObjects.h"
+#include "GeoLib/Polygon.h"
+#include "GeoLib/Surface.h"
 #include "GeoLib/Triangle.h"
 #include "GeoLib/IO/XmlIO/Qt/XmlGmlInterface.h"
 
 #include <QCoreApplication>
 
-std::unique_ptr<std::vector<GeoLib::Polyline *>> copyPolylinesVector(
-	std::vector<GeoLib::Polyline *> const &polylines,
-	std::vector<GeoLib::Point *> const &points)
+
+std::vector<double> readHeightValuesFromFile(std::string const& filename)
 {
-	std::size_t n_plys = polylines.size();
-	std::unique_ptr<std::vector<GeoLib::Polyline *>> new_lines(new std::vector<GeoLib::Polyline *>(n_plys, nullptr));
-
-	for (std::size_t i = 0; i < n_plys; ++i)
-	{
-		if (polylines[i] == nullptr)
-			continue;
-		(*new_lines)[i] = new GeoLib::Polyline(points);
-		std::size_t nLinePnts(polylines[i]->getNumberOfPoints());
-		for (std::size_t j = 0; j < nLinePnts; ++j)
-			(*new_lines)[i]->addPoint(polylines[i]->getPointID(j));
-	}
-	return new_lines;
+    auto const str_vec = BaseLib::IO::readStringListFromFile(filename);
+    std::vector<double> dbl_vec;
+    std::transform(str_vec.begin(), str_vec.end(), std::back_inserter(dbl_vec),
+        [](const std::string& str) { return std::stod(str); });
+    return dbl_vec;
 }
 
-std::unique_ptr<std::vector<GeoLib::Surface *>> copySurfacesVector(
-	std::vector<GeoLib::Surface *> const &surfaces,
-	std::vector<GeoLib::Point *> const &points)
+std::unique_ptr<std::vector<GeoLib::Polyline*>> copyPolylinesVector(
+    std::vector<GeoLib::Polyline*> const& polylines,
+    std::vector<GeoLib::Point*> const& points)
 {
-	std::size_t n_sfc = surfaces.size();
-	std::unique_ptr<std::vector<GeoLib::Surface *>> new_surfaces(new std::vector<GeoLib::Surface *>(n_sfc, nullptr));
-
-	for (std::size_t i = 0; i < n_sfc; ++i)
-	{
-		if (surfaces[i] == nullptr)
-			continue;
-		(*new_surfaces)[i] = new GeoLib::Surface(points);
-
-		std::size_t n_tris(surfaces[i]->getNumberOfTriangles());
-		for (std::size_t j = 0; j < n_tris; ++j)
-		{
-			GeoLib::Triangle const *t = (*surfaces[i])[j];
-			(*new_surfaces)[i]->addTriangle(t->getPoint(0)->getID(), t->getPoint(1)->getID(), t->getPoint(2)->getID());
-		}
-	}
-	return new_surfaces;
+    std::size_t n_plys = polylines.size();
+    std::unique_ptr<std::vector<GeoLib::Polyline*>> new_lines(
+        new std::vector<GeoLib::Polyline*>(n_plys, nullptr));
+
+    for (std::size_t i = 0; i < n_plys; ++i)
+    {
+        if (polylines[i] == nullptr)
+            continue;
+        (*new_lines)[i] = new GeoLib::Polyline(points);
+        std::size_t nLinePnts(polylines[i]->getNumberOfPoints());
+        for (std::size_t j = 0; j < nLinePnts/*-1*/; ++j)
+            (*new_lines)[i]->addPoint(polylines[i]->getPointID(j));
+        //(*new_lines)[i]->addPoint(polylines[i]->getPointID(0));
+    }
+    return new_lines;
 }
 
-void makeBuildings(GeoLib::GEOObjects &geo_objects, std::string const &geo_name, std::string &output_name, double height)
+std::unique_ptr<std::vector<GeoLib::Surface*>> copySurfacesVector(
+    std::vector<GeoLib::Surface*> const& surfaces,
+    std::vector<GeoLib::Point*> const& points)
 {
-	std::vector<GeoLib::Point *> const *pnts(geo_objects.getPointVec(geo_name));
-	std::vector<GeoLib::Polyline *> const *plys(geo_objects.getPolylineVec(geo_name));
-	std::vector<GeoLib::Surface *> const *sfcs(geo_objects.getSurfaceVec(geo_name));
-	std::size_t const n_pnts(pnts->size());
-	std::size_t const n_plys(plys->size());
-	std::size_t const n_sfcs(sfcs->size());
-
-	std::unique_ptr<std::vector<GeoLib::Point *>> new_pnts(new std::vector<GeoLib::Point *>);
-	for (GeoLib::Point *point : *pnts)
-		if (point)
-			new_pnts->push_back(new GeoLib::Point(*point));
-
-	std::unique_ptr<std::vector<GeoLib::Polyline *>> new_plys = copyPolylinesVector(*plys, *new_pnts);
-	std::unique_ptr<std::vector<GeoLib::Surface *>> new_sfcs = copySurfacesVector(*sfcs, *new_pnts);
-
-	for (GeoLib::Point *point : *pnts)
-	{
-		if (point)
-			new_pnts->push_back(new GeoLib::Point((*point)[0], (*point)[1], (*point)[2] + height, point->getID() + n_pnts));
-		else
-			new_pnts->push_back(nullptr);
-	}
-
-	for (GeoLib::Polyline *p : *plys)
-	{
-		if (p == nullptr)
-			continue;
-		std::size_t const np(p->getNumberOfPoints());
-		GeoLib::Surface *s = new GeoLib::Surface(*new_pnts);
-		for (std::size_t i = 1; i < np; ++i)
-		{
-			s->addTriangle(p->getPoint(i)->getID(), p->getPoint(i - 1)->getID(), p->getPoint(i - 1)->getID() + n_pnts);
-			s->addTriangle(p->getPoint(i)->getID(), p->getPoint(i - 1)->getID() + n_pnts, p->getPoint(i)->getID() + n_pnts);
-		}
-		new_sfcs->push_back(s);
-	}
-
-	for (std::size_t j = 0; j < n_sfcs; ++j)
-	{
-		if ((*sfcs)[j] == nullptr)
-			continue;
-		std::size_t const n_tris((*sfcs)[j]->getNumberOfTriangles());
-		GeoLib::Surface *s = new GeoLib::Surface(*new_pnts);
-		for (std::size_t i = 0; i < n_tris; i++)
-		{
-			GeoLib::Triangle const *t = (*(*sfcs)[j])[i];
-			s->addTriangle(t->getPoint(0)->getID() + n_pnts, t->getPoint(1)->getID() + n_pnts, t->getPoint(2)->getID() + n_pnts);
-		}
-		new_sfcs->push_back(s);
-	}
-
-	geo_objects.addPointVec(std::move(new_pnts), output_name);
-	geo_objects.addPolylineVec(std::move(new_plys), output_name);
-	geo_objects.addSurfaceVec(std::move(new_sfcs), output_name);
+    std::size_t n_sfc = surfaces.size();
+    std::unique_ptr<std::vector<GeoLib::Surface*>> new_surfaces(
+        new std::vector<GeoLib::Surface*>(n_sfc, nullptr));
+
+    for (std::size_t i = 0; i < n_sfc; ++i)
+    {
+        if (surfaces[i] == nullptr)
+            continue;
+        (*new_surfaces)[i] = new GeoLib::Surface(points);
+
+        std::size_t n_tris(surfaces[i]->getNumberOfTriangles());
+        for (std::size_t j = 0; j < n_tris; ++j)
+        {
+            GeoLib::Triangle const* t = (*surfaces[i])[j];
+            (*new_surfaces)[i]->addTriangle(t->getPoint(0)->getID(),
+                t->getPoint(1)->getID(),
+                t->getPoint(2)->getID());
+        }
+    }
+    return new_surfaces;
 }
 
-int main(int argc, char *argv[])
+void makeBuildings(GeoLib::GEOObjects& geo_objects, std::string const& geo_name,
+    std::string& output_name, double const height,
+    std::vector<double> const& height_vec)
 {
-	QCoreApplication app(argc, argv);
-
-	TCLAP::CmdLine cmd("Uses polygons from building plans to create 3d objects.", ' ', "0.1");
-
-	TCLAP::ValueArg<double> height("s", "size", "height of the 3d objects (buildings) in metres", true,
-								   1.0, "height of objects");
-	cmd.add(height);
-	TCLAP::ValueArg<std::string> geo_out("o", "geo-output-file",
-										 "the name of the file the 3d geometry will be written to", true,
-										 "", "file name of output geometry");
-	cmd.add(geo_out);
-	TCLAP::ValueArg<std::string> geo_in("i", "geo-input-file",
-										"the name of the file containing the input geometry", true,
-										"", "file name of input geometry");
-	cmd.add(geo_in);
-
-	cmd.parse(argc, argv);
-
-	INFO("Reading geometry {:s}.", geo_in.getValue().c_str());
-
-	GeoLib::GEOObjects geo_objects;
-	GeoLib::IO::XmlGmlInterface xml(geo_objects);
-	if (!xml.readFile(geo_in.getValue()))
-	{
-		ERR("Error reading geometry.");
-		return 1;
-	}
-	auto const geo_names = geo_objects.getGeometryNames();
-	std::string output_name("output");
-
-	makeBuildings(geo_objects, geo_names[0], output_name, height.getValue());
-
-	xml.export_name = output_name;
-	BaseLib::IO::writeStringToFile(xml.writeToString(), geo_out.getValue());
-
-	return 0;
+    bool const set_height = (height_vec.empty()) ? false : true;
+
+    auto pnts = geo_objects.getPointVec(geo_name);
+    std::size_t const n_pnts(pnts->size());
+    std::unique_ptr<std::vector<GeoLib::Point*>> new_pnts(
+        new std::vector<GeoLib::Point*>);
+    for (GeoLib::Point* point : *pnts)
+    {
+        if (point != nullptr)
+        {
+            new_pnts->push_back(new GeoLib::Point(*point));
+        }
+    }
+
+    auto plys = geo_objects.getPolylineVec(geo_name);
+    std::size_t const n_plys(plys->size());
+    std::unique_ptr<std::vector<GeoLib::Polyline*>> new_plys =
+        copyPolylinesVector(*plys, *new_pnts);
+
+    auto sfcs = const_cast<std::vector<GeoLib::Surface*>*>(
+        geo_objects.getSurfaceVec(geo_name));
+    if (sfcs == nullptr)
+    {
+        sfcs = new std::vector<GeoLib::Surface*>;
+        for (auto line : *new_plys)
+        {
+            auto sfc = FileIO::createSurfaceWithEarClipping(*line);
+            if (sfc != nullptr)
+            {
+
+                sfcs->push_back(sfc.release());
+            }
+        }
+    }
+    std::size_t const n_sfcs = sfcs->size();
+    std::unique_ptr<std::vector<GeoLib::Surface*>> new_sfcs =
+        copySurfacesVector(*sfcs, *new_pnts);
+
+    for (std::size_t ply_id = 0; ply_id < n_plys; ++ply_id)
+    {
+        GeoLib::Polyline* p = (*new_plys)[ply_id];
+        if (p == nullptr)
+        {
+            continue;
+        }
+
+        std::size_t const np(p->getNumberOfPoints());
+        GeoLib::Surface* s = new GeoLib::Surface(*new_pnts);
+        /**/
+        GeoLib::Point const* pnt0 = p->getPoint(0);
+        new_pnts->push_back(new GeoLib::Point((*pnt0)[0], (*pnt0)[1],
+            (*pnt0)[2] + height_vec[ply_id],
+            new_pnts->size()));
+        GeoLib::Polyline tmp(*new_pnts);
+        tmp.addPoint(new_pnts->size() - 1);
+        /**/
+        for (std::size_t i = 1; i < np - 1; ++i)
+        {
+            /*
+            std::size_t const pnt_idx1 = p->getPoint(i - 1)->getID();
+            std::size_t const pnt_idx2 = p->getPoint(i)->getID();
+            if (set_height)
+            {
+                (*(*new_pnts)[pnt_idx1 + n_pnts])[2] =
+                    (*(*pnts)[pnt_idx1])[2] + height_vec[ply_id];
+                (*(*new_pnts)[pnt_idx2 + n_pnts])[2] =
+                    (*(*pnts)[pnt_idx2])[2] + height_vec[ply_id];
+            }
+
+            s->addTriangle(p->getPoint(i)->getID(), p->getPoint(i - 1)->getID(),
+                           p->getPoint(i - 1)->getID() + n_pnts);
+            s->addTriangle(p->getPoint(i)->getID(),
+                           p->getPoint(i - 1)->getID() + n_pnts,
+                           p->getPoint(i)->getID() + n_pnts);
+            */
+            GeoLib::Point const* pnti = p->getPoint(i);
+            new_pnts->push_back(new GeoLib::Point(
+                (*pnti)[0], (*pnti)[1], (*pnti)[2] + height_vec[ply_id],
+                new_pnts->size()));
+            tmp.addPoint(new_pnts->size() - 1);
+
+            s->addTriangle(p->getPoint(i)->getID(), p->getPoint(i - 1)->getID(),
+                new_pnts->size() - 2);
+            s->addTriangle(p->getPoint(i)->getID(),
+                new_pnts->size() - 2,
+                new_pnts->size() - 1);
+        }
+        tmp.addPoint(tmp.getPointID(0));
+
+        s->addTriangle(p->getPoint(np - 2)->getID(), p->getPoint(0)->getID(),
+            new_pnts->size() - 1);
+        s->addTriangle(p->getPoint(0)->getID(),
+            new_pnts->size() - 1,
+            tmp.getPointID(0));
+        new_sfcs->push_back(s);
+
+        auto new_sfc = FileIO::createSurfaceWithEarClipping(tmp);
+        new_sfcs->push_back(new_sfc.release());
+    }
+    /*
+    for (std::size_t j = 0; j < n_sfcs; ++j)
+    {
+        if ((*sfcs)[j] == nullptr)
+            continue;
+        std::size_t const n_tris((*sfcs)[j]->getNumberOfTriangles());
+        GeoLib::Surface* s = new GeoLib::Surface(*new_pnts);
+        for (std::size_t i = 0; i < n_tris; i++)
+        {
+            GeoLib::Triangle const* t = (*(*sfcs)[j])[i];
+            s->addTriangle(t->getPoint(0)->getID() + n_pnts,
+                           t->getPoint(1)->getID() + n_pnts,
+                           t->getPoint(2)->getID() + n_pnts);
+        }
+        new_sfcs->push_back(s);
+    }
+    */
+    geo_objects.addPointVec(std::move(*new_pnts), output_name);
+    geo_objects.addPolylineVec(std::move(*new_plys), output_name,
+        GeoLib::PolylineVec::NameIdMap{});
+    geo_objects.addSurfaceVec(std::move(*new_sfcs), output_name,
+        GeoLib::SurfaceVec::NameIdMap{});
+}
+
+int main(int argc, char* argv[])
+{
+    QCoreApplication app(argc, argv);
+
+    TCLAP::CmdLine cmd(
+        "Uses polygons from building plans to create 3d objects.", ' ', "0.1");
+
+    TCLAP::ValueArg<std::string> file_height_arg(
+        "", "FileSize",
+        "file containing the height for each polygone in metres", false, "",
+        "height of objects");
+    cmd.add(file_height_arg);
+    TCLAP::ValueArg<double> fixed_height_arg(
+        "s", "FixedSize",
+        "fixed height of the 3d objects (buildings) in metres", false, 0,
+        "height of objects");
+    cmd.add(fixed_height_arg);
+    TCLAP::ValueArg<std::string> geo_out(
+        "o", "geo-output-file",
+        "the name of the file the 3d geometry will be written to", true, "",
+        "file name of output geometry");
+    cmd.add(geo_out);
+    TCLAP::ValueArg<std::string> geo_in(
+        "i", "geo-input-file",
+        "the name of the file containing the input geometry", true, "",
+        "file name of input geometry");
+    cmd.add(geo_in);
+
+    cmd.parse(argc, argv);
+
+    if (!fixed_height_arg.isSet() && !file_height_arg.isSet())
+    {
+        ERR("Either a fixed height value or a file containing variable heights "
+            "must be specified.");
+        return 3;
+    }
+
+    if (fixed_height_arg.isSet() && file_height_arg.isSet())
+    {
+        WARN(
+            "Building height will be based on given file values. Fixed height "
+            "value will be ignored.");
+    }
+
+    INFO("Reading geometry {:s}.", geo_in.getValue().c_str());
+    GeoLib::GEOObjects geo_objects;
+    GeoLib::IO::XmlGmlInterface xml(geo_objects);
+    if (!xml.readFile(geo_in.getValue()))
+    {
+        ERR("Error reading geometry.");
+        return 1;
+    }
+    std::vector<std::string> geo_names = geo_objects.getGeometryNames();
+
+    std::vector<double> height_vec;
+    if (file_height_arg.isSet())
+    {
+        height_vec = readHeightValuesFromFile(file_height_arg.getValue());
+
+        if (geo_objects.getPolylineVec(geo_names[0])->size() !=
+            height_vec.size())
+        {
+            ERR("Length of polyline vector ({:d}) is not equal to length of "
+                "height vector ({:d}).",
+                geo_objects.getPolylineVec(geo_names[0])->size(),
+                height_vec.size());
+            return 2;
+        }
+    }
+
+    std::string output_name("output");
+    makeBuildings(geo_objects, geo_names[0], output_name,
+        fixed_height_arg.getValue(), height_vec);
+
+    xml.export_name = output_name;
+    // xml.writeToFile(geo_out.getValue());
+    BaseLib::IO::writeStringToFile(xml.writeToString(), geo_out.getValue());
+
+    return 0;
 }
-- 
GitLab