diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,3 @@ +0.1.0.0 +------- +* initial release diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Stéphane Laurent (c) 2022 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Author name here nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/LICENSE_CDT.txt b/LICENSE_CDT.txt new file mode 100644 --- /dev/null +++ b/LICENSE_CDT.txt @@ -0,0 +1,384 @@ +The 'hcdt' package as a whole is distributed under BSD-3-Clause. + +It includes an open source C++ library: 'CDT'. + +Full copy of the license agreement used by this component is included +below. + + + +### 'CDT' license + +Mozilla Public License Version 2.0 +================================== + +1. Definitions +-------------- + +1.1. "Contributor" + means each individual or legal entity that creates, contributes to + the creation of, or owns Covered Software. + +1.2. "Contributor Version" + means the combination of the Contributions of others (if any) used + by a Contributor and that particular Contributor's Contribution. + +1.3. "Contribution" + means Covered Software of a particular Contributor. + +1.4. "Covered Software" + means Source Code Form to which the initial Contributor has attached + the notice in Exhibit A, the Executable Form of such Source Code + Form, and Modifications of such Source Code Form, in each case + including portions thereof. + +1.5. "Incompatible With Secondary Licenses" + means + + (a) that the initial Contributor has attached the notice described + in Exhibit B to the Covered Software; or + + (b) that the Covered Software was made available under the terms of + version 1.1 or earlier of the License, but not also under the + terms of a Secondary License. + +1.6. "Executable Form" + means any form of the work other than Source Code Form. + +1.7. "Larger Work" + means a work that combines Covered Software with other material, in + a separate file or files, that is not Covered Software. + +1.8. "License" + means this document. + +1.9. "Licensable" + means having the right to grant, to the maximum extent possible, + whether at the time of the initial grant or subsequently, any and + all of the rights conveyed by this License. + +1.10. "Modifications" + means any of the following: + + (a) any file in Source Code Form that results from an addition to, + deletion from, or modification of the contents of Covered + Software; or + + (b) any new file in Source Code Form that contains any Covered + Software. + +1.11. "Patent Claims" of a Contributor + means any patent claim(s), including without limitation, method, + process, and apparatus claims, in any patent Licensable by such + Contributor that would be infringed, but for the grant of the + License, by the making, using, selling, offering for sale, having + made, import, or transfer of either its Contributions or its + Contributor Version. + +1.12. "Secondary License" + means either the GNU General Public License, Version 2.0, the GNU + Lesser General Public License, Version 2.1, the GNU Affero General + Public License, Version 3.0, or any later versions of those + licenses. + +1.13. "Source Code Form" + means the form of the work preferred for making modifications. + +1.14. "You" (or "Your") + means an individual or a legal entity exercising rights under this + License. For legal entities, "You" includes any entity that + controls, is controlled by, or is under common control with You. For + purposes of this definition, "control" means (a) the power, direct + or indirect, to cause the direction or management of such entity, + whether by contract or otherwise, or (b) ownership of more than + fifty percent (50%) of the outstanding shares or beneficial + ownership of such entity. + +2. License Grants and Conditions +-------------------------------- + +2.1. Grants + +Each Contributor hereby grants You a world-wide, royalty-free, +non-exclusive license: + +(a) under intellectual property rights (other than patent or trademark) + Licensable by such Contributor to use, reproduce, make available, + modify, display, perform, distribute, and otherwise exploit its + Contributions, either on an unmodified basis, with Modifications, or + as part of a Larger Work; and + +(b) under Patent Claims of such Contributor to make, use, sell, offer + for sale, have made, import, and otherwise transfer either its + Contributions or its Contributor Version. + +2.2. Effective Date + +The licenses granted in Section 2.1 with respect to any Contribution +become effective for each Contribution on the date the Contributor first +distributes such Contribution. + +2.3. Limitations on Grant Scope + +The licenses granted in this Section 2 are the only rights granted under +this License. No additional rights or licenses will be implied from the +distribution or licensing of Covered Software under this License. +Notwithstanding Section 2.1(b) above, no patent license is granted by a +Contributor: + +(a) for any code that a Contributor has removed from Covered Software; + or + +(b) for infringements caused by: (i) Your and any other third party's + modifications of Covered Software, or (ii) the combination of its + Contributions with other software (except as part of its Contributor + Version); or + +(c) under Patent Claims infringed by Covered Software in the absence of + its Contributions. + +This License does not grant any rights in the trademarks, service marks, +or logos of any Contributor (except as may be necessary to comply with +the notice requirements in Section 3.4). + +2.4. Subsequent Licenses + +No Contributor makes additional grants as a result of Your choice to +distribute the Covered Software under a subsequent version of this +License (see Section 10.2) or under the terms of a Secondary License (if +permitted under the terms of Section 3.3). + +2.5. Representation + +Each Contributor represents that the Contributor believes its +Contributions are its original creation(s) or it has sufficient rights +to grant the rights to its Contributions conveyed by this License. + +2.6. Fair Use + +This License is not intended to limit any rights You have under +applicable copyright doctrines of fair use, fair dealing, or other +equivalents. + +2.7. Conditions + +Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted +in Section 2.1. + +3. Responsibilities +------------------- + +3.1. Distribution of Source Form + +All distribution of Covered Software in Source Code Form, including any +Modifications that You create or to which You contribute, must be under +the terms of this License. You must inform recipients that the Source +Code Form of the Covered Software is governed by the terms of this +License, and how they can obtain a copy of this License. You may not +attempt to alter or restrict the recipients' rights in the Source Code +Form. + +3.2. Distribution of Executable Form + +If You distribute Covered Software in Executable Form then: + +(a) such Covered Software must also be made available in Source Code + Form, as described in Section 3.1, and You must inform recipients of + the Executable Form how they can obtain a copy of such Source Code + Form by reasonable means in a timely manner, at a charge no more + than the cost of distribution to the recipient; and + +(b) You may distribute such Executable Form under the terms of this + License, or sublicense it under different terms, provided that the + license for the Executable Form does not attempt to limit or alter + the recipients' rights in the Source Code Form under this License. + +3.3. Distribution of a Larger Work + +You may create and distribute a Larger Work under terms of Your choice, +provided that You also comply with the requirements of this License for +the Covered Software. If the Larger Work is a combination of Covered +Software with a work governed by one or more Secondary Licenses, and the +Covered Software is not Incompatible With Secondary Licenses, this +License permits You to additionally distribute such Covered Software +under the terms of such Secondary License(s), so that the recipient of +the Larger Work may, at their option, further distribute the Covered +Software under the terms of either this License or such Secondary +License(s). + +3.4. Notices + +You may not remove or alter the substance of any license notices +(including copyright notices, patent notices, disclaimers of warranty, +or limitations of liability) contained within the Source Code Form of +the Covered Software, except that You may alter any license notices to +the extent required to remedy known factual inaccuracies. + +3.5. Application of Additional Terms + +You may choose to offer, and to charge a fee for, warranty, support, +indemnity or liability obligations to one or more recipients of Covered +Software. However, You may do so only on Your own behalf, and not on +behalf of any Contributor. You must make it absolutely clear that any +such warranty, support, indemnity, or liability obligation is offered by +You alone, and You hereby agree to indemnify every Contributor for any +liability incurred by such Contributor as a result of warranty, support, +indemnity or liability terms You offer. You may include additional +disclaimers of warranty and limitations of liability specific to any +jurisdiction. + +4. Inability to Comply Due to Statute or Regulation +--------------------------------------------------- + +If it is impossible for You to comply with any of the terms of this +License with respect to some or all of the Covered Software due to +statute, judicial order, or regulation then You must: (a) comply with +the terms of this License to the maximum extent possible; and (b) +describe the limitations and the code they affect. Such description must +be placed in a text file included with all distributions of the Covered +Software under this License. Except to the extent prohibited by statute +or regulation, such description must be sufficiently detailed for a +recipient of ordinary skill to be able to understand it. + +5. Termination +-------------- + +5.1. The rights granted under this License will terminate automatically +if You fail to comply with any of its terms. However, if You become +compliant, then the rights granted under this License from a particular +Contributor are reinstated (a) provisionally, unless and until such +Contributor explicitly and finally terminates Your grants, and (b) on an +ongoing basis, if such Contributor fails to notify You of the +non-compliance by some reasonable means prior to 60 days after You have +come back into compliance. Moreover, Your grants from a particular +Contributor are reinstated on an ongoing basis if such Contributor +notifies You of the non-compliance by some reasonable means, this is the +first time You have received notice of non-compliance with this License +from such Contributor, and You become compliant prior to 30 days after +Your receipt of the notice. + +5.2. If You initiate litigation against any entity by asserting a patent +infringement claim (excluding declaratory judgment actions, +counter-claims, and cross-claims) alleging that a Contributor Version +directly or indirectly infringes any patent, then the rights granted to +You by any and all Contributors for the Covered Software under Section +2.1 of this License shall terminate. + +5.3. In the event of termination under Sections 5.1 or 5.2 above, all +end user license agreements (excluding distributors and resellers) which +have been validly granted by You or Your distributors under this License +prior to termination shall survive termination. + +************************************************************************ +* * +* 6. Disclaimer of Warranty * +* ------------------------- * +* * +* Covered Software is provided under this License on an "as is" * +* basis, without warranty of any kind, either expressed, implied, or * +* statutory, including, without limitation, warranties that the * +* Covered Software is free of defects, merchantable, fit for a * +* particular purpose or non-infringing. The entire risk as to the * +* quality and performance of the Covered Software is with You. * +* Should any Covered Software prove defective in any respect, You * +* (not any Contributor) assume the cost of any necessary servicing, * +* repair, or correction. This disclaimer of warranty constitutes an * +* essential part of this License. No use of any Covered Software is * +* authorized under this License except under this disclaimer. * +* * +************************************************************************ + +************************************************************************ +* * +* 7. Limitation of Liability * +* -------------------------- * +* * +* Under no circumstances and under no legal theory, whether tort * +* (including negligence), contract, or otherwise, shall any * +* Contributor, or anyone who distributes Covered Software as * +* permitted above, be liable to You for any direct, indirect, * +* special, incidental, or consequential damages of any character * +* including, without limitation, damages for lost profits, loss of * +* goodwill, work stoppage, computer failure or malfunction, or any * +* and all other commercial damages or losses, even if such party * +* shall have been informed of the possibility of such damages. This * +* limitation of liability shall not apply to liability for death or * +* personal injury resulting from such party's negligence to the * +* extent applicable law prohibits such limitation. Some * +* jurisdictions do not allow the exclusion or limitation of * +* incidental or consequential damages, so this exclusion and * +* limitation may not apply to You. * +* * +************************************************************************ + +8. Litigation +------------- + +Any litigation relating to this License may be brought only in the +courts of a jurisdiction where the defendant maintains its principal +place of business and such litigation shall be governed by laws of that +jurisdiction, without reference to its conflict-of-law provisions. +Nothing in this Section shall prevent a party's ability to bring +cross-claims or counter-claims. + +9. Miscellaneous +---------------- + +This License represents the complete agreement concerning the subject +matter hereof. If any provision of this License is held to be +unenforceable, such provision shall be reformed only to the extent +necessary to make it enforceable. Any law or regulation which provides +that the language of a contract shall be construed against the drafter +shall not be used to construe this License against a Contributor. + +10. Versions of the License +--------------------------- + +10.1. New Versions + +Mozilla Foundation is the license steward. Except as provided in Section +10.3, no one other than the license steward has the right to modify or +publish new versions of this License. Each version will be given a +distinguishing version number. + +10.2. Effect of New Versions + +You may distribute the Covered Software under the terms of the version +of the License under which You originally received the Covered Software, +or under the terms of any subsequent version published by the license +steward. + +10.3. Modified Versions + +If you create software not governed by this License, and you want to +create a new license for such software, you may create and use a +modified version of this License if you rename the license and remove +any references to the name of the license steward (except to note that +such modified license differs from this License). + +10.4. Distributing Source Code Form that is Incompatible With Secondary +Licenses + +If You choose to distribute Source Code Form that is Incompatible With +Secondary Licenses under the terms of this version of the License, the +notice described in Exhibit B of this License must be attached. + +Exhibit A - Source Code Form License Notice +------------------------------------------- + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + +If it is not possible or desirable to put the notice in a particular +file, then You may include the notice in a location (such as a LICENSE +file in a relevant directory) where a recipient would be likely to look +for such a notice. + +You may add additional accurate notices of copyright ownership. + +Exhibit B - "Incompatible With Secondary Licenses" Notice +--------------------------------------------------------- + + This Source Code Form is "Incompatible With Secondary Licenses", as + defined by the Mozilla Public License, v. 2.0. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,75 @@ +# hcdt + +### Delaunay triangulation + +```haskell +ghci> import Text.Show.Pretty +ghci> import Geometry.HCDT +ghci> vertices = [Vertex 0 0, Vertex 0 1, Vertex 1 1, Vertex 1 0] +ghci> triangulation <- delaunay vertices +ghci> pPrint triangulation +Triangulation + { _vertices = + fromList + [ ( 0 , Vertex 0.0 0.0 ) + , ( 1 , Vertex 0.0 1.0 ) + , ( 2 , Vertex 1.0 1.0 ) + , ( 3 , Vertex 1.0 0.0 ) + ] + , _triangles = [ Triangle 1 0 2 , Triangle 2 0 3 ] + , _edges = [ Edge 2 3 , Edge 0 3 , Edge 1 2 , Edge 0 2 , Edge 0 1 ] + } +ghci> pPrint $ borderEdges triangulation +fromList [ Edge 1 0 , Edge 1 2 , Edge 0 3 , Edge 2 3 ] +``` + + +### Constrained Delaunay triangulation + +```haskell +ghci> import Text.Show.Pretty +ghci> import Geometry.HCDT +ghci> vertices = [Vertex 1 1, Vertex 3 1, Vertex 2 2, Vertex 0 0, Vertex 4 0, Vertex 2 5] +ghci> edges = [Edge 0 1, Edge 0 2, Edge 1 2, Edge 3 4, Edge 3 5, Edge 4 5] +ghci> triangulation <- cdelaunay vertices edges +ntriangles: 6 +ghci> pPrint triangulation +ConstrainedTriangulation + { _triangulation = + Triangulation + { _vertices = + fromList + [ ( 0 , Vertex 1.0 1.0 ) + , ( 1 , Vertex 3.0 1.0 ) + , ( 2 , Vertex 2.0 2.0 ) + , ( 3 , Vertex 0.0 0.0 ) + , ( 4 , Vertex 4.0 0.0 ) + , ( 5 , Vertex 2.0 5.0 ) + ] + , _triangles = + [ Triangle 1 0 3 + , Triangle 2 1 5 + , Triangle 2 5 0 + , Triangle 0 5 3 + , Triangle 1 4 5 + , Triangle 1 3 4 + ] + , _edges = + [ Edge 3 4 + , Edge 1 4 + , Edge 3 5 + , Edge 0 5 + , Edge 2 5 + , Edge 0 2 + , Edge 1 5 + , Edge 1 2 + , Edge 1 3 + , Edge 4 5 + , Edge 0 3 + , Edge 0 1 + ] + } + , _fixedEdges = + [ Edge 0 1 , Edge 0 2 , Edge 4 5 , Edge 1 2 , Edge 3 4 , Edge 3 5 ] + } +``` diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cpp/hcdt.cpp b/cpp/hcdt.cpp new file mode 100644 --- /dev/null +++ b/cpp/hcdt.cpp @@ -0,0 +1,133 @@ +#include // to use malloc +#include // to use printf +#include "hcdt.hpp" +#include "CDT.h" + +extern "C" +{ + +TriangulationT* delaunay(VertexT* points, size_t npoints){ + CDT::Triangulation cdt(CDT::VertexInsertionOrder::AsProvided); + // insert vertices + std::vector> vertices(npoints); + for(size_t k = 0; k < npoints; ++k) { + const VertexT v = points[k]; + vertices[k] = CDT::V2d::make(v.x, v.y); + } + cdt.insertVertices(vertices); + cdt.eraseSuperTriangle(); + //// output + // triangles + const CDT::TriangleVec triangles = cdt.triangles; + const size_t ntriangles = triangles.size(); + TriangleT* out_triangles = (TriangleT*)malloc(ntriangles * sizeof(TriangleT)); + for(size_t k = 0; k < ntriangles; ++k){ + const CDT::VerticesArr3 trgl = triangles[k].vertices; + out_triangles[k].i1 = trgl[0]; + out_triangles[k].i2 = trgl[1]; + out_triangles[k].i3 = trgl[2]; + } + // all edges + CDT::EdgeUSet allEdges = CDT::extractEdgesFromTriangles(triangles); + const size_t nedges = allEdges.size(); + EdgeT* out_alledges = (EdgeT*)malloc(nedges * sizeof(EdgeT)); + std::unordered_set :: iterator it; + size_t k = 0; + for(it = allEdges.begin(); it != allEdges.end(); it++){ + const CDT::Edge edge = *it; + out_alledges[k].i = CDT::edge_get_v1(edge); + out_alledges[k].j = CDT::edge_get_v2(edge); + k++; + } + /* output mesh */ + TriangulationT* out = (TriangulationT*)malloc(sizeof(TriangulationT)); + out->vertices = points; + out->nvertices = npoints; + out->triangles = out_triangles; + out->ntriangles = ntriangles; + out->edges = out_alledges; + out->nedges = nedges; + return out; +} + +CTriangulationT* cdelaunay(VertexT* points, size_t npoints, EdgeT* fedges, size_t nfedges){ + CDT::Triangulation cdt( + CDT::VertexInsertionOrder::AsProvided, + CDT::IntersectingConstraintEdges::Resolve, + 0.0 + ); + // insert vertices + std::vector> vertices(npoints); + for(size_t k = 0; k < npoints; ++k) { + const VertexT v = points[k]; + vertices[k] = CDT::V2d::make(v.x, v.y); + } + // insert edges + std::vector Edges; + Edges.reserve(nfedges); + for (size_t k = 0; k < nfedges; ++k) { + const EdgeT e = fedges[k]; + Edges.push_back(CDT::Edge(e.i, e.j)); + } + cdt.insertVertices(vertices); + cdt.insertEdges(Edges); + cdt.eraseOuterTrianglesAndHoles(); + //// output + // vertices + const std::vector> cdt_vertices = cdt.vertices; + const size_t nvertices = cdt_vertices.size(); + VertexT* out_vertices = (VertexT*)malloc(nvertices * sizeof(VertexT)); + for(size_t k = 0; k < nvertices; ++k){ + const CDT::V2d v = cdt_vertices[k]; + out_vertices[k].x = v.x; + out_vertices[k].y = v.y; + } + // triangles + const CDT::TriangleVec triangles = cdt.triangles; + const size_t ntriangles = triangles.size(); + //printf("ntriangles: %u\n", ntriangles); + TriangleT* out_triangles = (TriangleT*)malloc(ntriangles * sizeof(TriangleT)); + for(size_t k = 0; k < ntriangles; ++k){ + const CDT::VerticesArr3 trgl = triangles[k].vertices; + out_triangles[k].i1 = trgl[0]; + out_triangles[k].i2 = trgl[1]; + out_triangles[k].i3 = trgl[2]; + } + // all edges + CDT::EdgeUSet allEdges = CDT::extractEdgesFromTriangles(triangles); + const size_t nedges = allEdges.size(); + EdgeT* out_alledges = (EdgeT*)malloc(nedges * sizeof(EdgeT)); + std::unordered_set :: iterator it; + size_t k = 0; + for(it = allEdges.begin(); it != allEdges.end(); it++){ + const CDT::Edge edge = *it; + out_alledges[k].i = CDT::edge_get_v1(edge); + out_alledges[k].j = CDT::edge_get_v2(edge); + k++; + } + // fixed edges + CDT::EdgeUSet fixedEdges = cdt.fixedEdges; + const size_t nfixededges = fixedEdges.size(); + EdgeT* out_fixededges = (EdgeT*)malloc(nfixededges * sizeof(EdgeT)); + std::unordered_set :: iterator itedge; + size_t l = 0; + for(itedge = fixedEdges.begin(); itedge != fixedEdges.end(); itedge++){ + const CDT::Edge edge = *itedge; + out_fixededges[l].i = CDT::edge_get_v1(edge); + out_fixededges[l].j = CDT::edge_get_v2(edge); + l++; + } + /* output mesh */ + CTriangulationT* out = (CTriangulationT*)malloc(sizeof(CTriangulationT)); + out->vertices = out_vertices; + out->nvertices = nvertices; + out->triangles = out_triangles; + out->ntriangles = ntriangles; + out->edges = out_alledges; + out->nedges = nedges; + out->fixededges = out_fixededges; + out->nfixededges = nfixededges; + return out; +} + +} diff --git a/hcdt.cabal b/hcdt.cabal new file mode 100644 --- /dev/null +++ b/hcdt.cabal @@ -0,0 +1,44 @@ +cabal-version: 2.2 +name: hcdt +version: 0.1.0.0 +synopsis: 2d Delaunay triangulation +description: + This library performs the constrained or unconstrained 2d Delaunay triangulation. + . + It uses the C++ library 'CDT'. + . + For examples, look the README file. +homepage: https://github.com/githubuser/hcdt#readme +license: BSD-3-Clause +license-file: LICENSE +author: Stéphane Laurent +maintainer: laurent_step@outlook.fr +copyright: 2022 Stéphane Laurent +category: Geometry +build-type: Simple +extra-source-files: README.md + CHANGELOG.md + LICENSE_CDT.txt + +library + hs-source-dirs: src + exposed-modules: Geometry.HCDT + other-modules: Geometry.HCDT.Types + , Geometry.HCDT.Triangulation + , Geometry.HCDT.Delaunay + , Geometry.HCDT.Utils + build-depends: base >= 4.7 && < 5 + , containers >= 0.6.4.1 && < 0.7 + , indexed-traversable >= 0.1.2 && < 0.2 + include-dirs: cpp + C-sources: cpp/hcdt.cpp + extra-libraries: stdc++ + ghc-options: -Wall -optcxx-std=c++11 + other-extensions: CPP + , ForeignFunctionInterface + cxx-options: -fPIC -std=c++11 + default-language: Haskell2010 + +source-repository head + type: git + location: https://github.com/githubuser/hcdt diff --git a/src/Geometry/HCDT.hs b/src/Geometry/HCDT.hs new file mode 100644 --- /dev/null +++ b/src/Geometry/HCDT.hs @@ -0,0 +1,6 @@ +module Geometry.HCDT + (module X) + where +import Geometry.HCDT.Delaunay as X +import Geometry.HCDT.Types as X +import Geometry.HCDT.Utils as X diff --git a/src/Geometry/HCDT/Delaunay.hs b/src/Geometry/HCDT/Delaunay.hs new file mode 100644 --- /dev/null +++ b/src/Geometry/HCDT/Delaunay.hs @@ -0,0 +1,42 @@ +module Geometry.HCDT.Delaunay + ( delaunay, + cdelaunay + ) +where + +import Foreign.Marshal.Alloc (free, mallocBytes) +import Foreign.Marshal.Array (pokeArray) +import Foreign.Storable (peek, sizeOf) +import Geometry.HCDT.Triangulation +import Geometry.HCDT.Types (Triangulation, Vertex, ConstrainedTriangulation, Edge) + +-- | Unconstrained 2d Delaunay triangulation. +delaunay :: [Vertex] -> IO Triangulation +delaunay vertices = do + let nvertices = length vertices + verticesPtr <- mallocBytes (nvertices * sizeOf (undefined :: CVertex)) + cvertices <- mapM vertexToCVertex vertices + pokeArray verticesPtr cvertices + ctriangulationPtr <- c_delaunay verticesPtr (fromIntegral nvertices) + ctriangulation <- peek ctriangulationPtr + free verticesPtr + free ctriangulationPtr + cTriangulationToTriangulation ctriangulation + +-- | Constrained 2d Delaunay triangulation. +cdelaunay :: [Vertex] -> [Edge] -> IO ConstrainedTriangulation +cdelaunay vertices edges = do + let nvertices = length vertices + verticesPtr <- mallocBytes (nvertices * sizeOf (undefined :: CVertex)) + cvertices <- mapM vertexToCVertex vertices + pokeArray verticesPtr cvertices + let nedges = length edges + edgesPtr <- mallocBytes (nedges * sizeOf (undefined :: CEdge)) + cedges <- mapM edgeToCEdge edges + pokeArray edgesPtr cedges + cctriangulationPtr <- c_cdelaunay verticesPtr (fromIntegral nvertices) edgesPtr (fromIntegral nedges) + cctriangulation <- peek cctriangulationPtr + free verticesPtr + free edgesPtr + free cctriangulationPtr + cCTriangulationToConstrainedTriangulation cctriangulation diff --git a/src/Geometry/HCDT/Triangulation.hs b/src/Geometry/HCDT/Triangulation.hs new file mode 100644 --- /dev/null +++ b/src/Geometry/HCDT/Triangulation.hs @@ -0,0 +1,277 @@ +{-# LINE 1 "Triangulation.hsc" #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE CPP #-} +module Geometry.HCDT.Triangulation +( cTriangulationToTriangulation +, cCTriangulationToConstrainedTriangulation +, vertexToCVertex +, edgeToCEdge +, c_delaunay +, c_cdelaunay +, CTriangulation (..) +, CVertex (..) +, CCTriangulation (..) +, CEdge (..) ) + where +import Geometry.HCDT.Types +import Foreign +import Foreign.C.Types +import qualified Data.IntMap.Strict as IM + + + +data CVertex = CVertex { + __x :: CDouble + , __y :: CDouble +} + +instance Storable CVertex where + sizeOf __ = (16) +{-# LINE 29 "Triangulation.hsc" #-} + alignment __ = 8 +{-# LINE 30 "Triangulation.hsc" #-} + peek ptr = do + x' <- (`peekByteOff` 0) ptr +{-# LINE 32 "Triangulation.hsc" #-} + y' <- (`peekByteOff` 8) ptr +{-# LINE 33 "Triangulation.hsc" #-} + return CVertex { __x = x', __y = y' } + poke ptr (CVertex r1 r2) + = do + (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1 +{-# LINE 37 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2 +{-# LINE 38 "Triangulation.hsc" #-} + +cVertexToVertex :: CVertex -> IO Vertex +cVertexToVertex cvertex = do + let x = realToFrac $ __x cvertex + let y = realToFrac $ __y cvertex + return $ Vertex x y + +vertexToCVertex :: Vertex -> IO CVertex +vertexToCVertex (Vertex x y) = do + return $ CVertex { __x = realToFrac x, __y = realToFrac y } + +data CEdge = CEdge { + __i :: CUInt + , __j :: CUInt +} + +instance Storable CEdge where + sizeOf __ = (8) +{-# LINE 56 "Triangulation.hsc" #-} + alignment __ = 4 +{-# LINE 57 "Triangulation.hsc" #-} + peek ptr = do + i' <- (`peekByteOff` 0) ptr +{-# LINE 59 "Triangulation.hsc" #-} + j' <- (`peekByteOff` 4) ptr +{-# LINE 60 "Triangulation.hsc" #-} + return CEdge { __i = i' + , __j = j' } + poke ptr (CEdge r1 r2) + = do + (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1 +{-# LINE 65 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 4) ptr r2 +{-# LINE 66 "Triangulation.hsc" #-} + +cEdgeToEdge :: CEdge -> IO Edge +cEdgeToEdge cedge = do + let i = fromIntegral $ __i cedge + let j = fromIntegral $ __j cedge + return $ Edge i j + +edgeToCEdge :: Edge -> IO CEdge +edgeToCEdge (Edge i j) = do + return $ CEdge { __i = fromIntegral i, __j = fromIntegral j } + +data CTriangle = CTriangle { + __i1 :: CUInt + , __i2 :: CUInt + , __i3 :: CUInt +} + +instance Storable CTriangle where + sizeOf __ = (12) +{-# LINE 85 "Triangulation.hsc" #-} + alignment __ = 4 +{-# LINE 86 "Triangulation.hsc" #-} + peek ptr = do + i1' <- (`peekByteOff` 0) ptr +{-# LINE 88 "Triangulation.hsc" #-} + i2' <- (`peekByteOff` 4) ptr +{-# LINE 89 "Triangulation.hsc" #-} + i3' <- (`peekByteOff` 8) ptr +{-# LINE 90 "Triangulation.hsc" #-} + return CTriangle { __i1 = i1' + , __i2 = i2' + , __i3 = i3' } + poke ptr (CTriangle r1 r2 r3) + = do + (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1 +{-# LINE 96 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 4) ptr r2 +{-# LINE 97 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r3 +{-# LINE 98 "Triangulation.hsc" #-} + +cTriangleToTriangle :: CTriangle -> IO Triangle +cTriangleToTriangle ctriangle = do + let i1 = fromIntegral $ __i1 ctriangle + let i2 = fromIntegral $ __i2 ctriangle + let i3 = fromIntegral $ __i3 ctriangle + return $ Triangle i1 i2 i3 + +data CTriangulation = CTriangulation { + __vertices :: Ptr CVertex + , __nvertices :: CSize + , __triangles :: Ptr CTriangle + , __ntriangles :: CSize + , __edges :: Ptr CEdge + , __nedges :: CSize +} + +instance Storable CTriangulation where + sizeOf __ = (48) +{-# LINE 117 "Triangulation.hsc" #-} + alignment __ = 8 +{-# LINE 118 "Triangulation.hsc" #-} + peek ptr = do + vs <- (`peekByteOff` 0) ptr +{-# LINE 120 "Triangulation.hsc" #-} + nvs <- (`peekByteOff` 8) ptr +{-# LINE 121 "Triangulation.hsc" #-} + ts <- (`peekByteOff` 16) ptr +{-# LINE 122 "Triangulation.hsc" #-} + nts <- (`peekByteOff` 24) ptr +{-# LINE 123 "Triangulation.hsc" #-} + es <- (`peekByteOff` 32) ptr +{-# LINE 124 "Triangulation.hsc" #-} + nes <- (`peekByteOff` 40) ptr +{-# LINE 125 "Triangulation.hsc" #-} + return CTriangulation { __vertices = vs + , __nvertices = nvs + , __triangles = ts + , __ntriangles = nts + , __edges = es + , __nedges = nes } + poke ptr (CTriangulation r1 r2 r3 r4 r5 r6) + = do + (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1 +{-# LINE 134 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2 +{-# LINE 135 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 16) ptr r3 +{-# LINE 136 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 24) ptr r4 +{-# LINE 137 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 32) ptr r5 +{-# LINE 138 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 40) ptr r6 +{-# LINE 139 "Triangulation.hsc" #-} + +data CCTriangulation = CCTriangulation { + __vertices' :: Ptr CVertex + , __nvertices' :: CSize + , __triangles' :: Ptr CTriangle + , __ntriangles' :: CSize + , __edges' :: Ptr CEdge + , __nedges' :: CSize + , __fixededges' :: Ptr CEdge + , __nfixededges' :: CSize +} + +instance Storable CCTriangulation where + sizeOf __ = (64) +{-# LINE 153 "Triangulation.hsc" #-} + alignment __ = 8 +{-# LINE 154 "Triangulation.hsc" #-} + peek ptr = do + vs <- (`peekByteOff` 0) ptr +{-# LINE 156 "Triangulation.hsc" #-} + nvs <- (`peekByteOff` 8) ptr +{-# LINE 157 "Triangulation.hsc" #-} + ts <- (`peekByteOff` 16) ptr +{-# LINE 158 "Triangulation.hsc" #-} + nts <- (`peekByteOff` 24) ptr +{-# LINE 159 "Triangulation.hsc" #-} + es <- (`peekByteOff` 32) ptr +{-# LINE 160 "Triangulation.hsc" #-} + nes <- (`peekByteOff` 40) ptr +{-# LINE 161 "Triangulation.hsc" #-} + fes <- (`peekByteOff` 48) ptr +{-# LINE 162 "Triangulation.hsc" #-} + nfes <- (`peekByteOff` 56) ptr +{-# LINE 163 "Triangulation.hsc" #-} + return CCTriangulation { __vertices' = vs + , __nvertices' = nvs + , __triangles' = ts + , __ntriangles' = nts + , __edges' = es + , __nedges' = nes + , __fixededges' = fes + , __nfixededges' = nfes } + poke ptr (CCTriangulation r1 r2 r3 r4 r5 r6 r7 r8) + = do + (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1 +{-# LINE 174 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2 +{-# LINE 175 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 16) ptr r3 +{-# LINE 176 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 24) ptr r4 +{-# LINE 177 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 32) ptr r5 +{-# LINE 178 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 40) ptr r6 +{-# LINE 179 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 48) ptr r7 +{-# LINE 180 "Triangulation.hsc" #-} + (\hsc_ptr -> pokeByteOff hsc_ptr 56) ptr r8 +{-# LINE 181 "Triangulation.hsc" #-} + +cTriangulationToTriangulation :: CTriangulation -> IO Triangulation +cTriangulationToTriangulation ctriangulation = do + let nvertices = fromIntegral $ __nvertices ctriangulation + ntriangles = fromIntegral $ __ntriangles ctriangulation + nedges = fromIntegral $ __nedges ctriangulation + vertices <- peekArray nvertices (__vertices ctriangulation) + triangles <- peekArray ntriangles (__triangles ctriangulation) + edges <- peekArray nedges (__edges ctriangulation) + vertices' <- mapM cVertexToVertex vertices + triangles' <- mapM cTriangleToTriangle triangles + edges' <- mapM cEdgeToEdge edges + return $ Triangulation { _vertices = IM.fromAscList (zip [0 .. nvertices-1] vertices') + , _triangles = triangles' + , _edges = edges' } + +cCTriangulationToConstrainedTriangulation :: CCTriangulation -> IO ConstrainedTriangulation +cCTriangulationToConstrainedTriangulation cctriangulation = do + let nvertices = fromIntegral $ __nvertices' cctriangulation + ntriangles = fromIntegral $ __ntriangles' cctriangulation + nedges = fromIntegral $ __nedges' cctriangulation + nfedges = fromIntegral $ __nfixededges' cctriangulation + vertices <- peekArray nvertices (__vertices' cctriangulation) + triangles <- peekArray ntriangles (__triangles' cctriangulation) + edges <- peekArray nedges (__edges' cctriangulation) + fedges <- peekArray nfedges (__fixededges' cctriangulation) + vertices' <- mapM cVertexToVertex vertices + triangles' <- mapM cTriangleToTriangle triangles + edges' <- mapM cEdgeToEdge edges + fedges' <- mapM cEdgeToEdge fedges + let triangulation = Triangulation { + _vertices = IM.fromAscList (zip [0 .. nvertices-1] vertices') + , _triangles = triangles' + , _edges = edges' + } + return ConstrainedTriangulation { + _triangulation = triangulation + , _fixedEdges = fedges' } + +foreign import ccall unsafe "delaunay" c_delaunay + :: Ptr CVertex -> CSize -> IO (Ptr CTriangulation) + +foreign import ccall unsafe "cdelaunay" c_cdelaunay + :: Ptr CVertex -> CSize -> Ptr CEdge -> CSize -> IO (Ptr CCTriangulation) diff --git a/src/Geometry/HCDT/Types.hs b/src/Geometry/HCDT/Types.hs new file mode 100644 --- /dev/null +++ b/src/Geometry/HCDT/Types.hs @@ -0,0 +1,28 @@ +module Geometry.HCDT.Types + where +import Data.IntMap.Strict (IntMap) +import Data.List (sort) + +data Triangle = Triangle Int Int Int + deriving (Show, Read) +instance Eq Triangle where + Triangle i j k == Triangle i' j' k' = sort [i, j, k] == sort [i', j', k'] + +data Edge = Edge Int Int + deriving (Show, Read) +instance Eq Edge where + Edge i j == Edge i' j' = (i == i' && j == j') || (i == j' && j == i') + +data Vertex = Vertex Double Double + deriving (Show, Eq) + +data Triangulation = Triangulation { + _vertices :: IntMap Vertex + , _triangles :: [Triangle] + , _edges :: [Edge] +} deriving Show + +data ConstrainedTriangulation = ConstrainedTriangulation { + _triangulation :: Triangulation + , _fixedEdges :: [Edge] +} deriving Show diff --git a/src/Geometry/HCDT/Utils.hs b/src/Geometry/HCDT/Utils.hs new file mode 100644 --- /dev/null +++ b/src/Geometry/HCDT/Utils.hs @@ -0,0 +1,47 @@ +module Geometry.HCDT.Utils + ( borderEdges, + partitionEdges + ) +where + +import Geometry.HCDT.Types +import Data.Sequence as DS (Seq, fromList, (><), filter, drop) +import Data.Maybe (isNothing, fromJust) +import Data.Foldable (find) +import Data.Foldable.WithIndex (ifind) + +triangleEdges :: Triangle -> Seq Edge +triangleEdges (Triangle i j k) = fromList [Edge i j, Edge j k, Edge i k] + +-- isUnique :: Eq a => Seq a -> a -> Bool +-- isUnique xs x = length (DS.filter (== x) xs) == 1 + +allEdges :: Triangulation -> Seq Edge +allEdges triangulation = + let triplets = map triangleEdges (_triangles triangulation) in + foldl1 (><) triplets + +-- borderEdges :: Triangulation -> Seq Edge +-- borderEdges triangulation = DS.filter (isUnique edges) edges +-- where +-- edges = allEdges triangulation + +isUnique :: Eq a => Seq a -> a -> Bool +isUnique xs x = isNothing y + where + (i, _) = fromJust $ ifind (\_ x' -> x' == x) xs + y = find (== x) (DS.drop i xs) + +-- | Exterior edges of a Delaunay triangulation. +borderEdges :: Triangulation -> Seq Edge +borderEdges triangulation = DS.filter (isUnique edges) edges + where + edges = allEdges triangulation + +-- | Exterior and interior edges of a Delaunay triangulation. +partitionEdges :: Triangulation -> (Seq Edge, Seq Edge) +partitionEdges triangulation = (exEdges, inEdges) + where + edges = allEdges triangulation + exEdges = DS.filter (isUnique edges) edges + inEdges = DS.filter (not . isUnique edges) edges