diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d5fabb4 --- /dev/null +++ b/.gitignore @@ -0,0 +1,104 @@ +# PyCharm IDE +.idea/* + +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +.hypothesis/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# pyenv +.python-version + +# celery beat schedule file +celerybeat-schedule + +# SageMath parsed files +*.sage.py + +# dotenv +.env + +# virtualenv +.venv +venv/ +ENV/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..4d1b072 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "external/bih"] + path = external/bih + url = git@github.com:flow123d/bih.git +[submodule "external/matlab-intersections"] + path = external/matlab-intersections + url = git@github.com:jirikopal/Intersections.git diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..0c94ad8 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,8 @@ +nguage: python +python: + - "3.5" + - "3.6" +# command to install dependencies +install: "pip3 install -r requirements.txt" +# command to run tests +script: pytest diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..9cecc1d --- /dev/null +++ b/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + {one line to give the program's name and a brief idea of what it does.} + Copyright (C) {year} {name of author} + + 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, either version 3 of the License, or + (at your option) any later version. + + 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 . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + {project} Copyright (C) {year} {fullname} + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/README.md b/README.md new file mode 100644 index 0000000..3bbb984 --- /dev/null +++ b/README.md @@ -0,0 +1,27 @@ +Intersections +============== + +[![Build Status](https://travis-ci.org/GeoMop/Intersections.svg?branch=master)](https://travis-ci.org/GeoMop/Intersections) +[![Code Health](https://landscape.io/github/GeoMop/Intersections/master/landscape.svg?style=flat)](https://landscape.io/github/GeoMop/Intersections/master) +[![Code Climate](https://codeclimate.com/github/GeoMop/Intersections/badges/gpa.svg)](https://codeclimate.com/github/GeoMop/Intersections) +[![Test Coverage](https://codeclimate.com/github/GeoMop/Intersections/badges/coverage.svg)](https://codeclimate.com/github/GeoMop/Intersections/coverage) + + +Computing intersections of B-spline curves and surfaces. + +Library focus on fast intersection algorithms for non-degenerate intersections of B-spline curves and surfaces +of small degree (especially quadratic). + +Requirements +------------ + +* g++ 4.x or newer +* cmake 3.x + +In order to install BIH package locally for development just run the 'bih_install.sh' script. + +Theory +------ +[Patrikalakis-Maekawa-Cho](http://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/mathe.html) + + \ No newline at end of file diff --git a/basetestvec.m b/basetestvec.m deleted file mode 100644 index 64ec858..0000000 --- a/basetestvec.m +++ /dev/null @@ -1,21 +0,0 @@ -function [ ] = basetestvec( knots ) -n = 1000; -n_basf = length(knots)-3; -one = ones(n_basf,1); -y = sparse(zeros(n,n_basf)); -for i=1:n - y(i,:) = splinebasevec(knots,min(knots) + max(knots)*(i-1)/(n-1),0); -end -% spy(y) -% pause -x = linspace(min(knots),max(knots),n); -x -y -%spy(y) -plot(x,y) -%figure -if norm(y*one - ones(n,1))< n_basf*eps - disp('OK') -end -end - diff --git a/bih_install.sh b/bih_install.sh new file mode 100755 index 0000000..d293b8a --- /dev/null +++ b/bih_install.sh @@ -0,0 +1,17 @@ +#!/bin/bash + +echo "Installing BIH locally." +echo "g++ 4.x and cmake 3.x or newer are assumed." + +git submodule update --init --recursive +cd external/bih +mkdir build +cd build +cmake .. +make + +BIH_PY_PATH=`pwd` +echo "Modifiing your .bashrc:" +echo "Add "$BIH_PY_PATH" to PYTHONPATH" +echo "PYTHONPATH=\"\$PYTHONPATH:$BIH_PY_PATH\"" >>"${HOME}/.bashrc" + diff --git a/boundary_point.m b/boundary_point.m deleted file mode 100644 index 353766d..0000000 --- a/boundary_point.m +++ /dev/null @@ -1,55 +0,0 @@ -function [ bool ] = boundary_point(uv,iujv,u_knots,v_knots) -bool = [-1,-1]; - -uint = [u_knots(iujv(1)+2), u_knots(iujv(1)+3)]; -vint = [v_knots(iujv(2)+2), v_knots(iujv(2)+3)]; -uv; - - -if uint(1) == 0 - if abs(uv(1)-0)= mnXs) && (mnX <= mxXs)) || ... - ((mxX >= mnXs) && (mxX <= mxXs)) ) == 1 - - if (( (mnY >= mnYs) && (mnY <= mxYs)) || ... - ((mxY >= mnYs) && (mxY <= mxYs)) ) == 1 - if (( (mnZ >= mnZs) && (mnZ <= mxZs)) || ... - ((mxZ >= mnZs) && (mxZ <= mxZs)) ) == 1 - - - n_its(k,l) = n_its(k,l) + 1; - its(k,l,n_its(k,l)) = sp_i; - - end - end - end - end - end - end -end - - -end - diff --git a/build_LS_matrix.m b/build_LS_matrix.m deleted file mode 100644 index fcd2360..0000000 --- a/build_LS_matrix.m +++ /dev/null @@ -1,14 +0,0 @@ -function [ B,Interv ] = build_LS_matrix( u_knots,v_knots, Xp ) -u_n_basf = length(u_knots)-3; -v_n_basf = length(v_knots)-3; -[np k] = size(Xp); % k unused -B =spalloc(np,u_n_basf*v_n_basf,9*np); -Interv = zeros(np,2); -for j = 1:np - [uf, k] = splinebasevec(u_knots,Xp(j,1),0); - [vf, i] = splinebasevec(v_knots,Xp(j,2),0); - Interv(j,:) = [k,i]; - B(j,:) = kron(vf',uf'); -end -end - diff --git a/build_reg_matrix.m b/build_reg_matrix.m deleted file mode 100644 index 6f79078..0000000 --- a/build_reg_matrix.m +++ /dev/null @@ -1,82 +0,0 @@ -function [ A ] = build_reg_matrix( u_knots,v_knots, P0,P1,P2,P3,nnzA ) - -u_n_basf = length(u_knots)-3; -v_n_basf = length(v_knots)-3; -u_n_inter = length(u_knots) - 5 -v_n_inter = length(v_knots) - 5 - - a = P3 - P2; - b = P0 - P1; - c = P1 - P2; - d = P0 - P3; - -A =spalloc(u_n_basf*v_n_basf,u_n_basf*v_n_basf,nnzA); - -% Qpoints = [0 (1/2 - 1/sqrt(20)) (1/2 + 1/sqrt(20)) 1]; -% weights = [1/6 5/6 5/6 1/6]; - -Qpoints = [-0.90618 -0.538469 0 0.538469 0.90618] /2 + 0.5; -weights = [0.236927 0.478629 0.568889 0.478629 0.236927]; - - - -n_points = length(Qpoints); - -u_point_val = spalloc(u_n_basf,u_n_inter*n_points,u_n_inter*n_points*3); -ud_point_val = spalloc(u_n_basf,u_n_inter*n_points,u_n_inter*n_points*3); -q_u_point = zeros(u_n_inter*n_points,1); - -n = 0; -for i = 1:u_n_inter - us = u_knots(i+2); - uil = u_knots(i+3)- u_knots(i+2); - for k = 1:n_points - up = us + uil*Qpoints(k); - n = n+1; - q_u_point(n) = up; - u_point_val(:,n) = splinebasevec(u_knots,up,0,i); - ud_point_val(:,n) = splinebasevec(u_knots,up,1,i); - end -end - -%%% - -v_point_val = spalloc(v_n_basf,v_n_inter*n_points,v_n_inter*n_points*3); -vd_point_val = spalloc(v_n_basf,v_n_inter*n_points,v_n_inter*n_points*3); -q_v_point = zeros(v_n_inter*n_points,1); - -n = 0; -for i = 1:v_n_inter - vs = v_knots(i+2); - vil = v_knots(i+3)- v_knots(i+2); - for k = 1:n_points - vp = vs + vil*Qpoints(k); - n = n+1; - q_v_point(n) = vp; - v_point_val(:,n) = splinebasevec(v_knots,vp,0,i); - vd_point_val(:,n) = splinebasevec(v_knots,vp,1,i); - end -end - -%%% - -for i= 1:v_n_inter - for k = 1:n_points - v_point = v_point_val(:,(i-1)*n_points+k); - vd_point =vd_point_val(:,(i-1)*n_points+k); - for l =1:u_n_inter - for m = 1:n_points - u_point = u_point_val(:,(l-1)*n_points+m); - ud_point = ud_point_val(:,(l-1)*n_points+m); - vd = kron(vd_point,u_point); - ud = kron(v_point,ud_point); - v = q_v_point((i-1)*n_points +k); - u = q_u_point((l-1)*n_points +m); - J = det([v * a + (1 - v) * b, u * c + (1 - u) * d]); - A = A + J * weights(m)*weights(k)*(kron(ud,ud') +kron(vd,vd')); - end - end - end -end - -end diff --git a/compute_control_points.m b/compute_control_points.m deleted file mode 100644 index e93861b..0000000 --- a/compute_control_points.m +++ /dev/null @@ -1,23 +0,0 @@ -function [ X_coor,Y_coor ] = compute_control_points( u_knots, v_knots, P0,P1,P2,P3) -%UNTITLED Summary of this function goes here -% Detailed explanation goes here -u_n_basf = length(u_knots)-3; -v_n_basf = length(v_knots)-3; - - -X_coor = zeros(u_n_basf,v_n_basf); -Y_coor = zeros(u_n_basf,v_n_basf); - -for i =1:u_n_basf - u = (u_knots(i+1)+u_knots(i+2))/2; - for j =1:v_n_basf - v = (v_knots(j+1)+v_knots(j+2))/2; - P = u * (v * P2 + (1-v) * P1) + (1-u) * ( v * P3 + (1-v) * P0) ; - X_coor(i,j) = P(1); - Y_coor(i,j) = P(2); - end -end - - -end - diff --git a/compute_grid_lines.m b/compute_grid_lines.m deleted file mode 100644 index 41deea4..0000000 --- a/compute_grid_lines.m +++ /dev/null @@ -1,21 +0,0 @@ -function [ GX,GY ] = compute_grid_lines(u_knots,v_knots, P0,P1,P2,P3 ) - -u_n_grid = length(u_knots)-2; -v_n_grid = length(v_knots)-2; - -GX = zeros(u_n_grid,v_n_grid); -GY = zeros(u_n_grid,v_n_grid); - - -for i =1:u_n_grid - u = u_knots(i+2); - for j =1:v_n_grid - v = v_knots(j+2); - P = u * (v * P2 + (1-v) * P1) + (1-u) * ( v * P3 + (1-v) * P0) ; - GX(i,j) = P(1); - GY(i,j) = P(2); - end -end - -end - diff --git a/compute_patch_edges.m b/compute_patch_edges.m deleted file mode 100644 index ff068f3..0000000 --- a/compute_patch_edges.m +++ /dev/null @@ -1,23 +0,0 @@ -function [ X_coor,Y_coor ] = compute_patch_edges( u_knots, v_knots, P0,P1,P2,P3) -%UNTITLED Summary of this function goes here -% Detailed explanation goes here -u_n_bound = length(u_knots)-4; -v_n_bound = length(v_knots)-4; - - -X_coor = zeros(u_n_bound,v_n_bound); -Y_coor = zeros(u_n_bound,v_n_bound); - -for i =1:u_n_bound - u = u_knots(i+2); - for j =1:v_n_bound - v = v_knots(j+2); - P = u * (v * P2 + (1-v) * P1) + (1-u) * ( v * P3 + (1-v) * P0) ; - X_coor(i,j) = P(1); - Y_coor(i,j) = P(2); - end -end - - -end - diff --git a/external/bih b/external/bih new file mode 160000 index 0000000..2d7b774 --- /dev/null +++ b/external/bih @@ -0,0 +1 @@ +Subproject commit 2d7b774bae113c23e07b39b59253f9c49cce075e diff --git a/external/matlab-intersections b/external/matlab-intersections new file mode 160000 index 0000000..2bf73a4 --- /dev/null +++ b/external/matlab-intersections @@ -0,0 +1 @@ +Subproject commit 2bf73a4a4f65aefdc17c955c84471897b67bc32a diff --git a/find_int.m b/find_int.m deleted file mode 100644 index bfdd725..0000000 --- a/find_int.m +++ /dev/null @@ -1,59 +0,0 @@ -function [ k ] = find_int( T,t ) - -n = length(T); -mn = 3; -mx = n-2; -est = 3+floor(t*(n -5)); - -if t == T(mn) - k = mn-2; -elseif t == T(mx) - k = mx-3; -end - -if t>= T(est) - mn =est; -elseif t < T(est) % = - mx = est; -end - -s = max(ceil(log2(mx-mn)),1); - -for p=1:s+1 - if t < T(mn+1) % = - k = mn-2; - break - - elseif t>T(mx-1) % = - k = mx-3; - break - end - mid = mn + floor((mx-mn)/2); - %mn - %mx - if mid ~= mn - if t< T(mid) - mx = mid; - elseif t>= T(mid) - mn = mid; - end - else - k = mn-2; - break - end -end - -% T -% [T(mn) T(mx)] -% if (t~=0) && (t~=1) -% if (t>=T(k) && t<=T(k+1) ) -% [T T >t] -% t -% k -% disp('problem') -% pause -% end -% end - -end - diff --git a/get_errors.m b/get_errors.m deleted file mode 100644 index 266a59a..0000000 --- a/get_errors.m +++ /dev/null @@ -1,14 +0,0 @@ -function [ Err ] = get_errors(err,Interv,u_n,v_n ) - -n = length(err); -Err = zeros(v_n-2,u_n-2); -for j =1:n - Err(Interv(j,2),Interv(j,1)) = Err(Interv(j,2),Interv(j,1)) + err(j); -end - -%UNTITLED2 Summary of this function goes here -% Detailed explanation goes here - - -end - diff --git a/get_knot_vector.m b/get_knot_vector.m deleted file mode 100644 index a7bc984..0000000 --- a/get_knot_vector.m +++ /dev/null @@ -1,10 +0,0 @@ -function [ knots ] = get_knot_vector( n_basf ) - - -knots = zeros(n_basf+3,1); -knots(3:n_basf+1) = linspace(0,1,n_basf-1); -knots(n_basf+2:n_basf+3) = 1; - - -end - diff --git a/getpoints.m b/getpoints.m deleted file mode 100644 index 80211c3..0000000 --- a/getpoints.m +++ /dev/null @@ -1,23 +0,0 @@ -function [ X ] = getpoints() - -nu = 36; -nv = 26; -h = 0; -X = zeros(nu*nv,4); - -for i=1:nu - for j = 1:nv - h = h+1; - x = 2*(i-1)/(nu-1); - y = 2*(j-1)/(nv-1); - X(h,1) = x; - X(h,2) = y; - %X(h,3) = 0.3*cos(3*pi*exp(-x))*cos(3*pi*y^1.2)*cos(3*pi*sqrt(y^2+x^2)); - %X(h,3) = x^2+y^2 + 1; - X(h,3) = (2-x^2)+(2-y^2) + 1; - X(h,4) = 1; - end -end - -end - diff --git a/getpoints2.m b/getpoints2.m deleted file mode 100644 index 3f23f9a..0000000 --- a/getpoints2.m +++ /dev/null @@ -1,23 +0,0 @@ -function [ X ] = getpoints2() - -nu = 20; -nv = 20; -h = 0; -X = zeros(nu*nv,4); - -for i=1:nu - for j = 1:nv - h = h+1; - x = 2*(i-1)/(nu-1); - y = 2*(j-1)/(nv-1); - X(h,1) = x; - X(h,2) = y; - %X(h,3) = cos(10*pi*x)*cos(3*pi*y)+exp(2*sin(x*y)); - %X(h,3) = cos(3*pi*exp(-x))*cos(3*pi*y^2)+exp(2*sin(x*y)); - X(h,3) = 3*x*y; %- 0.5*rand(1,1) +10;%+exp(2*sin(x*y)); + 0.01*exp(x*y)+ - X(h,4) = 1; - end -end - -end - diff --git a/interpolate_intersections.m b/interpolate_intersections.m deleted file mode 100644 index 6089bbb..0000000 --- a/interpolate_intersections.m +++ /dev/null @@ -1,255 +0,0 @@ -%close all - -spline_surf_vec -intersections - - plot3(GXs,GYs,5*ones(us_n_basf+1,vs_n_basf+1)); - hold on - plot3(GYs,GXs,5*ones(us_n_basf+1,vs_n_basf+1)); - - for i=1:m - plot3(point(i,5),point(i,6),point(i,7),'.','MarkerSize',50); - end - -%pointx = point; -pointb = point; - - -%%% sort intersections by patches (1st patch) - - sp_i = (pointb(:,10)-1)*(us_n_intervs-1) + pointb(:,11); - -[ssp_i idx] = sort(sp_i); - -m = length(sp_i); %m-out - - - -pointb = pointb(idx,:); - -%%% Compute intersectioned patches -patch_int(1,2) = 1; -patch_int(1,1) = ssp_i(1); -%a = sp_i(1); -different = 1; -for i =2:m - if ssp_i(i) == patch_int(different,1)%a - patch_int(different,2) = patch_int(different,2) +1; - continue - else - %a = sp_i(i); - different = different +1; - patch_int(different,1) = ssp_i(i); - patch_int(different,2) = 1; - end -end - - -different; - - - -%return - - - -coinf = zeros(m,2); - - -a = 1; %?? - -%%% detect intersection point types -% -1 - interion -% 0 - global boundary -% 1 - patch boundary - -for j=1:m - coi = boundary_point(pointb(j,3:4),pointb(j,10:11),us_knots,vs_knots); - coinf(j,:) = coi; % type of interrsection (u,v) % 1D -end - -out = zeros(different,1); -point_type = zeros(m,1); % 2D - -offset = 0; -for j=1:different - for i = 1:patch_int(j,2) - - if coinf(i+offset,1) > -1 || coinf(i+offset,2) > -1 % number of boundary points for patch - out(j) = out(j) +1; - end - - % define intersection type - if coinf(i+offset,1) == -1 && coinf(i+offset,2) == -1 % internal - point_type(i+offset) = -1; - elseif coinf(i+offset,1) == 0 || coinf(i+offset,2) == 0 % global boundary - point_type(i+offset) = 0; - elseif coinf(i+offset,1) == 1 || coinf(i+offset,2) == 1 % patch internal boundary - point_type(i+offset) = 1; - end - - - - end - offset = offset + patch_int(j,2); -end - - -% number of outputs - - - -% % Sort points in patch -% offset = 0; -% for j=1:different -% -% -% patch_points= offset:offset+patch_int(j,2); -% -% boundg = find(point_type(patch_points) == 0) -% boundi = find(point_type(patch_points) == 1) -% bound = [boundg , boundi]; -% -% l = lenght(bound) -% if l >2 -% disp('more then 2 boundary points') -% end -% -% dist = zeros(patch_int(j,2)); -% -% for k = 1: offset+patch_int(j,2) -% for i = 1: offset+patch_int(j,2) -% -% dist(i) = norm(pointsb(patch_points(k),1:2)- pointsb(patch_points(i),1:2)) -% i+offset -% -% end -% -% end -% offset = offset + patch_int(j,2); -% end -% - -% Sort points in surface (create components) - - offset = 0; - bound = find(point_type == 0); - pointb = [pointb coinf point_type]; - [a,b] = size(pointb); - splinepoint = pointb; - splinepoint = swaplines(splinepoint,1,bound(1)); - - for j=1:m-1 - dist = 1/eps * ones(m,1); - for k=j+1:m - dist(k) = norm(splinepoint(j,1:2) - splinepoint(k,1:2)); - end - [a b] = min(dist); - splinepoint = swaplines(splinepoint,b,j+1); - end - - % Sort intersection on patches - boundg = find(point_type == 0); - boundl = find(point_type == 1); - bound = [boundg;boundl] - -% offset =0; -% for k=1:different -% boundg = find(point_type(1+offset:patch_int(different,2)+offset)== 0); -% boundl = find(point_type(1+offset:patch_int(different,2)+offset)== 1); -% bound = [boundg;boundl] -% %pause -% pointb = swaplines(pointb,bound(1)+offset,1+offset) -% for l=1:patch_int(k,2)-1 -% dist = 1/eps * ones(patch_int(k,2)-l,1); -% for j=1:patch_int(k,2)-l -% dist(j) = norm(splinepoint(l,1:2) - splinepoint(l+j,1:2)); -% end -% [a b] = min(dist); -% splinepoint = swaplines(splinepoint,b+offset,l+1+offset); -% -% end -% offset = offset + patch_int(k,2) -% end - - - - - %%%% Remove duplicite points () - - matrixpoints = ones(m,1); - a = 0; - for i=1:m -% if splinepoint(i,14) == 1 -% a = 1; -% end - if splinepoint(i,14)*a == 1 - matrixpoints(i) = 0; - a = 0; - elseif splinepoint(i,14) == 1 - a =1; - else - a = 0; - end - end - - -% Interpolate points -% -% -% pointb = pointb(1:m-out,:) -% -% mm = m - out; -% -% %return -% -mr = sum(matrixpoints); -XYZ = zeros(mr,3); -deltas = zeros(mr,1); -j = 0; - -for i=1:m - if matrixpoints(i)==1 - j = j+1; - XYZ(j,:) =splinepoint(i,5:7); - if j~=1 - deltas(j) = norm(XYZ(j,:)-XYZ(j-1,:)); - end - end -end - -Deltas = zeros(mr,1); - -for i=2:mr - Deltas(i)=sum(deltas(1:i)); -end - -pointspersegment=4; -nbasf = ceil(mr/pointspersegment); -t_knots = get_knot_vector(nbasf+2); -t = Deltas/Deltas(mr); -X = zeros(mr,nbasf+2); -for i=1:mr - - [tf, ~] = splinebasevec(t_knots,t(i),0); - X(i,:) = tf; -end - -sol = X\XYZ; -X*sol - XYZ; - -%X'*X\X'*XYZ - -ns = 100; -td = linspace(0,1,ns); -for i=1:ns - PT(i,:)=splinebasevec(t_knots,td(i),0)'*sol; - -end -plot3(PT(:,1),PT(:,2),PT(:,3),'r','LineWidth',3); -% -% -% -% -% -% \ No newline at end of file diff --git a/intersections.m b/intersections.m deleted file mode 100644 index 82bf21b..0000000 --- a/intersections.m +++ /dev/null @@ -1,192 +0,0 @@ - -% spline_surf_vec - %close all -% % % -% plotresultsvec(u_knots, v_knots,P0,P1,P2,P3,X,z,Err) -% hold on -% plotresultsvec(us_knots, vs_knots,P0s,P1s,P2s,P3s,Xs,zs,Errs) -hold on - - -%%%%%%%%%%%%%%%%%%%%%%%%% -%%% Compute intersections -%%%%%%%%%%%%%%%%%%%%%%%%% - -% u_n_basf = length(u_knots)-3; -% v_n_basf = length(v_knots)-3; -% us_n_basf = length(us_knots)-3; -% vs_n_basf = length(vs_knots)-3; - -u_n_intervs = u_n_basf - 2; -v_n_intervs = v_n_basf - 2; -us_n_intervs = us_n_basf - 2; -vs_n_intervs = vs_n_basf - 2; - -u_n_grid = u_n_basf+1; -v_n_grid = v_n_basf+1; -us_n_grid = us_n_basf+1; -vs_n_grid = vs_n_basf+1; - -% u_n_intervs = u_n_basf - 2; -% v_n_intervs = v_n_basf - 2; -% us_n_intervs = us_n_basf - 2; -% vs_n_intervs = vs_n_basf - 2; - -% Compute X & Y grid intersections - -%%% Grid Boundary -[ GX,GY ] = compute_grid_lines(u_knots,v_knots, P0,P1,P2,P3 ); -[ GXs,GYs ] = compute_grid_lines(us_knots,vs_knots, P0s,P1s,P2s,P3s ); - -%%% Grid centers -[ X_coor,Y_coor] = compute_control_points( u_knots, v_knots, P0,P1,P2,P3); -[ Xs_coor,Ys_coor] = compute_control_points( us_knots, vs_knots, P0s,P1s,P2s,P3s); - -%%% Compute Bounding Boxes -Z_coor = vec2mat( z,u_n_basf,v_n_basf ); -Zs_coor = vec2mat( zs,us_n_basf,vs_n_basf ); - -[patch_bound_X,patch_bound_Y] = compute_patch_edges( u_knots, v_knots, P0,P1,P2,P3); -[patch_bound_Xs,patch_bound_Ys] = compute_patch_edges( us_knots, vs_knots, P0s,P1s,P2s,P3s); - -% [ BB_X,BB_Y,BB_Z ] = compute_bounding_box( X_coor,Y_coor,Z_coor, u_n_intervs,v_n_intervs); -% [ BB_Xs,BB_Ys,BB_Zs ] = compute_bounding_box( Xs_coor,Ys_coor,Zs_coor, us_n_intervs,vs_n_intervs); - -%%% Bonding boxes intersections -[isec,n_isec] = bounding_boxes_intersection( patch_bound_X,patch_bound_Y,Z_coor,patch_bound_Xs,patch_bound_Ys,Zs_coor); - - -x = X_coor(:); -y = Y_coor(:); -xs = Xs_coor(:); -ys = Ys_coor(:); - -Xs = [xs ys zs]; - -X = [x y z]; -nt = 6; % number of main lines -n_points = 0; -ninter = zeros(u_n_intervs,v_n_intervs); -for k=1:u_n_intervs - us1 = u_knots(k+2); - ue1 = u_knots(k+3); - u1_c =(us1 + ue1)/2; - ui = [us1 ue1] ; - for l=1:v_n_intervs - vs1 = v_knots(l+2); - ve1 = v_knots(l+3); - v1_c = (vs1 + ve1)/2; - vi = [ vs1 ve1]; - s=0; - if n_isec(k,l) ~= 0 - for p =1: n_isec(k,l) - - m = ceil(isec(k,l,p) / us_n_intervs); - o = isec(k,l,p) - (m-1)*us_n_intervs; - sp_i = (m-1)*(us_n_intervs) + o; - % v2 fixed - u2i = [us_knots(m+2) us_knots(m+3)]; - v_knot = linspace(vs_knots(o+2),vs_knots(o+3),nt); - for h =1:length(v_knot) - u2_c = (us_knots(m+2) + us_knots(m+3))/2; - v2i = v_knot(h); - nit = 10; - uvu2v2 = [u1_c;v1_c;u2_c]; - [ uvu2v2,ptl,conv ] = patch_patch_intersection( uvu2v2, ui,vi,u2i,v2i, u_knots, v_knots, X,us_knots, vs_knots, Xs, nit ,m,o,k,l); - if conv ~= 0 - s = s+1; - plot3(ptl(1),ptl(2),ptl(3),'.','MarkerSize',50) - n_points = n_points +1; - point(n_points,:) = [uvu2v2',v_knot(h),ptl(1),ptl(2),ptl(3),k,l,m,o]; - end - end - - % u2 fixed - v2i = [vs_knots(o+2) vs_knots(o+3)]; - u_knot = linspace(us_knots(m+2),us_knots(m+3),nt); - for h =1:length(u_knot) - v2_c = (vs_knots(o+2) + vs_knots(o+3))/2; - u2i = u_knot(h); - nit = 10; - uvu2v2 = [u1_c;v1_c;v2_c]; - [ uvu2v2,ptl,conv ] = patch_patch_intersection( uvu2v2, ui,vi,u2i,v2i, u_knots, v_knots, X,us_knots, vs_knots, Xs, nit ,m,o,k,l); - if conv ~= 0 - s = s+1; - plot3(ptl(1),ptl(2),ptl(3),'.','MarkerSize',50) - n_points = n_points +1; - point(n_points,:) = [uvu2v2(1:2)',u_knot(h),uvu2v2(3),ptl(1),ptl(2),ptl(3),k,l,m,o]; - end - end - ninter(k,l) = s; - end - end - end -end - - ninter -%[isec,n_isec] = bounding_boxes_intersection2( patch_bound_X,patch_bound_Y,Z_coor,patch_bound_Xs,patch_bound_Ys,Zs_coor); -%[isec2,n_isec2] = bounding_boxes_intersection2( patch_bound_Xs,patch_bound_Ys,Zs_coor,patch_bound_X,patch_bound_Y,Z_coor); - -% for m=1:us_n_intervs -% us2 = us_knots(m+2); -% ue2 = us_knots(m+3); -% u2_c = (us2+ue2)/2; -% ui2 = [us2 ue2] ; -% for o=1:vs_n_intervs -% vs2 = vs_knots(o+2); -% ve2 = vs_knots(o+3); -% v2_c = (vs2+ve2)/2; -% vi2 = [ vs2 ve2]; -% s = 0; -% if n_isec2(m,o) ~= 0 -% for p =1: n_isec2(m,o) -% k = ceil(isec2(m,o,p) / u_n_intervs); -% l = isec2(m,o,p) - (k-1)*u_n_intervs; -% sp_i = (k-1)*(u_n_intervs) + l; -% % [sp_i isec2(m,o,p)] -% % pause -% t0 = 0.5; -% -% % v fixed -% v_knot = linspace(v_knots(l+2),v_knots(l+3),nt); -% for h =1:length(v_knot) -% u_val = linspace(u_knots(k+2),u_knots(k+3),3); -% v_val = v_knot(h); -% XYZ = get_span(u_knots,v_knots,X,u_val,v_val); -% nit = 6; -% uvt0 = [u2_c;v2_c;t0]; -% [ uvtl, ptl ,conv] = compute_initial_guess( uvt0, ui2,vi2,us_knots, vs_knots, Xs, XYZ, nit,m,o); -% if conv ~= 0 -% s = s+1; -% plot3(ptl(1),ptl(2),ptl(3),'.','MarkerSize',50) -% end -% end -% -% % u fixed -% u_knot = linspace(u_knots(k+2),u_knots(k+3),nt); -% for h =1:length(u_knot) -% u_val = u_knot(h); -% v_val = linspace(v_knots(l+2),v_knots(l+3),3); -% XYZ = get_span(u_knots,v_knots,X,u_val,v_val); -% nit = 6; -% uvt0 = [u2_c;v2_c;t0]; -% [ uvtl, ptl ,conv] = compute_initial_guess( uvt0, ui2,vi2,us_knots, vs_knots, Xs, XYZ, nit,m,o); -% if conv ~= 0 -% s = s+1; -% plot3(ptl(1),ptl(2),ptl(3),'.','MarkerSize',50) -% end -% end -% -% ninter(m,o) = ninter(m,o) + s; -% end -% end -% end -% end - -%t = toc - -%cas = sum(sum(ninter))/t - -ninter - plot3(GXs,GYs,5*ones(us_n_basf+1,vs_n_basf+1)) - plot3(GYs,GXs,5*ones(us_n_basf+1,vs_n_basf+1)) diff --git a/patch_patch_intersection.m b/patch_patch_intersection.m deleted file mode 100644 index b39ea02..0000000 --- a/patch_patch_intersection.m +++ /dev/null @@ -1,71 +0,0 @@ -function [ uvt,pt,conv ] = patch_patch_intersection( uvt, ui,vi,u2i,v2i,u_knots, v_knots, X,us_knots, vs_knots, Xs, nit ,m,o,k,l) - -pt =zeros(3,1); -conv =0; -tol = 1e-4; -tol2 = 1e-6; - -if length(u2i) == 1 - [u2f, ~] = splinebasevec(us_knots,u2i,0,m); -end -if length(v2i) == 1 - [v2f, ~] = splinebasevec(vs_knots,v2i,0,o); -end - - -for i=1:nit - [uf, ~] = splinebasevec(u_knots,uvt(1),0,k); - [vf, ~] = splinebasevec(v_knots,uvt(2),0,l); - [ufd, ~] = splinebasevec(u_knots,uvt(1),1,k); - [vfd, ~] = splinebasevec(v_knots,uvt(2),1,l); - - if length(u2i) == 1 - [v2f, ~] = splinebasevec(vs_knots,uvt(3),0,o); - [v2fd, ~] = splinebasevec(vs_knots,uvt(3),1,o); - dXYZp2 = (kron(v2fd',u2f')*Xs)'; - end - if length(v2i) == 1 - [u2f, ~] = splinebasevec(us_knots,uvt(3),0,m); - [u2fd, ~] = splinebasevec(us_knots,uvt(3),1,m); - dXYZp2 = (kron(v2f',u2fd')*Xs)'; - end - - dXYZu1 = (kron(vf',ufd')*X)'; % - dXYZv1 = (kron(vfd',uf')*X)'; % - - J = [dXYZu1 dXYZv1 -dXYZp2]; - - deltaXYZ = (kron(vf',uf') * X)' - (kron(v2f',u2f') * Xs)'; - uvt = uvt- J\deltaXYZ; - %if i~=nit - [test,uvt] = rangetest(uvt,ui,vi,u2i,v2i,0.0); - %end -end - -[test,uvt] = rangetest(uvt,ui,vi,u2i,v2i,tol2); -if test == 1 - [uf, ~] = splinebasevec(u_knots,uvt(1),0,k); - [vf, ~] = splinebasevec(v_knots,uvt(2),0,l); - - if length(u2i) == 1 - [v2f, ~] = splinebasevec(vs_knots,uvt(3),0,o); - end - if length(v2i) == 1 - [u2f, ~] = splinebasevec(us_knots,uvt(3),0,m); - end - - deltaXYZ = (kron(vf',uf') * X)' - (kron(v2f',u2f') * Xs)'; - dist = norm(deltaXYZ); -end - -if test == 1 - if dist <= tol - pt =kron(vf',uf')*X; - conv =1; - end -else - uvt = zeros(3,1); -end - -end - diff --git a/plotresultsvec.m b/plotresultsvec.m deleted file mode 100644 index 7d510b8..0000000 --- a/plotresultsvec.m +++ /dev/null @@ -1,73 +0,0 @@ -function [ ] = plotresultsvec( u_knots, v_knots,P0,P1,P2,P3,X,z, Err) - -[np k] = size(X); - -u_n_basf = length(u_knots)-3; -v_n_basf = length(v_knots)-3; - - -nu = 60; -nv = 60; - -Zsurf = zeros(nu,nv); -Xsurf = zeros(nu,nv); -Ysurf = zeros(nu,nv); - -% Compute X & Y control points - - -[ X_coor,Y_coor] = compute_control_points( u_knots, v_knots, P0,P1,P2,P3); - -x = X_coor(:); -y = Y_coor(:); - - -% Compute fine grid in X & Y & Z for draw - -uf = spalloc(u_n_basf,nu,3*nu); -vf = spalloc(v_n_basf,nv,3*nv); - -for j =1:nu - u = (j-1)/(nu-1); - uf(:,j) = splinebasevec(u_knots,u,0); -end - -for k =1:nv - v = (k-1)/(nv-1); - vf(:,k) = splinebasevec(v_knots,v,0); -end - -for k =1:nv - Zsurf(:,k) = kron(vf',uf(:,k)') * z; - Xsurf(:,k) = kron(vf',uf(:,k)') * x; - Ysurf(:,k) = kron(vf',uf(:,k)') * y; -end - - - -surf(Xsurf,Ysurf,Zsurf); - - -%%plot original points - -% hold on -% for k=1:np -% if X(k,4) ~=0 -% plot3(X(k,1),X(k,2),X(k,3),'.k','MarkerSize',25); -% end -% end - -%% - -% hold off -% -% figure -% -% Xsurferr = kron(xv(2:u_n_basf-1)',ones(v_n_basf-2,1)); -% Ysurferr = kron(yv(2:v_n_basf-1),ones(u_n_basf-2,1)'); -% -% -% surf(Xsurferr,Ysurferr,Err); - -end - diff --git a/rangetest.m b/rangetest.m deleted file mode 100644 index af8ef7c..0000000 --- a/rangetest.m +++ /dev/null @@ -1,51 +0,0 @@ -function [ test, uvt ] = rangetest( uvt,ui,vi,u2i,v2i,tol ) - - -test = 0; - -du = [uvt(1)-ui(1), ui(2)- uvt(1)]; -dv = [uvt(2)-vi(1), vi(2)- uvt(2)]; - -if length(v2i) == 1 -d2p = [uvt(3)-u2i(1), u2i(2)- uvt(3)]; -pi = u2i; -end - -if length(u2i) == 1 -d2p = [uvt(3)-v2i(1), v2i(2)- uvt(3)]; -pi = v2i; -end - -%dv = [uvt(2)-vi(1), vi(2)- uvt(2)]; - - -for i =1:2 - if (du(i) < -tol) - uvt(1) = ui(i); - end -end - -for i =1:2 - if (dv(i) < -tol) - uvt(2) = vi(i); - end -end - -for i =1:2 - if (d2p(i) < -tol) - uvt(3) = pi(i); - end -end - - -if (uvt(1)>=ui(1)) && (uvt(1)<=ui(2)) - if (uvt(2)>=vi(1)) && (uvt(2)<=vi(2)) - if ((uvt(3)>=pi(1)) && (uvt(3)<=pi(2))) - test = 1; - end - end -end - - -end - diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..e69de29 diff --git a/solve.m b/solve.m deleted file mode 100644 index 7f1d395..0000000 --- a/solve.m +++ /dev/null @@ -1,38 +0,0 @@ -function [ Xp ] = solve( Xp,P0,P1,P2,P3 ) - -[np k] = size(Xp); -Pi = [P0 P1 P2 P3]; - -%%% Compute local coordinates and drop all - -A = P3 - P2; -B = P0 - P1; -C = P1 - P2; -D = P0 - P3; - -for j=1:np - - P = [Xp(j,5); Xp(j,6)]; - - u = Xp(j,1); - v = Xp(j,2); - - uv = [u;v]; - - iters = 10; - - for i = 1:iters - u = uv(1); - v = uv(2); - - % J = [ v * A(1) + (1-v) * B(1), u * C(1) + (1 - u) * D(1); - % v * A(2) + (1-v) * B(2), u * C(2) + (1 - u) * D(2)]; - J = [ v * A + (1-v) * B, u * C + (1 - u) * D]; - Fxy = P - u * (v * P2 + (1-v) * P1) - (1-u) * ( v * P3 + (1-v) * P0); - uv = [u;v] - inv(J) * Fxy; - end - Xp(j,1) = uv(1); - Xp(j,2) = uv(2); -end -end - diff --git a/spline_surf_vec.m b/spline_surf_vec.m deleted file mode 100644 index 5d067a8..0000000 --- a/spline_surf_vec.m +++ /dev/null @@ -1,160 +0,0 @@ - -clc -clear all -close all - - -u_n_basf = 15; -v_n_basf = 15; -u_knots = get_knot_vector(u_n_basf); -v_knots = get_knot_vector(v_n_basf); - -%%% base test - -% basetestvec(v_knots); -% return -%%% patch boundary - -% P0 = [0.2; 0.3]; -% P1 = [1.2; 0.5]; -% P2 = [1.5; 2]; -% P3 = [0.9; 1.7]; - -P0 = [0.0; 0.0]; -P1 = [2.0; 0.0]; -P2 = [2.0; 2.0]; -P3 = [0.0; 2.0]; - -%%% Reduce points to be approximatex (i.e., points which are in convex hull of the points P0,P1,P2,P3) - -X = getpoints(); -[Xp X] = transform( X,P0,P1,P2,P3 ); -[np k] = size(Xp); -[ Xp ] = solve( Xp,P0,P1,P2,P3 ); - -%%% Construction of the matrix - - -%interv ?? -[B, Interv ] = build_LS_matrix( u_knots,v_knots, Xp); - -W = sparse(diag(Xp(:,4))); - -g = Xp(:,3); -b = B'*W*g; - -C = B'*W*B; -nnzC = nnz(C); - - -A = build_reg_matrix( u_knots,v_knots, P0,P1,P2,P3,nnzC); - -nC = norm(full(C)); -nA = norm(full(A)); -r = norm(full(C))/ norm(full(A)); -%r = 0.0; - -% [q r] = qr(B); -% -% -% z = r\(q'*g) - -S = C+0.01*r*A; - -z = pcg(S,b,1e-12,500); - -%%% Solution - -% Direct Solver -% C = B'*W*B;%+A; -% b = B'*W*g; -% z = C\b; - - -% Errors -Err = get_errors(abs(W*B*z-g),Interv,u_n_basf,v_n_basf); - - -%%% PLOT results - -plotresultsvec(u_knots, v_knots,P0,P1,P2,P3,X,z,Err) - -%return - -%%%%%%%%%%%%%%%%%% -%%% Second Surface -%%%%%%%%%%%%%%%%%% - - -us_n_basf = 10; -vs_n_basf = 10; -us_knots = get_knot_vector(us_n_basf); -vs_knots = get_knot_vector(vs_n_basf); - -%%% patch boundary - -% P0s = [0.4; 0.2]; -% P1s = [1.5; 0.1]; -% P2s = [1.1; 1.8]; -% P3s = [0.6; 1.7]; - -P0s = [0.0; 0.0]; -P1s = [2.0; 0.0]; -P2s = [2.0; 2.0]; -P3s = [0.0; 2.0]; - - -%%% Reduce points to be approximatex (i.e., points which are in convex hull of the points P0,P1,P2,P3) - -Xs = getpoints2(); -[Xps Xs] = transform( Xs,P0s,P1s,P2s,P3s ); -[nps ks] = size(Xps); -[ Xps ] = solve( Xps,P0s,P1s,P2s,P3s ); - -%%% Construction of the matrix - -%interv ?? -[Bs, Intervs ] = build_LS_matrix( us_knots,vs_knots, Xps); - -Ws = sparse(diag(Xps(:,4))); - -gs = Xps(:,3); -bs = Bs'*Ws*gs; - -Cs = Bs'*Ws*Bs; -nnzCs = nnz(Cs); - - -As = build_reg_matrix( us_knots,vs_knots, P0s,P1s,P2s,P3s,nnzCs); - -nCs= norm(full(Cs)); -nAs = norm(full(As)); -rs = norm(full(Cs))/ norm(full(As)); -%rs = 0.0; - -Ss = Cs+0.0010*rs*As; - -zs = pcg(Ss,bs,1e-14,500); -%zs = Ss\bs; - -%%% Solution - -% Direct Solver -% C = B'*W*B;%+A; -% b = B'*W*g; -% z = C\b; - - -% Errors -Errs = get_errors(abs(Ws*Bs*zs-gs),Intervs,us_n_basf,vs_n_basf); - - -%%% PLOT results -hold on -plotresultsvec(us_knots, vs_knots,P0s,P1s,P2s,P3s,Xs,zs,Errs) -hold off - - -%intersections - -% \ No newline at end of file diff --git a/splinebasevec.m b/splinebasevec.m deleted file mode 100644 index 007e303..0000000 --- a/splinebasevec.m +++ /dev/null @@ -1,80 +0,0 @@ -function [ f, k ] = splinebasevec( T,t,order,varargin) - -% T - knot vector -% k - index of basis function, k = 1,...,length(T)-3 -% t - parameter t \in [0,1] -% f - vector of basis functions values -% k - kth interval - -n = length(T); - -f = spalloc(n-3,1,3); - -optargin = length(varargin); - -if optargin == 0 - k = find_int( T,t ); -elseif optargin == 1 - k = varargin{1}; -% if (T(k) <= t+2) && (T(k+3) >= t) -% disp('OK') -% [T(k+2), t,T(k+3)] -% else -% disp('PROBLEM') -% return -% end -if (T(k) <= t+2) && (T(k+3) >= t) - % disp('OK'); - [T(k+2), t,T(k+3)]; -else - disp('PROBLEM') - pause - return -end -end - - -% for k = 1:n-5 -% if(t>=T(k+2) && t<=T(k+3)) -% % k_int = k; -% break; -% end -% end -% -% if k ~=k2-2 -% [k k2-2] -% pause -% end - -tk1 = T(k+1); -tk2 = T(k+2); -tk3 = T(k+3); -tk4 = T(k+4); - -d31 = tk3-tk1; -d32 = tk3-tk2; -d42 = tk4-tk2; - -d3t = tk3-t; -dt1 = t-tk1; -dt2 = t-tk2; -d4t = tk4-t; - -d31d32 = d31*d32; -d42d32 = d42*d32; - -% f(k) = (tk3-t)^2/((tk3-tk1)*(tk3-tk2)); -% f(k+1)= (t-tk1)*(tk3 -t)/((tk3-tk1)*(tk3-tk2)) + (t-tk2)*(tk4 -t)/((tk4-tk2)*(tk3-tk2)); -% f(k+2) = (t-tk2)^2/((tk4 - tk2)*(tk3-tk2)); - -if order == 0 - f(k) = d3t^2/d31d32; - f(k+1)= dt1*d3t/d31d32 + dt2*d4t/d42d32; - f(k+2) = dt2^2/d42d32; -elseif order ==1 - f(k) = -2*d3t/d31d32; - f(k+1)= (d3t - dt1)/d31d32 + (d4t - dt2)/d42d32; - f(k+2) = 2*dt2/d42d32; -end - -end \ No newline at end of file diff --git a/src/brep_writer.py b/src/brep_writer.py new file mode 100644 index 0000000..7deb770 --- /dev/null +++ b/src/brep_writer.py @@ -0,0 +1,1091 @@ +import enum +import numpy as np + + +''' +TODO: +- For Solid make auto conversion from Faces similar to Face from Edges +- For Solid make test that Shells are closed. +- Implement closed test for shells (similar to wire) +- Improve test of slosed (Wire, Shell) to check also orientation of (Edges, Faces). ?? May be both if holes are allowed. +- Rename attributes and methods to more_words_are_separated_by_underscore. +- rename _writeformat to _brep_output +- remove (back) groups parameter from _brepo_output, make checks of IDs at main level (only in DEBUG mode) +- document public methods +''' + + +class ParamError(Exception): + pass + +def check_matrix(mat, shape, values, idx=[]): + ''' + Check shape and type of scalar, vector or matrix. + :param mat: Scalar, vector, or vector of vectors (i.e. matrix). Vector may be list or other iterable. + :param shape: List of dimensions: [] for scalar, [ n ] for vector, [n_rows, n_cols] for matrix. + If a value in this list is None, the dimension can be arbitrary. The shape list is set fo actual dimensions + of the matrix. + :param values: Type or tuple of allowed types of elements of the matrix. E.g. ( int, float ) + :param idx: Internal. Used to pass actual index in the matrix for possible error messages. + :return: + ''' + try: + + if len(shape) == 0: + if not isinstance(mat, values): + raise ParamError("Element at index {} of type {}, expected instance of {}.".format(idx, type(mat), values)) + else: + + if shape[0] is None: + shape[0] = len(mat) + l=None + if not hasattr(mat, '__len__'): + l=0 + elif len(mat) != shape[0]: + l=len(mat) + if not l is None: + raise ParamError("Wrong len {} of element {}, should be {}.".format(l, idx, shape[0])) + for i, item in enumerate(mat): + sub_shape = shape[1:] + check_matrix(item, sub_shape, values, idx = [i] + idx) + shape[1:] = sub_shape + return shape + except ParamError: + raise + except Exception as e: + raise ParamError(e) + + +class Location: + """ + Location defines an affine transformation in 3D space. Corresponds to the in the BREP file. + BREP format allows to use different transformations for individual shapes. + Location are numberd from 1. Zero index means identity location. + """ + def __init__(self, matrix=None): + """ + Constructor for elementary afine transformation. + :param matrix: Transformation matrix 3x4. First three columns forms the linear transformation matrix. + Last column is the translation vector. + matrix==None means identity location (ID=0). + """ + if matrix is None: + self.matrix = None + self.id = 0 + return + + # checks + check_matrix(matrix, [3, 4], (int, float)) + self.matrix=matrix + + def _dfs(self, groups): + """ + Deep first search that assign numbers to shapes + :param groups: dict(locations=[], curves_2d=[], curves_3d=[], surfaces=[], shapes=[]) + :return: None + """ + if not hasattr(self, 'id'): + id=len(groups['locations'])+1 + self.id = id + groups['locations'].append(self) + + def _brep_output(self, stream, groups): + stream.write("1\n") + for row in self.matrix: + for number in row: + stream.write(" {}".format(number)) + stream.write("\n") + +class ComposedLocation(Location): + """ + Defines an affine transformation as a composition of othr transformations. Corresponds to the in the BREP file. + BREP format allows to use different transformations for individual shapes. + """ + def __init__(self, location_powers=[]): + """ + + :param location_powers: List of pairs (location, power) where location is instance of Location and power is float. + """ + locs, pows = zip(*location_powers) + l = len(locs) + check_matrix(locs, [ l ], (Location, ComposedLocation) ) + check_matrix(pows, [ l ], int) + + self.location_powers=location_powers + + + def _dfs(self, groups): + + for location,power in self.location_powers: + location._dfs(groups) + Location._dfs(self, groups) + + def _brep_output(self, stream, groups): + stream.write("2 ") + for loc, pow in self.location_powers: + stream.write("{} {} ".format(loc.id, pow)) + stream.write("0\n") + +def check_knots(deg, knots, N): + total_multiplicity = 0 + for knot, mult in knots: + # This condition must hold if we assume only (0,1) interval of curve or surface parameters. + #assert float(knot) >= 0.0 and float(knot) <= 1.0 + total_multiplicity += mult + assert total_multiplicity == deg + N + 1 + + +# TODO: perform explicit conversion to np.float64 in order to avoid problems on different arch +# should be unified in bspline as well, convert to np.arrays as soon as posible +scalar_types = (int, float, np.int64, np.float64) + + +def curve_from_bs( curve): + """ + Make BREP writer Curve (2d or 3d) from bspline curve. + :param curve: bs.Curve object + :return: + """ + dim = curve.dim + if dim == 2: + curve_dim = Curve2D + elif dim == 3: + curve_dim = Curve3D + else: + assert False + c = curve_dim(curve.poles, curve.basis.pack_knots(), curve.rational, curve.basis.degree) + c._bs_curve = curve + return c + +class Curve3D: + """ + Defines a 3D curve as B-spline. We shall work only with B-splines of degree 2. + Corresponds to "B-spline Curve - <3D curve record 7>" from BREP format description. + """ + def __init__(self, poles, knots, rational=False, degree=2): + """ + Construct a B-spline in 3d space. + :param poles: List of poles (control points) ( X, Y, Z ) or weighted points (X,Y,Z, w). X,Y,Z,w are floats. + Weighted points are used only for rational B-splines (i.e. nurbs) + :param knots: List of tuples (knot, multiplicity), where knot is float, t-parameter on the curve of the knot + and multiplicity is positive int. Total number of knots, i.e. sum of their multiplicities, must be + degree + N + 1, where N is number of poles. + :param rational: True for rational B-spline, i.e. NURB. Use weighted poles. + :param degree: Positive int. + """ + + if rational: + check_matrix(poles, [None, 4], scalar_types ) + else: + check_matrix(poles, [None, 3], scalar_types) + N = len(poles) + check_knots(degree, knots, N) + + self.poles=poles + self.knots=knots + self.rational=rational + self.degree=degree + + def _eval_check(self, t, point): + if hasattr(self, '_bs_curve'): + repr_pt = self._bs_curve.eval(t) + if not np.allclose(np.array(point), repr_pt , rtol = 1.0e-3): + raise Exception("Point: {} far from curve repr: {}".format(point, repr_pt)) + + + def _dfs(self, groups): + if not hasattr(self, 'id'): + id = len(groups['curves_3d']) + 1 + self.id = id + groups['curves_3d'].append(self) + + + def _brep_output(self, stream, groups): + # writes b-spline curve + stream.write("7 {} 0 {} {} {} ".format(int(self.rational), self.degree, len(self.poles), len(self.knots))) + for pole in self.poles: + for value in pole: + stream.write(" {}".format(value)) + stream.write(" ") + for knot in self.knots: + for value in knot: + stream.write(" {}".format(value)) + stream.write(" ") + stream.write("\n") + +class Curve2D: + """ + Defines a 2D curve as B-spline. We shall work only with B-splines of degree 2. + Corresponds to "B-spline Curve - <2D curve record 7>" from BREP format description. + """ + def __init__(self, poles, knots, rational=False, degree=2): + """ + Construct a B-spline in 2d space. + :param poles: List of points ( X, Y ) or weighted points (X,Y, w). X,Y,w are floats. + Weighted points are used only for rational B-splines (i.e. nurbs) + :param knots: List of tuples (knot, multiplicity), where knot is float, t-parameter on the curve of the knot + and multiplicity is positive int. Total number of knots, i.e. sum of their multiplicities, must be + degree + N + 1, where N is number of poles. + :param rational: True for rational B-spline, i.e. NURB. Use weighted poles. + :param degree: Positive int. + """ + + N = len(poles) + if rational: + check_matrix(poles, [N, 3], scalar_types ) + else: + check_matrix(poles, [N, 2], scalar_types) + check_knots(degree, knots, N) + + self.poles=poles + self.knots=knots + self.rational=rational + self.degree=degree + + def _eval_check(self, t, surface, point): + if hasattr(self, '_bs_curve'): + u, v = self._bs_curve.eval(t) + surface._eval_check(u, v, point) + + + def _dfs(self, groups): + if not hasattr(self, 'id'): + id = len(groups['curves_2d']) + 1 + self.id = id + groups['curves_2d'].append(self) + + def _brep_output(self, stream, groups): + # writes b-spline curve + stream.write("7 {} 0 {} {} {} ".format(int(self.rational), self.degree, len(self.poles), len(self.knots))) + for pole in self.poles: + for value in pole: + stream.write(" {}".format(value)) + stream.write(" ") + for knot in self.knots: + for value in knot: + stream.write(" {}".format(value)) + stream.write(" ") + stream.write("\n") + + + +def surface_from_bs(surf): + """ + Make BREP writer Surface from bspline surface. + :param surf: bs.Surface object + :return: + """ + s = Surface(surf.poles, (surf.u_basis.pack_knots(), surf.v_basis.pack_knots()), + surf.rational, (surf.u_basis.degree, surf.v_basis.degree) ) + s._bs_surface = surf + return s + +class Surface: + """ + Defines a B-spline surface in 3d space. We shall work only with B-splines of degree 2. + Corresponds to "B-spline Surface - < surface record 9 >" from BREP format description. + """ + def __init__(self, poles, knots, rational=False, degree=(2,2)): + """ + Construct a B-spline in 3d space. + :param poles: Matrix (list of lists) of Nu times Nv poles (control points). + Single pole is a points ( X, Y, Z ) or weighted point (X,Y,Z, w). X,Y,Z,w are floats. + Weighted points are used only for rational B-splines (i.e. nurbs) + :param knots: Tuple (u_knots, v_knots). Both u_knots and v_knots are lists of tuples + (knot, multiplicity), where knot is float, t-parameter on the curve of the knot + and multiplicity is positive int. For both U and V knot vector the total number of knots, + i.e. sum of their multiplicities, must be degree + N + 1, where N is number of poles. + :param rational: True for rational B-spline, i.e. NURB. Use weighted poles. BREP format have two independent flags + for U and V parametr, but only choices 0,0 and 1,1 have sense. + :param degree: (u_degree, v_degree) Both positive ints. + """ + + if rational: + check_matrix(poles, [None, None, 4], scalar_types ) + else: + check_matrix(poles, [None, None, 3], scalar_types) + + assert len(poles) > 0 + assert len(poles[0]) > 0 + self.Nu = len(poles) + self.Nv = len(poles[0]) + for row in poles: + assert len(row) == self.Nv + + assert (not rational and len(poles[0][0]) == 3) or (rational and len(poles[0][0]) == 4) + + (u_knots, v_knots) = knots + check_knots(degree[0], u_knots, self.Nu) + check_knots(degree[1], v_knots, self.Nv) + + self.poles=poles + self.knots=knots + self.rational=rational + self.degree=degree + + def _eval_check(self, u, v, point): + if hasattr(self, '_bs_surface'): + repr_pt = self._bs_surface.eval(u, v) + if not np.allclose(np.array(point), repr_pt, rtol = 1.0e-3): + raise Exception("Point: {} far from curve repr: {}".format(point, repr_pt)) + + + def _dfs(self, groups): + if not hasattr(self, 'id'): + id = len(groups['surfaces']) + 1 + self.id = id + groups['surfaces'].append(self) + + def _brep_output(self, stream, groups): + #writes b-spline surface + stream.write("9 {} {} 0 0 ".format(int(self.rational),int(self.rational))) #prints B-spline surface u or v rational flag - both same + for i in self.degree: #prints <_> + stream.write(" {}".format(i)) + (u_knots, v_knots) = self.knots + stream.write(" {} {} {} {} ".format(self.Nu, self.Nv, len(u_knots), len(v_knots))) + #prints <_> <_> <_> +# stream.write(" {}".format(self.poles)) #TODO: tohle smaz, koukam na format poles a chci: B-spline surface weight poles + for pole in self.poles: #TODO: check, takovy pokus o poles + for vector in pole: + for value in vector: + stream.write(" {}".format(value)) + stream.write(" ") + stream.write(" ") + for knot in u_knots: #prints B-spline surface u multiplicity knots + for value in knot: + stream.write(" {}".format(value)) + stream.write(" ") + for knot in v_knots: #prints B-spline surface v multiplicity knots + for value in knot: + stream.write(" {}".format(value)) + stream.write(" ") + stream.write("\n") + +class Approx: + """ + Approximation methods for B/splines of degree 2. + + """ + @classmethod + def plane(cls, vtxs): + """ + Returns B-spline surface of a plane given by 3 points. + We retun also list of UV coordinates of the given points. + :param vtxs: List of tuples (X,Y,Z) + :return: ( Surface, vtxs_uv ) + """ + assert len(vtxs) == 3, "n vtx: {}".format(len(vtxs)) + vtxs.append( (0,0,0) ) + vtxs = np.array(vtxs) + vv = vtxs[1] + vtxs[2] - vtxs[0] + vtx4 = [ vtxs[0], vtxs[1], vv, vtxs[2]] + (surf, vtxs_uv) = cls.bilinear(vtx4) + return (surf, [ vtxs_uv[0], vtxs_uv[1], vtxs_uv[3] ]) + + @classmethod + def bilinear(cls, vtxs): + """ + Returns B-spline surface of a bilinear surface given by 4 corner points. + We retun also list of UV coordinates of the given points. + :param vtxs: List of tuples (X,Y,Z) + :return: ( Surface, vtxs_uv ) + """ + assert len(vtxs) == 4, "n vtx: {}".format(len(vtxs)) + vtxs = np.array(vtxs) + def mid(*idx): + return np.mean( vtxs[list(idx)], axis=0) + + # v - direction v0 -> v2 + # u - direction v0 -> v1 + poles = [ [vtxs[0], mid(0, 3), vtxs[3]], + [mid(0,1), mid(0,1,2,3), mid(2,3)], + [vtxs[1], mid(1,2), vtxs[2]] + ] + knots = [(0.0, 3), (1.0, 3)] + surface = Surface(poles, (knots, knots)) + vtxs_uv = [ (0, 0), (1, 0), (1, 1), (0, 1) ] + return (surface, vtxs_uv) + + + + + @classmethod + def _line(cls, vtxs, overhang=0.0): + ''' + :param vtxs: List of tuples (X,Y) or (X,Y,Z) + :return: + ''' + assert len(vtxs) == 2 + vtxs = np.array(vtxs) + mid = np.mean(vtxs, axis=0) + poles = [ vtxs[0], mid, vtxs[1] ] + knots = [(0.0+overhang, 3), (1.0-overhang, 3)] + return (poles, knots) + + @classmethod + def line_2d(cls, vtxs): + """ + Return B-spline approximation of line from two 2d points + :param vtxs: [ (X0, Y0), (X1, Y1) ] + :return: Curve2D + """ + return Curve2D( *cls._line(vtxs) ) + + @classmethod + def line_3d(cls, vtxs): + """ + Return B-spline approximation of line from two 3d points + :param vtxs: [ (X0, Y0, Z0), (X1, Y1, Z0) ] + :return: Curve2D + """ + return Curve3D(*cls._line(vtxs)) + + +class Orient(enum.IntEnum): + Forward=1 + Reversed=2 + Internal=3 + External=4 + +#op=Orient.Forward +#om=Orient.Reversed +#oi=Orient.Internal +#oe=Orient.External + +class ShapeRef: + """ + Auxiliary data class to store an object with its orientation + and possibly location. Meaning of location in this context is not clear yet. + Identity location (0) in all BREPs produced by OCC. + All methods accept the tuple (shape, orient, location) and + construct the ShapeRef object automatically. + """ + + orient_chars = ['+', '-', 'i', 'e'] + + def __init__(self, shape, orient=Orient.Forward, location=Location()): + """ + :param shape: referenced shape + :param orient: orientation of the shape, value is enum Orient + :param location: A Location object. Default is None = identity location. + """ + if not issubclass(type(shape), Shape): + raise ParamError("Expected Shape, get: {}.".format(shape)) + assert isinstance(orient, Orient) + assert issubclass(type(location), Location) + + self.shape=shape + self.orientation=orient + self.location=location + + def _writeformat(self, stream, groups): + + stream.write("{}{} {} ".format(self.orient_chars[self.orientation-1], self.shape.id, self.location.id)) + + def __repr__(self): + return "{}{} ".format(self.orient_chars[self.orientation-1], self.shape.id) + + +class ShapeFlag(dict): + """ + Auxiliary data class representing the shape flag word of BREP shapes. + All methods set the flags automatically, but it can be overwritten. + + Free - Seems to indicate a top level shapes. + Modified - ?? + Checked - for format version 2 may indicate that shape topology is already checked + Orientable - ?? + Closed - used to indicate closed Wires and Shells + Infinite - ?? may indicate shapes extending to infinite, not our case + Convex - ?? may indicate convexity of the shape, not clear how this is combined with geometry + """ + flag_names = ['free', 'modified', 'checked', 'orientable', 'closed', 'infinite', 'convex'] + + def __init__(self, *args): + for k, f in zip(self.flag_names, args): + assert f in [0, 1] + self[k]=f + + def set(self, key, value=1): + if value: + value =1 + else: + value =0 + self[key] = value + + def unset(self, key): + self[key] = 0 + + def _brep_output(self, stream): + for k in self.flag_names: + stream.write(str(self[k])) + +class Shape: + def __init__(self, childs): + """ + Construct base Shape object. + Examples: + Wire([ edge_1, edge_2.m(), edge_3]) # recommended + Wire(edge_1, ShapeRef(edge_2, Orient.Reversed, some_location), edge_3) + ... not recommended since it is bad idea to reference same shape with different Locations. + + :param childs: List of ShapeRefs or child objects. + """ + + # self.subtypes - List of allowed types of childs. + assert hasattr(self, 'sub_types'), self + + # convert list of shape reference tuples to ShapeRef objects + # automaticaly wrap naked shapes into tuple. + self.childs=[] + for child in childs: + self.append(child) # append convert to ShapeRef + + # Thes flags are usualy produced by OCC for all other shapes safe vertices. + self.flags=ShapeFlag(0,1,0,1,0,0,0) + + # self.shpname: Shape name, defined in childs + assert hasattr(self, 'shpname'), self + + + """ + Methods to simplify ceration of oriented references. + """ + def p(self): + return ShapeRef(self, Orient.Forward) + + def m(self): + return ShapeRef(self, Orient.Reversed) + + def i(self): + return ShapeRef(self, Orient.Internal) + + def e(self): + return ShapeRef(self, Orient.External) + + def subshapes(self): + # Return list of subshapes stored in child ShapeRefs. + return [chld.shape for chld in self.childs] + + def append(self, shape_ref): + """ + Append a reference to shild + :param shape_ref: Either ShapeRef or child shape. + :return: None + """ + if type(shape_ref) != ShapeRef: + shape_ref=ShapeRef(shape_ref) + if not isinstance(shape_ref.shape, tuple(self.sub_types)): + raise ParamError("Wrong child type: {}, allowed: {}".format(type(shape_ref.shape), self.sub_types)) + self.childs.append(shape_ref) + + #def _convert_to_shaperefs(self, childs): + + + def set_flags(self, flags): + """ + Set flags given as tuple. + :param flags: Tuple of 7 flags. + :return: + """ + self.flags = ShapeFlag(*flags) + + + def is_closed(self): + return self.flags['closed'] + + + def _dfs(self, groups): + """ + Deep first search that assign numbers to shapes + :param groups: dict(locations=[], curves_2d=[], curves_3d=[], surfaces=[], shapes=[]) + :return: None + """ + if hasattr(self, 'id'): + return + for sub_ref in self.childs: + sub_ref.location._dfs(groups) + sub_ref.shape._dfs(groups) + groups['shapes'].append(self) + self.id = len(groups['shapes']) + + def _brep_output(self, stream, groups): + stream.write("{}\n".format(self.shpname)) + self._subrecordoutput(stream) + self.flags._brep_output(stream) + stream.write("\n") +# stream.write("{}".format(self.childs)) + for child in self.childs: + child._writeformat(stream, groups) + stream.write("*\n") + #subshape, tj. childs + + def _subrecordoutput(self, stream): + stream.write("\n") + + def __repr__(self): + if not hasattr(self, 'id'): + self.index_all() + if len(self.childs)==0: + return "" + repr = "" + repr+=self.shpname + " " + str(self.id) + " : " + for child in self.childs: + repr+=child.__repr__() + repr+="\n" + for child in self.childs: + repr+=child.shape.__repr__() + repr+="\n" + return repr + + + def index_all(self, location=Location()): + # print("Index") + # print(compound.__class__.__name__) #prints class name + + groups = dict(locations=[], curves_2d=[], curves_3d=[], surfaces=[], shapes=[]) + self._dfs(groups) # pridej jako parametr dictionary listu jednotlivych grup. v listech primo objekty + location._dfs(groups) + # print(groups) + return groups + + +""" +Shapes with no special parameters, only flags and subshapes. +Writer can be generic implemented in bas class Shape. +""" + +class Compound(Shape): + def __init__(self, shapes=[]): + self.sub_types = [CompoundSolid, Solid, Shell, Wire, Face, Edge, Vertex] + self.shpname = 'Co' + super().__init__(shapes) + #flags: free, modified, IGNORED, orientable, closed, infinite, convex + self.set_flags( (1, 1, 0, 0, 0, 0, 0) ) # free, modified + + def set_free_shapes(self): + """ + Set 'free' attributes to all shapes of the compound. + :return: + """ + for shape in self.subshapes(): + shape.flags.set('free', True) + + +class CompoundSolid(Shape): + def __init__(self, solids=[]): + self.sub_types = [Solid] + self.shpname = 'Cs' + super().__init__(solids) + + +class Solid(Shape): + def __init__(self, shells=[]): + self.sub_types = [Shell] + self.shpname='So' + super().__init__(shells) + self.set_flags((0, 1, 0, 0, 0, 0, 0)) # modified + +class Shell(Shape): + def __init__(self, faces=[]): + self.sub_types = [Face] + self.shpname='Sh' + super().__init__(faces) + self.set_flags((0, 1, 0, 1, 0, 0, 0)) # modified, orientable + + +class Wire(Shape): + def __init__(self, edges=[]): + self.sub_types = [Edge] + self.shpname='Wi' + super().__init__(edges) + self.set_flags((0, 1, 0, 1, 0, 0, 0)) # modified, orientable + self._set_closed() + + def _set_closed(self): + ''' + Return true for the even parity of vertices. + :return: REtrun true if wire is closed. + ''' + vtx_set = {} + for edge in self.subshapes(): + for vtx in edge.subshapes(): + vtx_set[vtx] = 0 + vtx.n_edges += 1 + closed = True + for vtx in vtx_set.keys(): + if vtx.n_edges % 2 != 0: + closed = False + vtx.n_edges = 0 + self.flags.set('closed', closed) + + +""" +Shapes with special parameters. +Specific writers are necessary. +""" + +class Face(Shape): + """ + Face class. + Like vertex and edge have some additional parameters in the BREP format. + """ + + def __init__(self, wires, surface=None, location=Location(), tolerance=1.0e-3): + """ + :param wires: List of wires, or list of edges, or list of ShapeRef tuples of Edges to construct a Wire. + :param surface: Representation of the face, surface on which face lies. + :param location: Location of the surface. + :param tolerance: Tolerance of the representation. + """ + self.sub_types = [Wire, Edge] + self.tol=tolerance + self.restriction_flag =0 + self.shpname = 'Fa' + + if type(wires) != list: + wires = [ wires ] + assert(len(wires) > 0) + super().__init__(wires) + + # auto convert list of edges into wire + shape_type = type(self.childs[0].shape) + for shape in self.subshapes(): + assert type(shape) == shape_type + if shape_type == Edge: + wire = Wire(self.childs) + self.childs = [] + self.append(wire) + + # check that wires are closed + for wire in self.subshapes(): + if not wire.is_closed(): + raise Exception("Trying to make face from non-closed wire.") + + if surface is None: + self.repr=[] + else: + assert type(surface) == Surface + self.repr=[(surface, location)] + + + + def _dfs(self, groups): + Shape._dfs(self,groups) + if not self.repr: + self.implicit_surface() + assert len(self.repr) == 1 + for repr, loc in self.repr: + repr._dfs(groups) + loc._dfs(groups) + + # update geometry representation of edges (add 2D curves) + for wire in self.subshapes(): + for edge in wire.subshapes(): + edge._dfs(groups) + + + def implicit_surface(self): + """ + Construct a surface if surface is None. Works only for + 3 and 4 vertices (plane or bilinear surface) + Should be called in _dfs just after all child shapes are passed. + :return: None + """ + edges = {} + vtxs = [] + for wire in self.subshapes(): + for edge in wire.childs: + edges[edge.shape.id] = edge.shape + e_vtxs = edge.shape.subshapes() + if edge.orientation == Orient.Reversed: + e_vtxs.reverse() + for vtx in e_vtxs: + vtxs.append( (vtx.id, vtx.point) ) + vtxs = vtxs[1:] + vtxs[:1] + odd_vtx = vtxs[1::2] + even_vtx = vtxs[0::2] + assert odd_vtx == even_vtx, "odd: {} even: {}".format(odd_vtx, even_vtx) + vtxs = odd_vtx + if len(vtxs) == 3: + constructor = Approx.plane + elif len(vtxs) == 4: + constructor = Approx.bilinear + else: + raise Exception("Too many vertices {} for implicit surface construction.".format(len(vtxs))) + (ids, points) = zip(*vtxs) + (surface, vtxs_uv) = constructor(list(points)) + self.repr = [(surface, Location())] + + # set representation of edges + assert len(ids) == len(vtxs_uv) + id_to_uv = dict(zip(ids, vtxs_uv)) + for edge in edges.values(): + e_vtxs = edge.subshapes() + v0_uv = id_to_uv[e_vtxs[0].id] + v1_uv = id_to_uv[e_vtxs[1].id] + edge.attach_to_plane( surface, v0_uv, v1_uv ) + + # TODO: Possibly more general attachment of edges to 2D curves for general surfaces, but it depends + # on organisation of intersection curves. + + def _subrecordoutput(self, stream): + assert len(self.repr) == 1 + surf,loc = self.repr[0] + stream.write("{} {} {} {}\n\n".format(self.restriction_flag,self.tol,surf.id,loc.id)) + + +class Edge(Shape): + """ + Edge class. Special edge flags have unclear meaning. + Allow setting representations of the edge, this is crucial for good mash generation. + """ + + class Repr(enum.IntEnum): + Curve3d = 1 + Curve2d = 2 + #Continuous2d=3 + + + def __init__(self, vertices, tolerance=1.0e-3): + """ + :param vertices: List of shape reference tuples, see ShapeRef class. + :param tolerance: Tolerance of the representation. + """ + self.sub_types = [Vertex] + self.shpname = 'Ed' + self.tol = tolerance + self.repr = [] + self.edge_flags=(1,1,0) # this is usual value + + assert(len(vertices) == 2) + + super().__init__(vertices) + # Overwrite vertex orientation + self.childs[0].orientation = Orient.Forward + self.childs[1].orientation = Orient.Reversed + + def set_edge_flags(self, same_parameter, same_range, degenerated): + """ + Edge flags with unclear meaning. + :param same_parameter: + :param same_range: + :param degenerated: + :return: + """ + self.edge_flags=(same_parameter,same_range, degenerated) + + def points(self): + ''' + :return: List of coordinates of the edge vertices. + ''' + return [ vtx.point for vtx in self.subshapes()] + + def attach_to_3d_curve(self, t_range, curve, location=Location()): + """ + Add vertex representation on a 3D curve. + :param t_range: Tuple (t_min, t_max). + :param curve: 3D curve object (Curve3d) + :param location: Location object. Default is None = identity location. + :return: None + """ + assert type(curve) == Curve3D + curve._eval_check(t_range[0], self.points()[0]) + curve._eval_check(t_range[1], self.points()[1]) + self.repr.append( (self.Repr.Curve3d, t_range, curve, location) ) + + def attach_to_2d_curve(self, t_range, curve, surface, location=Location()): + """ + Add vertex representation on a 2D curve. + :param t_range: Tuple (t_min, t_max). + :param curve: 2D curve object (Curve2d) + :param surface: Surface on which the curve lies. + :param location: Location object. Default is None = identity location. + :return: None + """ + assert type(surface) == Surface + assert type(curve) == Curve2D + curve._eval_check(t_range[0], surface, self.points()[0]) + curve._eval_check(t_range[1], surface, self.points()[1]) + self.repr.append( (self.Repr.Curve2d, t_range, curve, surface, location) ) + + def attach_to_plane(self, surface, v0, v1): + """ + Construct and attach 2D line in UV space of the 'surface' + :param surface: A Surface object. + :param v0: UV coordinate of the first edge point + :param v1: UV coordinate of the second edge point + :return: + """ + assert type(surface) == Surface + + self.attach_to_2d_curve((0.0, 1.0), Approx.line_2d([v0, v1]), surface) + + def implicit_curve(self): + """ + Construct a line 3d curve if there is no 3D representation. + Should be called in _dfs. + :return: + """ + vtx_points = self.points() + self.attach_to_3d_curve((0.0,1.0), Approx.line_3d( vtx_points )) + + #def attach_continuity(self): + + def _dfs(self, groups): + Shape._dfs(self,groups) + if not self.repr: + self.implicit_curve() + assert len(self.repr) > 0 + for i,repr in enumerate(self.repr): + if repr[0]==self.Repr.Curve2d: + repr[2]._dfs(groups) #curve + repr[3]._dfs(groups) #surface + repr[4]._dfs(groups) #location + elif repr[0]==self.Repr.Curve3d: + repr[2]._dfs(groups) #curve + repr[3]._dfs(groups) #location + + def _subrecordoutput(self, stream): #prints edge data #TODO: tisknu nekolik data representation + assert len(self.repr) > 0 + stream.write(" {} {} {} {}\n".format(self.tol,self.edge_flags[0],self.edge_flags[1],self.edge_flags[2])) + for i,repr in enumerate(self.repr): + if repr[0] == self.Repr.Curve2d: + curve_type, t_range, curve, surface, location = repr + stream.write("2 {} {} {} {} {}\n".format(curve.id, surface.id, location.id,t_range[0],t_range[1] )) #TODO: 2 <_> <_> + elif repr[0] == self.Repr.Curve3d: + curve_type, t_range, curve, location = repr + stream.write("1 {} {} {} {}\n".format(curve.id, location.id, t_range[0], t_range[1])) #TODO: 3 + stream.write("0\n") + +class Vertex(Shape): + """ + Vertex class. + Allow setting representations of the vertex but seems it is not used in BREPs produced by OCC. + """ + + class Repr(enum.IntEnum): + Curve3d = 1 + Curve2d = 2 + Surface = 3 + + def __init__(self, point, tolerance=1.0e-3): + """ + :param point: 3d point (X,Y,Z) + :param tolerance: Tolerance of the representation. + """ + check_matrix(point, [3], scalar_types) + + # These flags are produced by OCC for vertices. + self.flags = ShapeFlag(0, 1, 0, 1, 1, 0, 1) + # Coordinates in the 3D space. [X, Y, Z] + self.point=np.array(point) + # tolerance of representations. + self.tolerance=tolerance + # List of geometrical representations of the vertex. Possibly not necessary for meshing. + self.repr=[] + # Number of edges in which vertex is used. Used internally to check closed wires. + self.n_edges = 0 + self.shpname = 'Ve' + self.sub_types=[] + + super().__init__(childs=[]) + + def attach_to_3d_curve(self, t, curve, location=Location()): + """ + Add vertex representation on a 3D curve. + :param t: Parameter of the point on the curve. + :param curve: 3D curve object (Curve3d) + :param location: Location object. Default is None = identity location. + :return: None + """ + curve._eval_check(t, self.point) + self.repr.append( (self.Repr.Curve3d, t, curve, location) ) + + def attach_to_2d_curve(self, t, curve, surface, location=Location()): + """ + Add vertex representation on a 2D curve on a surface. + :param t: Parameter of the point on the curve. + :param curve: 2D curve object (Curve2d) + :param surface: Surface on which the curve lies. + :param location: Location object. Default is None = identity location. + :return: None + """ + curve._eval_check(t, surface, self.point) + self.repr.append( (self.Repr.Curve2d, t, curve, surface, location) ) + + def attach_to_surface(self, u, v, surface, location=Location()): + """ + Add vertex representation on a 3D curve. + :param u,v: Parameters u,v of the point on the surface. + :param surface: Surface object. + :param location: Location object. Default is None = identity location. + :return: None + """ + surface._eval_check(u, v, self.point) + self.repr.append( (self.Repr.Surface, u,v, surface, location) ) + + + def _dfs(self, groups): + Shape._dfs(self,groups) + for repr in self.repr: + if repr[0]==self.Repr.Curve2d: + repr[2]._dfs(groups) #curve + repr[3]._dfs(groups) #surface + repr[4]._dfs(groups) #location + elif repr[0]==self.Repr.Curve3d: + repr[2]._dfs(groups) #curve + repr[3]._dfs(groups) #location + + + def _subrecordoutput(self, stream): #prints vertex data + stream.write("{}\n".format(self.tolerance)) + for i in self.point: + stream.write("{} ".format(i)) + stream.write("\n0 0\n\n") #no added + + + + +def write_model(stream, compound, location): + + groups = compound.index_all(location=location) + + # fix shape IDs + n_shapes = len(groups['shapes']) + 1 + for shape in groups['shapes']: + shape.id = n_shapes - shape.id + + stream.write("DBRep_DrawableShape\n\n") + stream.write("CASCADE Topology V1, (c) Matra-Datavision\n") + stream.write("Locations {}\n".format(len(groups['locations']))) + for loc in groups['locations']: + loc._brep_output(stream, groups) + + stream.write("Curve2ds {}\n".format(len(groups['curves_2d']))) + for curve in groups['curves_2d']: + curve._brep_output(stream, groups) + + stream.write("Curves {}\n".format(len(groups['curves_3d']))) + for curve in groups['curves_3d']: + curve._brep_output(stream, groups) + + stream.write("Polygon3D 0\n") + + stream.write("PolygonOnTriangulations 0\n") + + stream.write("Surfaces {}\n".format(len(groups['surfaces']))) + for surface in groups['surfaces']: + surface._brep_output(stream, groups) + + stream.write("Triangulations 0\n") + + stream.write("\nTShapes {}\n".format(len(groups['shapes']))) + for shape in groups['shapes']: + #print("# {} id: {} childs: {}\n".format(shape.shpname, shape.id, + # [ch.id for ch in shape.subshapes()])) + shape._brep_output(stream, groups) + stream.write("\n+1 0") + #stream.write("0\n") + + diff --git a/src/bspline.py b/src/bspline.py new file mode 100644 index 0000000..7687026 --- /dev/null +++ b/src/bspline.py @@ -0,0 +1,1095 @@ +""" +Module with classes representing various B-spline and NURBS curves and surfaces. +These classes provide just basic functionality: +- storing the data +- evaluation of XYZ for UV +- evaluation and xy<->uv functions accepting np.arrays, +- evaluation of derivatives +In future: +- use de Boor algorithm for evaluation of curves and surfaces +- serialization and deserialization using JSONdata - must make it an installable module +- implement degree increasing and knot insertion +""" + + +import numpy as np +import numpy.linalg as la +import copy + +__author__ = 'Jan Brezina , Jiri Hnidek , Jiri Kopal ' + + +class SplineBasis: + """ + Represents a spline basis for a given knot vector and degree. + Provides canonical evaluation for the bases functions and their derivatives, knot vector lookup etc. + """ + + @classmethod + def make_equidistant(cls, degree, n_intervals, knot_range=[0.0, 1.0]): + """ + Returns spline basis for an eqidistant knot vector + having 'n_intervals' subintervals. + :param degree: degree of the spline basis + :param n_intervals: Number of subintervals. + :param knot_range: support of the spline, min and max valid 't' + :return: SplineBasis object. + """ + n = n_intervals + 2 * degree + 1 + knots = np.array((knot_range[0],) * n) + diff = (knot_range[1] - knot_range[0]) / n_intervals + for i in range(degree + 1, n - degree): + knots[i] = (i - degree) * diff + knot_range[0] + knots[-degree - 1:] = knot_range[1] + return cls(degree, knots) + + + @classmethod + def make_from_packed_knots(cls, degree, knots): + """ + Construct basis from the vector of packed knots. + :param degree: Degree of the basis. + :param knots: List of knots with their multiplicities, [ (knot, mult), ..] + :return: SplineBasis object. + """ + full_knots = [ q for q, mult in knots for i in range(mult) ] + return cls(degree, full_knots) + + + def __init__(self, degree, knots): + """ + Constructor of the spline basis. + :param degree: Degree of Bezier polynomials >=0. + :param knots: Numpy array of the knots including multiplicities. + """ + assert degree >=0 + self.degree = degree + + # check free ends + for i in range(self.degree): + assert knots[i] == knots[i+1] + assert knots[-i-1] == knots[-i-2] + self.knots = np.array(knots) + + self.size = len(self.knots) - self.degree -1 + # Number of basis functions. + + + self.knots_idx_range = [self.degree, len(self.knots) - self.degree - 1] + # Range of knot indices corrsponding to the basis domain. + + self.domain = self.knots[self.knots_idx_range] + # Support domain of the spline. + + self.domain_size = self.domain[1] - self.domain[0] + # Size of the domain. + + self.n_intervals = self.size - self.degree + # Number of subintervals ( assuming multiplicities only on ends. ) + + # Set optimized functions for specific degrees. + if self.degree == 2: + self.eval_base_vector = self._eval_vector_deg_2 + self.eval_diff_base_vector = self._eval_diff_vector_deg_2 + + + def pack_knots(self): + """ + :return: Packed knot vector, [ (knot, multiplicity), .. ] + """ + last, mult = self.knots[0], 0 + packed_knots = [] + for q in self.knots: + if q == last: + mult+=1 + else: + packed_knots.append( (last, mult) ) + last, mult = q, 1 + packed_knots.append((last,mult)) + return packed_knots + + + def check(self, t, rtol = 1e-10): + """ + Check that 't' is inside basis domain. Fix small perturbations. + :param t: + :return: + """ + if self.domain[0] <= t <= self.domain[1]: + return t + else: + tol = (np.abs(t) + 1e-4)*rtol + if not self.domain[0] - tol <= t <= self.domain[1] + tol: + raise IndexError("Evaluate spline, t={}, out of domain: {}.".format(t, self.domain)) + + if t < self.domain[0]: + return self.domain[0] + else: + return self.domain[1] + + self.domain[0] + + def find_knot_interval(self, t): + """ + Find the first non-empty knot interval containing the value 't'. + i.e. knots[i] <= t < knots[i+1], where knots[i] < knots[i+1] + Returns I = i - degree, which is the index of the first basis function + nonzero in 't'. + + :param t: float, must be within knots limits. + :return: I + """ + idx = np.searchsorted(self.knots[self.degree: -self.degree -1], [t], side='right')[0] - 1 + assert 0 <= idx <= self.n_intervals, "Evaluation out of spline domain; t: {} min: {} max: {}".format(t, self.knots[0], self.knots[-1]) + return min(idx, self.n_intervals - 1) # deals with t == self.knots[-1] + + + + def _basis(self, deg, idx, t): + """ + Recursive evaluation of basis function of given degree and index. + + :param deg: Degree of the basis function + :param idx: Index of the basis function to evaluate. + :param t: Point of evaluation. + :return Value of the basis function. + """ + + if deg == 0: + t_0 = self.knots[idx] + t_1 = self.knots[idx + 1] + value = 1.0 if t_0 <= t < t_1 else 0.0 + return value + else: + t_i = self.knots[idx] + t_ik = self.knots[idx + deg] + top = t - t_i + bottom = t_ik - t_i + if bottom != 0: + value = top / bottom * self._basis(deg-1, idx, t) + else: + value = 0.0 + + t_ik1 = self.knots[idx + deg + 1] + t_i1 = self.knots[idx + 1] + top = t_ik1 - t + bottom = t_ik1 - t_i1 + if bottom != 0: + value += top / bottom * self._basis(deg-1, idx+1, t) + + return value + + def fn_supp(self, i_base): + """ + Support of the base function 'i_base'. + :param i_base: + :return: (t_min, t_max) + """ + return (self.knots[i_base], self.knots[i_base + self.degree + 1]) + + + def eval(self, i_base, t): + """ + :param i_base: Index of base function to evaluate. + :param t: point in which evaluate + :return: b_i(y) + """ + assert 0 <= i_base < self.size + if i_base == self.size -1 and t == self.domain[1]: + return 1.0 + return self._basis(self.degree, i_base, t) + + + def eval_diff(self, i_base, t): + assert 0 <= i_base < self.size + deg = self.degree + i = i_base + if i_base == self.size - 2 and t == self.domain[1]: + return - deg / (self.knots[i + deg + 1] - self.knots[i + 1]) + + if i_base == self.size - 1 and t == self.domain[1]: + return deg / (self.knots[i + deg] - self.knots[i]) + + diff = 0.0 + if i > 0: + B = self._basis(deg-1, i, t) + diff += deg * B / (self.knots[i + deg] - self.knots[i]) + if i < self.size - 1: + B = self._basis(deg-1, i + 1, t) + diff -= deg * B / (self.knots[i + deg + 1] - self.knots[i + 1]) + + return diff + + def make_linear_poles(self): + """ + Return poles of basis functions to get a f(x) = x. + :return: + """ + poles= [ 0.0 ] + for i in range(self.size-1): + pp = poles[-1] + (self.knots[i + self.degree + 1] - self.knots[i + 1]) / float(self.degree) + poles.append(pp) + return poles + + + def eval_vector(self, i_base, t): + """ + This function compute base function of B-Spline curve on given subinterval. + :param i_int: Interval in which 't' belongs. Three nonzero basis functions on this interval are evaluated. + :param t: Where to evaluate. + :return: Numpy array of three values. + """ + values = [] + for ib in range(i_base, i_base + self.degree + 1): + values.append( self.eval(ib, t)) + return values + + + def eval_diff_vector(self, i_base, t): + """ + This function compute derivative of base function of B-Spline curve on given subinterval. + :param i_int: Interval in which 't' belongs. Derivatives of the 3 nonzero basis functions on this interval are evaluated. + :param t: Where to evaluate. + :return: Numpy array of three values. + """ + values = [] + for ib in range(i_base, i_base + self.degree + 1): + values.append( self.eval_diff(ib, t)) + return values + + + """ + Specializations. + TODO: + - Try usage of scipy evaluation, compare speed with optimized eval_vector and diff_eval_vector. + - Generalize optimized evaluation of eval_vector (If scipy is not faster), De Boor algortihm, Hoschek 4.3.3. + - Optimize eval and eval_diff - without recursion, based on combinatorGeneralize optimized evaluation of eval_vector (If scipy is not faster) + """ + def _eval_vector_deg_2(self, i_int, t): + """ + This function compute base function of B-Spline curve on given subinterval. + :param i_int: Interval in which 't' belongs. Three nonzero basis functions on this interval are evaluated. + :param t: Where to evaluate. + :return: Numpy array of three values. + Note: Keep code redundancy with 'diff' as optimization. + """ + + basis_values = np.zeros(3) + + tk1, tk2, tk3, tk4 = self.knots[i_int + 1 : i_int + 5] + + d31 = tk3 - tk1 + d32 = tk3 - tk2 + d42 = tk4 - tk2 + + dt1 = t - tk1 + dt2 = t - tk2 + d3t = tk3 - t + d4t = tk4 - t + + d31_d32 = d31 * d32 + d42_d32 = d42 * d32 + + basis_values[0] = (d3t * d3t) / d31_d32 + basis_values[1] = ((dt1 * d3t) / d31_d32) + ((dt2 * d4t) / d42_d32) + basis_values[2] = (dt2 * dt2) / d42_d32 + + return basis_values + + + def _eval_diff_vector_deg_2(self, i_int, t): + """ + This function compute derivative of base function of B-Spline curve on given subinterval. + :param i_int: Interval in which 't' belongs. Derivatives of the 3 nonzero basis functions on this interval are evaluated. + :param t: Where to evaluate. + :return: Numpy array of three values. + Note: Keep code redundancy with 'diff' as optimization. + """ + + basis_values = np.zeros(3) + + tk1, tk2, tk3, tk4 = self.knots[i_int + 1: i_int + 5] + + d31 = tk3 - tk1 + d32 = tk3 - tk2 + d42 = tk4 - tk2 + + dt1 = t - tk1 + dt2 = t - tk2 + d3t = tk3 - t + d4t = tk4 - t + + d31_d32 = d31 * d32 + d42_d32 = d42 * d32 + + basis_values[0] = -2 * d3t / d31_d32 + basis_values[1] = (d3t - dt1) / d31_d32 + (d4t - dt2) / d42_d32 + basis_values[2] = 2 * dt2 / d42_d32 + + return basis_values + + +class Curve: + """ + Defines a D-dim B-spline curve. + """ + + @classmethod + def make_raw(cls, poles, knots, rational=False, degree=2): + """ + Construct a B-spline curve. + :param poles: Numpy array N x (D+r) of poles (control points). N is number of poles, D is dimension of the curve, 'r' is 1 for rational curves. + For rational case, poles[:, D] are weights of the control points. + :param knots, degree: See SplineBasis. + :param rational: True for rational B-spline, i.e. NURB. Use weighted poles. + :param degree: Non-negative int + """ + basis = SplineBasis(degree, knots) + return cls(basis, poles, rational) + + def __init__(self, basis, poles, rational = False): + """ + Construct a B-spline curve. + :param poles: Numpy array N x (D+r) of poles (control points). N is number of poles, D is dimension of the curve, 'r' is 1 for rational curves. + For rational case, poles[:, D] are weights of the control points. + :param basis: SplineBasis object. + :param rational: True for rational B-spline, i.e. NURB. Use weighted poles. + """ + + self.basis = basis + # Spline basis. + + self.poles = np.array(poles, dtype=float) # N x D + assert self.poles.shape[0] == self.basis.size + # Spline poles. + + self.dim = len(poles[0]) - rational + # Dimension of the curve. + + self.rational = rational + # Indicator of rational B-spline (NURBS). + + if rational: + # precomputations + self._weights = poles[:, self.dim] + self._poles = (poles[:, 0:self.dim].T * self._weights ).T + + + def eval_local(self, t, it): + """ + Evaluate a B-spline curve for parameter 't' with given knotinterval 'it'. + :param t: Evaluation point. + :param it: Index of knot subinterval (see doc of 'find_knot_interval') + :return: D-dimensional numpy array. D - is dimension given by dimension of poles. + """ + dt = self.basis.degree + 1 + t_base_vec = self.basis.eval_vector(it, t) + + if self.rational: + top_value = np.inner(t_base_vec, self._poles[it: it + dt, :].T) + bot_value = np.inner(t_base_vec, self._weights[it: it + dt]) + return top_value / bot_value + else: + return np.inner(t_base_vec, self.poles[it: it + dt, :].T) + + + def eval(self, t): + """ + Evaluate a B-spline curve for paramater 't'. Check and fix range of 't'. + :param t: Evaluation point. + :return: D-dimensional numpy array. D - is dimension given by dimension of poles. + TODO: + - test evaluation for rational curves + """ + t = self.basis.check(t) + it = self.basis.find_knot_interval(t) + return self.eval_local(t, it) + + + def eval_array(self, t_points): + """ + Evaluate in array of t-points. + :param t_points: array N x float + :return: Numpy array N x D, D is dimension of the curve. + """ + return np.array( [ self.eval(t) for t in t_points] ) + + + def aabb(self): + """ + Return Axes Aligned Bounding Box of the poles, which should be also bounding box of the curve itself. + :return: ( min_corner, max_corner); Box corners are numpy arryas of dimension D. + """ + return np.array( [np.amin(self.poles, axis=0), np.amax(self.poles, axis=0)] ) + + +class Surface: + """ + Defines D-dim B-spline surface. + """ + + @classmethod + def make_raw(cls, poles, knots, rational=False, degree=(2,2)): + """ + Construct a B-spline curve. + :param poles: List of poles (control points) ( X, Y, Z ) or weighted points (X,Y,Z, w). X,Y,Z,w are floats. + Weighted points are used only for rational B-splines (i.e. nurbs) + :param knots: List of tuples (knot, multiplicity), where knot is float, t-parameter on the curve of the knot + and multiplicity is positive int. Total number of knots, i.e. sum of their multiplicities, must be + degree + N + 1, where N is number of poles. + :param rational: True for rational B-spline, i.e. NURB. Use weighted poles. + :param degree: Non-negative int + """ + u_basis = SplineBasis(degree[0], knots[0]) + v_basis = SplineBasis(degree[1], knots[1]) + return cls( (u_basis, v_basis), poles, rational) + + + def __init__(self, basis, poles, rational=False): + """ + Construct a B-spline surface. + :param poles: Numpy array Nu x Nv x (D+r) of poles (control points). + Nu and Nv are sizes of u_basis, v_basis respectively. + D is dimension of the surface, 'r' is 1 for rational surfaces. + For rational case, poles[:, :, D] are weights of the control points. + :param basis: (u_basis, v_basis) SplineBasis objects for U and V parameter axis. + :param rational: True for rational B-spline, i.e. NURB. Use weighted poles. + """ + self.u_basis, self.v_basis = basis + # Surface basis for U and V axis. + + self.poles = np.array(poles, dtype=float) + self.dim = len(self.poles[0,0,:]) - rational + # Surface dimension, D. + + + # Surface poles matrix: Nu x Nv x (D+r) + assert self.poles.shape == (self.u_basis.size, self.v_basis.size, self.dim + rational) + + self.rational = rational + # Rational surface indicator. + if rational: + # precomputations + self._weights = poles[:, :, self.dim] + self._poles = (poles[:, :, 0:self.dim].T * self._weights.T ).T + + + + + def eval_local(self, u, v, iu, iv): + """ + Evaluate a B-spline surface for paramaters u,v with given knot subintervals. + :param u, v: Evaluation point. + :param iu, iv: Knot subintervals of u, v. + :return: D-dimensional numpy array. D - is dimension given by dimension of poles. + """ + du = self.u_basis.degree + 1 + dv = self.v_basis.degree + 1 + u_base_vec = self.u_basis.eval_vector(iu, u) + v_base_vec = self.v_basis.eval_vector(iv, v) + + if self.rational: + top_value = np.inner(u_base_vec, self._poles[iu: iu + du, iv: iv + dv, :].T) + top_value = np.inner(top_value, v_base_vec) + bot_value = np.inner(u_base_vec, self._weights[iu: iu + du, iv: iv + dv].T) + bot_value = np.inner(bot_value, v_base_vec) + return top_value / bot_value + else: + #print("u: {} v: {} p: {}".format(u_base_vec.shape, v_base_vec.shape, self.poles[iu: iu + du, iv: iv + dv, :].shape)) + # inner product sums along the last axis of its parameters + value = np.inner(u_base_vec, self.poles[iu: iu + du, iv: iv + dv, :].T ) + #print("val: {}".format(value.shape)) + return np.inner( value, v_base_vec) + + + + def eval(self, u, v): + """ + Evaluate a B-spline surface for paramaters u,v. Check and fix range of 'u, v'. + :param u, v: Evaluation point. + :return: D-dimensional numpy array. D - is dimension given by dimension of poles. + TODO: + - test evaluation for rational curves + """ + u = self.u_basis.check(u) + v = self.v_basis.check(v) + iu = self.u_basis.find_knot_interval(u) + iv = self.v_basis.find_knot_interval(v) + return self.eval_local(u, v, iu, iv) + + + + + + + def deep_copy(self): + u_basis = copy.copy(self.u_basis) + v_basis = copy.copy(self.v_basis) + poles = copy.copy(self.poles) + return Surface( ( u_basis, v_basis), poles) + + def eval_array(self, uv_points): + """ + Evaluate in array of uv-points. + :param uv_points: numpy array N x [u, v] + :return: Numpy array N x D; D is dimension of the curve. + """ + assert uv_points.shape[1] == 2 + return np.array( [ self.eval(u, v) for u, v in uv_points] ) + + def aabb(self): + """ + Return Axes Aligned Bounding Box of the poles, which should be also bounding box of the curve itself. + :return: ( min_corner, max_corner); Box corners are numpy arryas of dimension D. + TODO: test + """ + return np.array( [np.amin(self.poles, axis=(0,1)), np.amax(self.poles, axis=(0,1))] ) + + +class Z_Surface: + """ + Simplified B-spline surface that use just linear or bilinear transform between XY and UV. + """ + def __init__(self, xy_quad, z_surface): + """ + Construct a surface given by the 1d surface for the Z coordinate and XY quadrilateral + for the bilinear UV <-> XY mapping. + :param xy_quad: np array N x 2 + Four or three points, determining bilinear or linear mapping, respectively. + Four points giving XY coordinates for the uv corners: (0,1), (0,0), (1,0), (1,1) + Three points giving XY coordinates for the uv corners: (0,1), (0,0), (1,0) + Linear case is also detected for the four points. + :param z_surface: 1D Surface object. + """ + assert z_surface.dim == 1 + self.dim = 3 + # Fixed surface dimension. + + self.z_surface = z_surface + # Underlaying 1d surface object for Z coord evaluation. + + self.u_basis = z_surface.u_basis + self.v_basis = z_surface.v_basis + # Basis for UV directions. + + self.orig_quad = xy_quad + self.quad = None + # Boundary quadrilateral. + + self.reset_transform(xy_quad) + # Set further private attributes, see comment there: + # _z_mat, _have_z_mat, _xy_shift, _mat_xy_to_uv, _mat_uv_to_xy + + def make_full_surface(self): + """ + Return representation of the surface by the 3d Surface object. + Compute redundant XY poles. + :return: Surface. + """ + basis = (self.z_surface.u_basis, self.z_surface.v_basis) + + u = basis[0].make_linear_poles() + v = basis[1].make_linear_poles() + V, U = np.meshgrid(v,u) + uv_poles_vec = np.stack([U.ravel(), V.ravel()], axis=1) + xy_poles = self.uv_to_xy(uv_poles_vec).reshape(U.shape[0], U.shape[1], 2) + z_poles = self.z_surface.poles.copy() + if self._have_z_mat: + z_poles *= self.z_mat[0] + z_poles += self.z_mat[1] + poles = np.concatenate( (xy_poles, z_poles), axis = 2 ) + + return Surface(basis, poles) + + def reset_transform(self, xy_quad=None): + """ + Set XY transform according to given domain quadrilateral (or triangle for linear mapping case). + :param xy_quad: np array N x 2 + Four or three points, determining bilinear or linear mapping, respectively. + Four points giving XY coordinates for the uv corners: (0,1), (0,0), (1,0), (1,1) + Three points giving XY coordinates for the uv corners: (0,1), (0,0), (1,0) + Linear case is also detected for the four points. + + If no xy_quad is given, we reset to the quad used in constructor. + :return: None + """ + # Build envelope quadrilateral polygon in XY plane. + if xy_quad is None: + xy_quad = self.orig_quad + self.quad = np.array(xy_quad, dtype=float) + assert self.quad.shape[0] in [3, 4], "Three or four points must be given." + assert self.quad.shape[1] == 2 + + v11 = self.quad[0] + self.quad[2] - self.quad[1] + if self.quad.shape[0] == 3: + self.quad = np.concatenate( (self.quad, v11[None, :]), axis = 0) + + if np.allclose(self.quad[3], v11): + # linear case + self._xy_shift = self.quad[1] + v_vec = self.quad[0] - self.quad[1] + u_vec = self.quad[2] - self.quad[1] + self._mat_uv_to_xy = np.column_stack((u_vec, v_vec)) + self._mat_xy_to_uv = la.inv(self._mat_uv_to_xy) + + self.xy_to_uv = self._linear_xy_to_uv + self.uv_to_xy = self._linear_uv_to_xy + + else: + # bilinear case + self.xy_to_uv = self._bilinear_xy_to_uv + self.uv_to_xy = self._bilinear_uv_to_xy + + self.z_mat = np.array( [1.0, 0.0] ) + # [ z_scale, z_shift ] + + self._have_z_mat = False + # Indicate that we have z_mat different from identity. + # Optimization for array evaluation methods. z_mat must be set anyway. + + def transform(self, xy_mat, z_mat = None ): + """ + Transform the Z-surface by arbitrary XY linear transform and Z linear transform. + :param xy_mat: np array, 2 rows 3 cols, last column is xy shift + :param z_shift: np.array, [ z_scale, z_shift] + :return: None + """ + if xy_mat is not None: + xy_mat = np.array(xy_mat) + assert xy_mat.shape == (2, 3) + self._mat_uv_to_xy = xy_mat[0:2,0:2].dot( self._mat_uv_to_xy ) + self._xy_shift = xy_mat[0:2,0:2].dot( self._xy_shift ) + xy_mat[0:2, 2] + self._mat_xy_to_uv = la.inv(self._mat_uv_to_xy) + + # transform quad + self.quad = np.dot(self.quad, xy_mat[0:2,0:2].T) + xy_mat[0:2, 2] + + # apply z-transfrom directly to the poles + if z_mat is not None: + z_mat = np.array(z_mat) + assert z_mat.shape == (2,) + self.z_mat[0] *= z_mat[0] + self.z_mat[1] = z_mat[0] * self.z_mat[1] + z_mat[1] + self._have_z_mat = True + + def get_xy_matrix(self): + """ + Return xy_mat of curent XY tranasform. + :return: + """ + return np.concatenate((self._mat_uv_to_xy, self._xy_shift[:, None]), axis=1) + + def apply_z_transform(self): + """ + Make copy of the bs.Surface for Z coordinates and apply current Z transform to the poles of the copy. + Reset the Z transform. + :return: + """ + if self._have_z_mat: + self.z_surface = self.z_surface.deep_copy() + self.z_surface.poles *= self.z_mat[0] + self.z_surface.poles += self.z_mat[1] + self.z_mat = np.array([1.0, 0.0]) + self._have_z_mat = False + + + def get_copy(self): + """ + Returns a copy of the Z_surface with own Z and XY transforms, but shared + bs.Surface for Z coordinates (Z poles). + :return: A copy of the Z_Surface object. + """ + return copy.copy(self) + + + """ + def uv_to_xy(self, uv_points): + + Abstract method. Converts array of uv points to array of xy points. + :param uv_points: numpy array N x [u,v] + :return: numpy array N x [x,y] + """ + def _linear_uv_to_xy(self, uv_points): + assert uv_points.shape[1] == 2, "Size: {}".format(uv_points.shape) + return ( np.dot(uv_points, self._mat_uv_to_xy.T) + self._xy_shift) + + + def _bilinear_uv_to_xy(self, uv_points): + assert uv_points.shape[1] == 2, "Size: {}".format(uv_points.shape) + xy_list = [] + for u,v in uv_points: + weights = np.array([ (1-u)*v, (1-u)*(1-v), u*(1-v), u*v ]) + xy_list.append( self.quad.T.dot(weights) ) + return np.array(xy_list) + + """ + def xy_to_uv(self, xy_points): + + Abstract method. Converts array of xy points to array of uv points. + :param xy_points: numpy array N x [x,y] + :return: numpy array N x [u,v] + """ + def _linear_xy_to_uv(self, xy_points): + # assert xy_points.shape[0] == 2 + assert xy_points.shape[1] == 2 + return np.dot((xy_points - self._xy_shift), self._mat_xy_to_uv.T) + + + def _bilinear_xy_to_uv(self, xy_points): + assert xy_points.shape[1] == 2 + assert False, "Not implemented yet." + + + def eval(self, u, v): + """ + Evaluate a B-spline surface for paramaters u,v. + :param u, v: Evaluation point. + :return: D-dimensional numpy array. D - is dimension given by dimension of poles. + """ + z = self.z_surface.eval(u, v) + z = self.z_mat[0] * z + self.z_mat[1] + uv_points = np.array([[u, v]]) + x, y = self.uv_to_xy( uv_points )[0] + return np.array( [x, y, z] ) + + + def eval_array(self, uv_points): + """ + Evaluate a B-spline surface in array of UV points. + :param uv_points: numpy array N x [u, v] + :return: array N x D; D - is dimension given by dimension of poles. + """ + assert uv_points.shape[1] == 2 + z_points = self.z_surface.eval_array(uv_points) + if self._have_z_mat: + z_points *= self.z_mat[0] + z_points += self.z_mat[1] + + xy_points = self.uv_to_xy(uv_points) + return np.concatenate( (xy_points, z_points), axis = 1) + + + def eval_xy_array(self, xy_points): + """ + Evaluate a B-spline surface in array of XY points. + :param xy_points: numpy array N x [x, y] + :return: array N x D; D - is dimension given by dimension of poles. + """ + uv_points = self.xy_to_uv(xy_points) + return self.eval_array(uv_points) + + + def z_eval_array(self, uv_points): + """ + Evaluate just Z coordinate for array of UV points. + :param uv_points: numpy array N x [u, v] + :return: array N x D; D - is dimension given by dimension of poles. + """ + assert uv_points.shape[1] == 2 + z_points = self.z_surface.eval_array(uv_points) + if self._have_z_mat: + z_points *= self.z_mat[0] + z_points += self.z_mat[1] + return z_points.reshape(-1) + + + def z_eval_xy_array(self, xy_points): + """ + Evaluate just Z coordinate for array of XY points. + :param uv_points: numpy array N x [x, y] + :return: array N x D; D - is dimension given by dimension of poles. + """ + uv_points = self.xy_to_uv(xy_points) + return self.z_eval_array(uv_points) + + + def center(self): + """ + Evaluate surface in center of UV square. + :return: (X, Y, Z) + """ + return self.eval_array( np.array([ [0.5, 0.5] ]))[0] + + + def aabb(self): + xyz_box = np.empty( (2, 3) ) + xyz_box[0, 0:2] = np.amin(self.quad, axis=0) + xyz_box[1, 0:2] = np.amax(self.quad, axis=0) + xyz_box[:, 2] = self.z_mat[0]*self.z_surface.aabb()[:,0] + self.z_mat[1] + return xyz_box + + + +class GridNotInShapeExc(Exception): + pass + +class IrregularGridExc(Exception): + pass + + +class GridSurface: + step_tolerance = 1e-5 + # relative step_tolerance of + + + @staticmethod + def load(filename): + """ + Load the grid surface from file + :param filename: + :return: GridSurface object. + """ + point_seq = np.loadtxt(filename) + assert min(point_seq.shape) > 1 + return GridSurface(point_seq) + + """ + Can load and check grid of XYZ points and construct a + Z-surface of degree 1 for them. + """ + # TODO: calling transform lose mapping between original point grid and unit square, approximation can not be performed + + + def __init__(self, point_seq, tolerance = 1e-6): + """ + Construct the GridSurface from sequence of points. + :param point_seq: N x 3 numpy array; organized as Nu x Nv grid. + Nu, Nv are detected automaticaly, both must be greater then 1. + :param step_tolerance: Tolerance between XY position given by envelope quad and actual point position. + Relative to the length of maximal side of the envelope quad. + """ + self.tolerance = tolerance + + assert point_seq.shape[1] == 3 + n_points = point_seq.shape[0] + point_seq_xy = point_seq[:, 0:2] + + self.quad = None + # Envelope quad - polygon, oriented counter clockwise. + + self.shape = None + # Number of points along axis: Nu x Nv + + self._get_grid_corners(point_seq_xy) + self._check_grid_regularity(point_seq_xy) + + self.z_scale = 1.0 + self.z_shift = 0.0 + + self._uv_step = (1.0 / float(self.shape[0] - 1), 1.0 / float(self.shape[1] - 1)) + # Grid step in u, v direction respectively. + + self.grid_uvz = None + # numpy array nu x nv x 3 with original XYZ values transformed to unit square. + + self._make_z_surface( point_seq) + self._check_map() + + + def _get_grid_corners(self, xy_points): + n_points = len(xy_points) + vtx_00 = xy_points[0, 0:2] + vtx_dv = xy_points[1, 0:2] + + # detect grid shape + v_step = vtx_dv - vtx_00 + step_tolerance = self.tolerance * la.norm(v_step, np.inf) + for i in range(2, n_points): + step = xy_points[i,:] - xy_points[i - 1, :] + if la.norm(v_step - step, np.inf) > step_tolerance: + break + else: + raise GridNotInShapeExc("End of the first row not detected.") + nv = i + nu = int(n_points / nv) + if not n_points == nu * nv: + raise GridNotInShapeExc("Not a Nu x Nv grid.") + self.shape = (nu, nv) + + # make envelope quad + vtx_01 = xy_points[nv - 1, :] + vtx_10 = xy_points[-nv, :] + vtx_11 = xy_points[-1, :] + self.quad = np.array( [ vtx_01, vtx_00, vtx_10, vtx_11 ], dtype = float ) + self._orig_quad = self.quad + + # check that quad is parallelogram. + diff = np.roll(self.quad, -1, axis = 0) - self.quad + if not la.norm(diff[0] + diff[2]) < self.tolerance * la.norm(diff[0]) or \ + not la.norm(diff[1] + diff[3]) < self.step_tolerance * la.norm(diff[1]): + raise GridNotInShapeExc("Grid XY envelope is not a parallelogram.") + + self._point_tol = self.tolerance * np.max( la.norm(diff, axis = 1) ) + self._u_step = diff[1] / (nu-1) # v10 - v00 + self._v_step = diff[2] / (nv-1) # v11 - v10 + + + def _check_grid_regularity(self, points_xy): + # check regularity of the grid + nu, nv = self.shape + for i in range(nu * nv): + pred_y = i - 1 + pred_x = i - nv + if i%nv == 0: + pred_y = -1 + if pred_x > 0 and not la.norm(points_xy[i, :] - points_xy[pred_x, :] - self._u_step) < self._point_tol: + raise IrregularGridExc("Irregular grid in X direction, point %d"%i) + if pred_y > 0 and not la.norm(points_xy[i, :] - points_xy[pred_y, :] - self._v_step) < self._point_tol: + raise IrregularGridExc("Irregular grid in Y direction, point %d"%i) + + def _make_z_surface(self, points): + nu, nv = self.shape + u_basis = SplineBasis.make_equidistant(1, nu - 1) + v_basis = SplineBasis.make_equidistant(1, nv - 1) + + + + poles_z = points[:, 2].reshape(nu, nv, 1) + self.z_surface = Z_Surface(self.quad[0:3], Surface((u_basis, v_basis), poles_z) ) + self.u_basis = self.z_surface.u_basis + self.v_basis = self.z_surface.v_basis + + uv_points = self.z_surface.xy_to_uv( points[:, 0:2] ) + grid_uv = uv_points.reshape(nu, nv, 2) + self.grid_uvz = np.concatenate((grid_uv, poles_z), axis=2) + + # related bilinear Z-surface, all evaluations just call this object. + + + def _check_map(self): + # check that xy_to_uv works fine + uv_quad = self.xy_to_uv(self.quad) + #print( "Check quad: ", uv_quad ) + assert np.allclose( uv_quad, np.array([[0, 1], [0, 0], [1, 0], [1, 1]]) ) + + def reset_transform(self): + """ + Set identify transform just as after construction. + :param xy_mat: np array, 2 rows 3 cols, last column is xy shift + :param z_mat: [ z_scale, z_shift] + :return: + """ + self.z_scale = 1.0 + self.z_shift = 0.0 + self.z_surface.reset_transform(self._orig_quad) + self.quad = self._orig_quad + self._check_map() + + + + def transform(self, xy_mat, z_mat): + """ + Transform the GridSurface by arbitrary XY linear transform and Z linear transform. + :param xy_mat: np array, 2 rows 3 cols, last column is xy shift + :param z_shift: [ z_scale, z_shift] + :return: None + Note: + """ + # just XY transfrom of underlaying z_surface + if xy_mat is not None: + self.z_surface.transform(xy_mat) + # transform quad + self.quad = self.z_surface.uv_to_xy( np.array([[0, 1], [0, 0], [1, 0], [1, 1]]) ) + + # transform z_scale + if z_mat is not None: + self.z_scale *= z_mat[0] + self.z_shift = z_mat[0]*self.z_shift + z_mat[1] + + + def xy_to_uv(self, xy_points): + """ + :param xy_points: numpy matrix 2 rows N cols + :return: matrix of UV coordinates + """ + return self.z_surface.xy_to_uv(xy_points) + + + def uv_to_xy(self, uv_points): + """ + :param xy_points: numpy matrix 2 rows N cols + :return: matrix of UV coordinates + """ + return self.z_surface.uv_to_xy(uv_points) + + + def eval_array(self, uv_points): + xyz = self.z_surface.eval_array(uv_points) + xyz[:, 2] *= self.z_scale + xyz[:, 2] += self.z_shift + return xyz + + + def z_eval_xy_array(self, xy_points): + """ + Return np array of z-values. Optimized version. + :param points: np array N x 2 of XY cooridinates + :return: array of Z values + """ + assert xy_points.shape[1] == 2, "Size: {}".format(xy_points.shape) + uv_points = self.z_surface.xy_to_uv(xy_points) + return self.z_eval_array(uv_points) + + + def z_eval_array(self, uv_points): + """ + Return np array of z-values. Optimized version. + :param points: np array N x 2 of XY cooridinates + :return: array of Z values + """ + + assert uv_points.shape[1] == 2 + + result = np.zeros(uv_points.shape[0]) + for i, uv in enumerate(uv_points): + iuv = np.floor(uv / self._uv_step) + iu = max(0, min(self.shape[0] - 2, int(iuv[0]))) + iv = max(0, min(self.shape[1] - 2, int(iuv[1]))) + iuv = np.array([iu, iv]) + + uv_loc = uv / self._uv_step - iuv + u_loc = np.array([1 - uv_loc[0], uv_loc[0]]) + v_loc = np.array([1 - uv_loc[1], uv_loc[1]]) + Z_mat = self.grid_uvz[iu: (iu + 2), iv: (iv + 2), 2] + result[i] = self.z_scale*(u_loc.dot(Z_mat).dot(v_loc)) + self.z_shift + return result + + + def center(self): + return self.eval_array( np.array([ [0.5, 0.5] ]))[0] + + + def aabb(self): + z_surf_box = self.z_surface.aabb() + z_surf_box[:,2] *= self.z_scale + z_surf_box[:, 2] += self.z_shift + return z_surf_box + + +def make_function_grid(fn, nu, nv): + """ + Make a grid of points on a graph of the function. + :param fn: fn( [x, y] ) -> z + :param nu: n-points in u-direction + :param nv: n-points in v-direction + :return: array of points: nu x nv x 3 + """ + X_grid = np.linspace(0, 1.0, nu) + Y_grid = np.linspace(0, 1.0, nv) + Y, X = np.meshgrid(Y_grid, X_grid) + + points_uv = np.stack([X.ravel(), Y.ravel()], 1) + Z = np.apply_along_axis(fn, 1, points_uv) + points = np.stack([X.ravel(), Y.ravel(), Z], 1) + + return points.reshape( (nu, nv, 3) ) + + + + + + + + + + diff --git a/src/bspline_approx.py b/src/bspline_approx.py new file mode 100644 index 0000000..1c25a96 --- /dev/null +++ b/src/bspline_approx.py @@ -0,0 +1,664 @@ +""" +Collection of functions to produce Bapline curves and surfaces as approximation of various analytical curves and surfaces. +""" + +import logging +import time +#import math + +import numpy as np +import numpy.linalg as la +import scipy.sparse +import scipy.sparse.linalg +import scipy.interpolate + +import bspline as bs +import csv + +#logging.basicConfig(level=logging.DEBUG) +#logging.info("Test info mesg.") +""" +Approximation methods for B/splines of degree 2. + +""" +def plane_surface(vtxs, overhang=0.0): + """ + Returns B-spline surface of a plane given by 3 points. + We retun also list of UV coordinates of the given points. + U direction v0 -> v1 + V direction v0 -> v2 + :param vtxs: List of tuples (X,Y,Z) + :return: ( Surface, vtxs_uv ) + """ + assert len(vtxs) == 3, "n vtx: {}".format(len(vtxs)) + vtxs = np.array(vtxs) + vv = vtxs[1] + vtxs[2] - vtxs[0] + vtx4 = [ vtxs[0], vtxs[1], vv, vtxs[2]] + return bilinear_surface(vtx4, overhang) + + + +def bilinear_surface(vtxs, overhang=0.0): + """ + Returns B-spline surface of a bilinear surface given by 4 corner points: + uv coords: + We retun also list of UV coordinates of the given points. + :param vtxs: List of tuples (X,Y,Z) + :return: ( Surface, vtxs_uv ) + """ + assert len(vtxs) == 4, "n vtx: {}".format(len(vtxs)) + vtxs = np.array(vtxs) + if overhang > 0.0: + dv = np.roll(vtxs, -1, axis=0) - vtxs + dv *= overhang + vtxs += np.roll(dv, 1, axis=0) - dv + + def mid(*idx): + return np.mean( vtxs[list(idx)], axis=0) + + # v - direction v0 -> v2 + # u - direction v0 -> v1 + poles = [ [vtxs[0], mid(0, 3), vtxs[3]], + [mid(0,1), mid(0,1,2,3), mid(2,3)], + [vtxs[1], mid(1,2), vtxs[2]] + ] + knots = 3 * [0.0 - overhang] + 3 * [1.0 + overhang] + basis = bs.SplineBasis(2, knots) + surface = bs.Surface((basis, basis), poles) + #vtxs_uv = [ (0, 0), (1, 0), (1, 1), (0, 1) ] + return surface + + + + +def line(vtxs, overhang = 0.0): + ''' + Return B-spline approximation of a line from two points + :param vtxs: [ X0, X1 ], Xn are point coordinates in arbitrary dimension D + :return: Curve2D + ''' + assert len(vtxs) == 2 + vtxs = np.array(vtxs) + if overhang > 0.0: + dv = overhang*(vtxs[1] - vtxs[0]) + vtxs[0] -= dv + vtxs[1] += dv + mid = np.mean(vtxs, axis=0) + poles = [ vtxs[0], mid, vtxs[1] ] + knots = 3*[0.0 - overhang] + 3*[1.0 + overhang] + basis = bs.SplineBasis(2, knots) + return bs.Curve(basis, poles) + + + + +def surface_from_grid(grid_surface, nuv): + """ + Make a Z_Surface of degree 2 as an approximation of the GridSurface. + :param grid_surface: grid surface to approximate + :param (nu, nv) Prescribed number of poles in u and v directions. + :return: Z_surface object. + """ + approx = SurfaceApprox.approx_from_grid_surface(grid_surface) + return approx.compute_approximation(nuv=nuv) + + + +def curve_from_grid(points, **kwargs): + """ + Make a Curve (of degree 3) as an approximation of a sequence of points. + :param points - N x D array, D is dimension + :param nt Prescribed number of poles of the resulting spline. + :return: Curve object. + + TODO: + - Measure efficiency. Estimate how good we can be. Do it our self if we can do at leas 10 times better. + - Find out which method is used. Hoschek (4.4.1) refers to several methods how to determine parametrization of + the curve, i.e. Find parameters t_i to the given approximation points P_i. + - Further on it is not clear what is the mening of the 's' parameter and how one cna influence tolerance and smoothness. + - Some sort of adaptivity is used. + + + """ + deg = kwargs.get('degree', 3) + tol = kwargs.get('tol', 0.01) + weights = np.ones(points.shape[0]) + weights[0] = weights[-1] = 1000.0 + tck = scipy.interpolate.splprep(points.T, k=deg, s=tol, w = weights)[0] + knots, poles, degree = tck + curve_poles=np.array(poles).T + curve_poles[0] = points[0] + curve_poles[-1] = points[-1] + basis = bs.SplineBasis(degree, knots) + curve = bs.Curve(basis, curve_poles) + return curve + + + +def convex_hull_2d(sample): + link = lambda a, b: np.concatenate((a, b[1:])) + edge = lambda a, b: np.concatenate(([a], [b])) + + def dome(sample, base): + """ + Return convex hull of the points on the right side from the base. + :param sample: Nx2 nupy array of points + :param base: A segment, np array [[x0,y0], [x1,y1]] + :return: np array of points Nx2 on forming the convex hull + """ + # End points of line. + h, t = base + normal = np.dot(((0,-1),(1,0)),(t-h)) + # Distances from the line. + dists = np.dot(sample-h, normal) + + outer = np.repeat(sample, dists>0, 0) + if len(outer): + pivot = sample[np.argmax(dists)] + return link(dome(outer, edge(h, pivot)), + dome(outer, edge(pivot, t))) + else: + return base + + if len(sample) > 2: + axis = sample[:,0] + # Get left most and right most points. + base = np.take(sample, [np.argmin(axis), np.argmax(axis)], 0) + return link(dome(sample, base), dome(sample, base[::-1])) + else: + return sample + + +def min_bounding_rect(hull): + """ + Compute minimal area bounding box from a convex hull. + Quadratic algorithm with respect to number of hull points is used, anyway calculation a convex hull + takes longer since number of hull points is about sqrt of all points. + :param hull: Nx2 numpy array of the convex hull points. First and last must be the same. + :return: Corners of the rectangle. + """ + # Compute edges (x2-x1,y2-y1) + edges = hull[1:, :] - hull[:-1, :] + + # Calculate edge angles atan2(y/x) + edge_angles = np.arctan2(edges[:, 1], edges[:, 0]) + + # Check for angles in 1st quadrant + edge_angles = np.abs( edge_angles%(np.pi/2)) + + # Remove duplicate angles + edge_angles = np.unique(edge_angles) + + # Test each angle to find bounding box with smallest area + min_bbox = (0, float("inf"), 0, 0, 0, 0, 0, 0) # rot_angle, area, width, height, min_x, max_x, min_y, max_y + for i in range( len(edge_angles) ): + + # Create rotation matrix to shift points to baseline + # R = [ cos(theta) , cos(theta-PI/2) + # cos(theta+PI/2) , cos(theta) ] + angle = edge_angles[i] + R = np.array([[np.cos(angle), np.cos(angle - (np.pi / 2))], + [np.cos(angle + (np.pi / 2)), np.cos(angle)]]) + + # Apply this rotation to convex hull points + rot_points = np.dot(R, np.transpose(hull)) # 2x2 * 2xn + + # Find min/max x,y points + min_x = np.nanmin(rot_points[0], axis=0) + max_x = np.nanmax(rot_points[0], axis=0) + min_y = np.nanmin(rot_points[1], axis=0) + max_y = np.nanmax(rot_points[1], axis=0) + + # Calculate height/width/area of this bounding rectangle + area = (max_x - min_x) * (max_y - min_y) + + # Store the smallest rect found first (a simple convex hull might have 2 answers with same area) + if (area < min_bbox[1]): + min_bbox = ( edge_angles[i], area, min_x, max_x, min_y, max_y ) + + # Re-create rotation matrix for smallest rect + angle = min_bbox[0] + R = np.array([[np.cos(angle), np.cos(angle- (np.pi / 2))], + [np.cos(angle + (np.pi / 2)), np.cos(angle)]]) + + + # min/max x,y points are against baseline + min_x = min_bbox[2] + max_x = min_bbox[3] + min_y = min_bbox[4] + max_y = min_bbox[5] + + # Calculate corner points and project onto rotated frame + corner_points = np.zeros( (4,2) ) # empty 2 column array + corner_points[0] = np.dot( [ min_x, max_y ], R ) + corner_points[1] = np.dot( [ min_x, min_y ], R ) + corner_points[2] = np.dot( [ max_x, min_y ], R ) + corner_points[3] = np.dot( [ max_x, max_y ], R ) + + return corner_points + + +class SurfaceApprox: + """ + Class to compute a Bspline surface approximation from given set of XYZ points. + TODO: + - Check efficiency of scipy methods, compare it to our approach assuming theoretical number of operations. + - Compute BtB directly during single assembly pass, local 9x9 matricies as in A matrix. + - In contradiction to some literature (Hoschek) solution of the LS system is fast as long as the basis is local ( + this is true for B-splines). + - Extensions to fitting X and Y as well - general Surface + + """ + + @staticmethod + def approx_from_file(filename): + """ + Load a sequence of XYZ points on a surface to be approximated. + Optionally points may have weights (i.e. four values per line: XYZW) + :param filename: Path to the input text file. + :return: The approximation object. + """ + with open(filename, 'r') as f: + point_seq = np.array([l for l in csv.reader(f, delimiter=' ')], dtype=float) + + # too slow: alternatives: loadtxt (16s), csv.reader (1.6s), pandas. read_csv (0.6s) + #point_seq = np.loadtxt(filename) + return SurfaceApprox(point_seq) + + + @staticmethod + def approx_from_grid_surface(grid_surface): + """ + Approximation from a GrodSurface object. Use grid of Z coords in + XY pozitions of poles. + :param grid_surface: GridSurface. + :return: + """ + u_basis, v_basis = grid_surface.u_basis, grid_surface.v_basis + + u_coord = u_basis.make_linear_poles() + v_coord = v_basis.make_linear_poles() + + U, V = np.meshgrid(u_coord, v_coord) + uv_points = np.stack( [U.ravel(), V.ravel()], axis = 1 ) + + xyz = grid_surface.eval_array(uv_points) + approx = SurfaceApprox(xyz) + approx.quad = grid_surface.quad + return approx + + def __init__(self, points): + """ + Initialize the approximation object with the points. + :param points: Nx3 (XYZ) or Nx4 (XYZW - points with weights) + """ + + # Degree of approximation in U anv V directions, fixed to 2. + self._degree = np.array((2, 2)) + + + assert( points.shape[1] >= 3 ) + # XYZ points + self._n_points = points.shape[0] + self._xy_points = points[:, 0:2] + self._z_points = points[:, 2] + + # point weights + if points.shape[1] > 3: + self._weights = points[:, 3] + else: + self._weights = None + + ## Approximation parameters. + + # Bounding quadrilateral of the approximation (currently only parallelograms are supported). + # Only first three points P0, P1, P2 are considered. V direction is P0 - P1, U direction is P2 - P1. + # I.e. points are sorted counter-clockwise. + self.quad = None + + # (nu, nv) number of subintervals of the BSSurface on U and V axis. + # Default is estimated from the number of input points N as nu=nv=sqrt(N)/3. + self.nuv = None + + # Weight of the regularizing term. + self.regularization_weight = 0.001 + + ## Approximation results + + # Approximationg BSSurface + self.surface = None + + # Error of the approximation + self.error = None + + def set_quad(self, quad = None): + if quad is None: + quad = np.array([[0,1], [0,0], [1,0], [1,1]]) + self.quad = quad + + def compute_default_quad(self): + """ + Compute and set boundary quad as a minimum area bounding box of the input XY point set. + :return: The quadrilateral vertices. + """ + hull = convex_hull_2d(self._xy_points) + self.quad = min_bounding_rect(hull) + return self.quad + + + def transformed_quad(self, xy_mat): + """ + Return actual quad transformed by given transform matrix. + Boudary quadrilateral of the approximation is not touched. + :param xy_mat: np array, 2 rows 3 cols, last column is xy shift + :return: transformed quad as 4x2 numpy array or None + """ + if self.quad is None: + return None + assert xy_mat.shape == (2, 3) + + # transform quad + return np.dot(self.quad, xy_mat[0:2, 0:2].T) + xy_mat[0:2, 2] + + + def compute_default_nuv(self): + """ + Compute default quad (if not set) filter points in the quad and estimate + nuv from their count. Set self.nuv + :return: nuv = (nu, nv) + """ + if self.quad is None: + self.quad = self.compute_default_quad() + self._compute_uv_points() + + nuv = self._compute_default_nuv(len(self._z_quad_points)) + self.nuv = nuv.astype(int) + if self.nuv[0] < 1 or self.nuv[1] < 1: + raise Exception("Two few points, {}, to make approximation, degree: {}".format(self._n_points, self._degree)) + return self.nuv + + + + + def compute_approximation(self, **kwargs): + """ + Compute approximation of the point set (given to constructor). + Approximation parameters can be passed in through kwargs or set in the object before the call. + :return: B-Spline surface + """ + + self.quad = kwargs.get("quad", self.quad) + self.nuv = kwargs.get("nuv", self.nuv) + self.regularization_weight = kwargs.get("regularization_weight", self.regularization_weight) + + logging.info('Transforming points (n={}) ...'.format(self._n_points)) + start_time = time.time() + if self.quad is None: + self.compute_default_quad() + if self.nuv is None: + self.compute_default_nuv() + + # TODO: better logic, since this has to be recomputed only if quad is changed. + self._compute_uv_points() + + logging.info("Using {} x {} B-spline approximation.".format(self.nuv[0], self.nuv[1])) + self._u_basis = bs.SplineBasis.make_equidistant(2, self.nuv[0]) + self._v_basis = bs.SplineBasis.make_equidistant(2, self.nuv[1]) + + end_time = time.time() + logging.info('Computed in: {} s'.format(end_time - start_time)) + + # Approximation itself + logging.info('Creating B matrix ...') + start_time = time.time() + b_mat, interval = self._build_ls_matrix() + end_time = time.time() + logging.info('Computed in: {} s'.format(end_time - start_time)) + + logging.info('Creating A matrix ...') + start_time = time.time() + a_mat = self._build_sparse_reg_matrix() + end_time = time.time() + logging.info('Computed in: {} s'.format(end_time - start_time)) + + logging.info('Scaling + B^T B ...') + start_time = time.time() + g_vec = self._z_quad_points[:] + if self._weights is not None: + W = scipy.sparse.diags(self._w_quad_points, 0) + wg_vec = np.dot(W, g_vec) + wb_mat = np.dot(W, b_mat) + else: + wg_vec = g_vec + wb_mat = b_mat + b_vec = b_mat.transpose().dot( wg_vec ) + bb_mat = b_mat.transpose().dot(wb_mat) + end_time = time.time() + logging.info('Computed in: {} s'.format(end_time - start_time)) + + logging.info('Computing A and B svds approximation ...') + start_time = time.time() + bb_norm = scipy.sparse.linalg.svds(bb_mat, k=1, ncv=10, tol=1e-2, which='LM', v0=None, + maxiter=300, return_singular_vectors=False) + a_norm = scipy.sparse.linalg.eigsh(a_mat, k=1, ncv=10, tol=1e-2, which='LM', + maxiter=300, return_eigenvectors=False) + c_mat = bb_mat + self.regularization_weight * (bb_norm[0] / a_norm[0]) * a_mat + end_time = time.time() + logging.info('Computed in: {} s'.format(end_time - start_time)) + + logging.info('Solving for Z coordinates ...') + start_time = time.time() + z_vec = scipy.sparse.linalg.spsolve(c_mat, b_vec) + assert not np.isnan(np.sum(z_vec)), "Singular matrix for approximation." + end_time = time.time() + logging.info('Computed in: {} s'.format(end_time - start_time)) + + logging.info('Computing error ...') + start_time = time.time() + diff = b_mat.dot(z_vec) - g_vec + self.error = max_diff = np.max(diff) + logging.info("Approximation error (max norm): {}".format(max_diff) ) + end_time = time.time() + logging.info('Computed in: {} s'.format(end_time - start_time)) + + # Construct Z-Surface + poles_z = z_vec.reshape(self._v_basis.size, self._u_basis.size).T + #poles_z *= self.grid_surf.z_scale + #poles_z += self.grid_surf.z_shift + surface_z = bs.Surface((self._u_basis, self._v_basis), poles_z[:,:,None]) + self.surface = bs.Z_Surface(self.quad[0:3], surface_z) + + return self.surface + + + def _compute_default_nuv(self, n_points): + """ + Default nu and nv for given number of points inside of quad. + :return: (nu, nv) + """ + assert(self.quad is not None) + + dv = la.norm(self.quad[0, :] - self.quad[1, :]) + du = la.norm(self.quad[2, :] - self.quad[1, :]) + + # try to make number of unknowns less then number of remaining points + # +1 to improve determination + nv = np.sqrt( n_points * dv / du ) + nu = nv * du / dv + nuv = np.array( [np.floor(nu / 3), np.floor(nv / 3)] ) - self._degree + self.nuv = np.maximum(1, nuv) + return self.nuv + + + + def _compute_uv_points(self): + """ + Map XY points to quad, remove points out of quad. + Results: self._uv_quad_points, self._z_quad_points, self._w_quad_points + :return: + """ + xy_shift = self.quad[1,:] + v_vec = self.quad[0, :] - self.quad[1, :] + u_vec = self.quad[2, :] - self.quad[1, :] + mat_uv_to_xy = np.column_stack((u_vec, v_vec)) + mat_xy_to_uv = la.inv(mat_uv_to_xy) + points_uv = np.dot((self._xy_points - xy_shift), mat_xy_to_uv.T) + + # remove points far from unit square + eps = 1.0e-15 + cut_min = np.array( [ -eps, -eps ]) + cut_max = np.array( [ 1+eps, 1+eps ]) + in_idx = np.all(np.logical_and(cut_min < points_uv, points_uv <= cut_max), axis=1) + points_uv = points_uv[in_idx] + + logging.debug("Number of points out of the grid domain: {}".format(len(points_uv) - np.sum(in_idx))) + + # snap to unit square + points_uv = np.maximum(points_uv, np.array([0.0, 0.0])) + self._uv_quad_points = np.minimum(points_uv, np.array([1.0, 1.0])) + self._z_quad_points = self._z_points[in_idx] + if self._weights is not None: + self._w_quad_points = self._weights[in_idx] + + + + + + + def _build_ls_matrix(self): + """ + Construction of the matrix (B) of the system of linear algebraic + equations for control points of the 2th order B-spline surface + :param u_knots: + :param v_knots: + :param terrain: + :param sparse: + :return: + """ + u_n_basf = self._u_basis.size + v_n_basf = self._v_basis.size + n_points = self._uv_quad_points.shape[0] + + n_uv_loc_nz = (self._u_basis.degree + 1) * (self._v_basis.degree + 1) + row = np.zeros(n_points * n_uv_loc_nz) + col = np.zeros(n_points * n_uv_loc_nz) + data = np.zeros(n_points * n_uv_loc_nz) + + nnz_b = 0 + + interval = np.empty((n_points, 2)) + + for idx in range(0, n_points): + u, v = self._uv_quad_points[idx, 0:2] + iu = self._u_basis.find_knot_interval(u) + iv = self._v_basis.find_knot_interval(v) + u_base_vec = self._u_basis.eval_base_vector(iu, u) + v_base_vec = self._v_basis.eval_base_vector(iv, v) + # Hard-coded Kronecker product (problem based) + for n in range(0, 3): + data[nnz_b + 3 * n:nnz_b + 3 * (n + 1)] = v_base_vec[n] * u_base_vec + for m in range(0, 3): + col_item = (iv + n) * u_n_basf + iu + m + col[nnz_b + (3 * n) + m] = col_item + row[nnz_b:nnz_b + 9] = idx + nnz_b += 9 + + interval[idx][0] = iu + interval[idx][1] = iv + + mat_b = scipy.sparse.csr_matrix((data, (row, col)), shape=(n_points, u_n_basf * v_n_basf)) + + return mat_b, interval + + def _basis_in_q_points(self, basis): + n_int = basis.n_intervals + nq_points = len(self._q_points) + q_point = np.zeros((n_int * nq_points, 1)) + point_val_outer = np.zeros((3, 3, n_int)) # "3" considers degree 2 + d_point_val_outer = np.zeros((3, 3, n_int)) # "3" considers degree 2 + + #TODO: use numpy functions for quadrature points + n = 0 + for i in range(n_int): + us = basis.knots[i + 2] + uil = basis.knots[i + 3] - basis.knots[i + 2] + for j in range(nq_points): + up = us + uil * self._q_points[j] + q_point[n] = up + u_base_vec = basis.eval_base_vector(i, up) + u_base_vec_diff = basis.eval_diff_base_vector(i, up) + point_val_outer[:, :, i] += self._q_weights[j] * np.outer(u_base_vec,u_base_vec) + d_point_val_outer[:, :, i] += self._q_weights[j] * np.outer(u_base_vec_diff,u_base_vec_diff) + n += 1 + + return point_val_outer, d_point_val_outer,q_point + + + def _build_sparse_reg_matrix(self): + """ + Construction of the regularization matrix (A) to decrease variation of the terrain + B z = b ---> (B^T B + A)z = B^T b + :param u_knots: vector of v-knots + :param v_knots: vector of u-knots + :param quad: points defining quadrangle area (array) + :return: matrix + + - + """ + + #a = quad[:, 3] - quad[:, 2] + #b = quad[:, 0] - quad[:, 1] + #c = quad[:, 1] - quad[:, 2] + #d = quad[:, 0] - quad[:, 3] + + u_n_basf = self._u_basis.size + v_n_basf = self._v_basis.size + u_n_inter = self._u_basis.n_intervals + v_n_inter = self._v_basis.n_intervals + n_uv_loc_nz = (self._u_basis.degree + 1) * (self._v_basis.degree + 1) + + # TODO: use Gauss quadrature from scipy + # in fact for general degrees we should use different quadrature for u and different for v + self._q_points = [0, (0.5 - 1 / np.sqrt(20)), (0.5 + 1 / np.sqrt(20)), 1] + self._q_weights = [1.0 / 6, 5.0 / 6, 5.0 / 6, 1.0 / 6] + nq_points = len(self._q_points) + + u_val_outer, u_diff_val_outer, q_u_point = self._basis_in_q_points(self._u_basis) + v_val_outer, v_diff_val_outer, q_v_point = self._basis_in_q_points(self._v_basis) + # xy_outer shape is (3, 3, n_inter) + + row_m = np.zeros((v_n_inter * u_n_inter * n_uv_loc_nz * n_uv_loc_nz)) + col_m = np.zeros((v_n_inter * u_n_inter * n_uv_loc_nz * n_uv_loc_nz)) + data_m = np.zeros((v_n_inter * u_n_inter * n_uv_loc_nz * n_uv_loc_nz)) + + nnz_a = 0 + #linsp = np.linspace(0, self._u_basis.degree, self._u_basis.degree+1) + #llinsp = np.tile(linsp, self._u_basis.degree+1) + #np.repeat((iv + linsp) * u_n_basf, self._u_basis.degree + 1) + llinsp + i_local = np.arange(self._u_basis.degree+1, dtype=int) + iuv_local = (u_n_basf * i_local[:, None] + i_local[None,:]).ravel() # 0,1,2, N+[0,1,2], 2*N+[0,1,2] + #print("vnint: {} unint: {} nqp: {} prod: {}".format(v_n_inter, u_n_inter, nq_points, v_n_inter* u_n_inter* nq_points*nq_points)) + jac = 1.0 / u_n_inter / v_n_inter + idx_range = n_uv_loc_nz * n_uv_loc_nz # 9 * 9 = 81 NZ per single bspline square + for iv in range(v_n_inter): + v_val_outer_loc = v_val_outer[:, :, iv] + dv_val_outer_loc = v_diff_val_outer[:, :, iv] + + for iu in range(u_n_inter): + u_val_outer_loc = u_val_outer[:, :, iu] + du_val_outer_loc = u_diff_val_outer[:, : , iu] + # xy_outer_loc have shape 3x3 + + v_du = np.kron(v_val_outer_loc, du_val_outer_loc) + dv_u = np.kron(dv_val_outer_loc, u_val_outer_loc) + data_m[nnz_a:nnz_a + idx_range] = jac * ( v_du + dv_u).ravel() # 9x9 values + + iuv = iu + iv * u_n_basf + colv = iuv + iuv_local + col_m[nnz_a:nnz_a + idx_range] = np.repeat(colv, n_uv_loc_nz) + row_m[nnz_a:nnz_a + idx_range] = np.tile(colv, n_uv_loc_nz) + nnz_a += idx_range + #print("Assembled") + mat_a = scipy.sparse.coo_matrix((data_m, (row_m, col_m)), + shape=(u_n_basf * v_n_basf, u_n_basf * v_n_basf)).tocsr() + return mat_a \ No newline at end of file diff --git a/src/bspline_plot.py b/src/bspline_plot.py new file mode 100644 index 0000000..06ecb42 --- /dev/null +++ b/src/bspline_plot.py @@ -0,0 +1,254 @@ +""" +Functions to plot Bspline curves and surfaces. +""" + +plot_lib = "plotly" + +import plotly.offline as pl +import plotly.graph_objs as go + +import matplotlib.pyplot as plt +from matplotlib import cm +from mpl_toolkits.mplot3d import Axes3D + +import numpy as np + + + + +class PlottingPlotly: + def __init__(self): + self.i_figure = -1 + self._reinit() + + def _reinit(self): + self.i_figure += 1 + self.data_3d = [] + self.data_2d = [] + + def add_curve_2d(self, X, Y, **kwargs): + self.data_2d.append( go.Scatter(x=X, y=Y, mode = 'lines') ) + + def add_points_2d(self, X, Y, **kwargs): + marker = dict( + size=10, + color='red', + ) + self.data_2d.append( go.Scatter(x=X, y=Y, + mode = 'markers', + marker=marker) ) + + + def add_surface_3d(self, X, Y, Z, **kwargs): + hue = (120.0*(len(self.data_3d)))%360 + colorscale = [[0.0, 'hsv({}, 50%, 10%)'.format(hue)], [1.0, 'hsv({}, 50%, 90%)'.format(hue)]] + self.data_3d.append( go.Surface(x=X, y=Y, z=Z, colorscale=colorscale)) + + + def add_points_3d(self, X, Y, Z, **kwargs): + marker = dict( + size=5, + color='red', + # line=dict( + # color='rgba(217, 217, 217, 0.14)', + # width=0.5 + # ), + opacity=0.6 + ) + self.data_3d.append( go.Scatter3d( + x=X, y=Y, z=Z, + mode='markers', + marker=marker + )) + + + def show(self): + """ + Show added plots and clear the list for other plotting. + :return: + """ + if self.data_3d: + fig_3d = go.Figure(data=self.data_3d) + pl.plot(fig_3d, filename='bc_plot_3d_%d.html'%(self.i_figure)) + if self.data_2d: + fig_2d = go.Figure(data=self.data_2d) + pl.plot(fig_2d, filename='bc_plot_2d_%d.html'%(self.i_figure)) + self._reinit() + + +class PlottingMatplot: + def __init__(self): + self.fig_2d = plt.figure(1) + self.fig_3d = plt.figure(2) + self.ax_3d = self.fig_3d.gca(projection='3d') + + def add_curve_2d(self, X, Y, **kwargs): + plt.figure(1) + plt.plot(X, Y, **kwargs) + + def add_points_2d(self, X, Y, **kwargs): + plt.figure(1) + plt.plot(X, Y, 'bo', color='red', **kwargs) + + def add_surface_3d(self, X, Y, Z, **kwargs): + plt.figure(2) + self.ax_3d.plot_surface(X, Y, Z, **kwargs) + + def add_points_3d(self, X, Y, Z, **kwargs): + plt.figure(2) + return self.ax_3d.scatter(X, Y, Z, color='red', **kwargs) + + def show(self): + """ + Show added plots and clear the list for other plotting. + :return: + """ + plt.show() + + +class Plotting: + """ + Debug plotting class. Several 2d and 3d plots can be added and finally displayed on common figure + calling self.show(). Matplotlib or plotly library is used as backend. + """ + def __init__(self, backend = PlottingPlotly()): + self.backend = backend + + def plot_2d(self, X, Y): + """ + Add line scatter plot. Every plot use automatically different color. + :param X: x-coords of points + :param Y: y-coords of points + """ + self.backend.add_curve_2d(X,Y) + + def scatter_2d(self, X, Y): + """ + Add point scatter plot. Every plot use automatically different color. + :param X: x-coords of points + :param Y: y-coords of points + """ + self.backend.add_points_2d(X,Y) + + def plot_surface(self, X, Y, Z): + """ + Add line scatter plot. Every plot use automatically different color. + :param X: x-coords of points + :param Y: y-coords of points + """ + self.backend.add_surface_3d(X, Y, Z) + + def plot_curve_2d(self, curve, n_points=100, poles=False): + """ + Add plot of a 2d Bspline curve. + :param curve: Curve t -> x,y + :param n_points: Number of evaluated points. + :param: kwargs: Additional parameters passed to the mtplotlib plot command. + """ + + basis = curve.basis + t_coord = np.linspace(basis.domain[0], basis.domain[1], n_points) + + coords = [curve.eval(t) for t in t_coord] + x_coord, y_coord = zip(*coords) + + self.backend.add_curve_2d(x_coord, y_coord) + if poles: + self.plot_curve_poles_2d(curve) + + + def plot_curve_poles_2d(self, curve): + """ + Plot poles of the B-spline curve. + :param curve: Curve t -> x,y + :return: Plot object. + """ + x_poles, y_poles = curve.poles.T[0:2, :] # remove weights + return self.backend.add_points_2d(x_poles, y_poles) + + def scatter_3d(self, X, Y, Z): + """ + Add point scatter plot. Every plot use automatically different color. + :param X: x-coords of points + :param Y: y-coords of points + """ + self.backend.add_points_3d(X, Y, Z) + + + def plot_surface_3d(self, surface, n_points=(100, 100), poles=False): + """ + Plot a surface in 3d. + Usage: + plotting=Plotting() + plotting.plot_surface_3d(surf_1) + plotting.plot_surface_3d(surf_2) + plotting.show() + + :param surface: Parametric surface in 3d. + :param n_points: (nu, nv), nu*nv - number of evaluation point + """ + + u_basis, v_basis = surface.u_basis, surface.v_basis + + u_coord = np.linspace(u_basis.domain[0], u_basis.domain[1], n_points[0]) + v_coord = np.linspace(v_basis.domain[0], v_basis.domain[1], n_points[1]) + + U, V = np.meshgrid(u_coord, v_coord) + points = np.stack( [U.ravel(), V.ravel()], axis = 1 ) + + xyz = surface.eval_array(points) + X, Y, Z = xyz.T + + X = X.reshape(U.shape) + Y = Y.reshape(U.shape) + Z = Z.reshape(U.shape) + + # Plot the surface. + self.backend.add_surface_3d(X, Y, Z) + + if poles: + self.plot_surface_poles_3d(surface) + + + + + def plot_grid_surface_3d(self, surface, n_points=(100, 100)): + """ + Plot a surface in 3d, on UV plane. + + :param surface: Parametric surface in 3d. + :param n_points: (nu, nv), nu*nv - number of evaluation point + """ + + u_coord = np.linspace(0, 1.0, n_points[0]) + v_coord = np.linspace(0, 1.0, n_points[1]) + + U, V = np.meshgrid(u_coord, v_coord) + points = np.stack( [U.ravel(), V.ravel()], axis = 1 ) + + xyz = surface.eval_array(points) + X, Y, Z = xyz.T + + Z = Z.reshape(U.shape) + + # Plot the surface. + self.backend.add_surface_3d(U, V, Z) + + + def plot_surface_poles_3d(self, surface, **kwargs): + """ + Plot poles of the B-spline curve. + :param curve: Curve t -> x,y + :param: kwargs: Additional parameters passed to the mtplotlib plot command. + :return: Plot object. + """ + x_poles, y_poles, z_poles = surface.poles[:, :, 0:3].reshape(-1, 3).T # remove weights and flatten nu, nv + return self.backend.add_points_3d(x_poles, y_poles, z_poles, **kwargs) + + + def show(self): + """ + Display added plots. Empty the queue. + :return: + """ + self.backend.show() diff --git a/src/isec_curv_surf.py b/src/isec_curv_surf.py new file mode 100644 index 0000000..980ba0c --- /dev/null +++ b/src/isec_curv_surf.py @@ -0,0 +1,8 @@ +class IsecCurvSurf: + ''' + Class for calculation and representation of the intersection of + B-spline curve and B-spline surface. + Currently only (multi)-point intersections are assumed. + ''' + def __init__(self, curve, surface): + pass \ No newline at end of file diff --git a/src/isec_surf_surf.py b/src/isec_surf_surf.py new file mode 100644 index 0000000..ee49352 --- /dev/null +++ b/src/isec_surf_surf.py @@ -0,0 +1,746 @@ +import sys + +build_path = "/home/jiri/Soft/Geomop/Intersections/external/bih/build" +sys.path += [build_path] +print(sys.path) + +import bih +import numpy as np +import numpy.linalg as la + +import bspline as bs + + +class IsecPoint: + """ + Point as the result of intersection with correspondind coordinates on both surfaces + """ + def __init__(self, surf1, iu1, iv1, uv1, boundary_flag, flag ,sum_idx, surf2, iu2, iv2, uv2,XYZ): + + self.surf1 = surf1 + self.iu1 = iu1 + self.iv1 = iv1 + self.uv1 = uv1 + self.boundary_flag = boundary_flag + + + self.surf2 = surf2 + self.iu2 = iu2 + self.iv2 = iv2 + self.uv2 = uv2 + + + self.tol = 1 # implement + self.R3_coor = XYZ + + self.connected = 0 + + + if boundary_flag == 1: + self.add_patches(sum_idx,flag[2]) + + if np.logical_or(flag[0]!=-1,flag[1]!=-1): + self.extend_patches( flag[0:2]) + + # flag + + def belongs_to(self,iu,iv,isurf): + + belongs = 0 + if isurf == 0: + for i in range(self.iu1.__len__): + if (np.logical_and(self.iu1[i] == iu),np.logical_and(self.iv1[i] == iv)): + belongs = 1 + break + + return belongs + + #def duplicite(self,point): + # duplicite = 0 + + # if point.iu1.__len__ + + + + + + + # return duplicite + + + def connect(self): + self.connected = 1 + + def extend_patches(self, flag): + + if flag[0]!=-1: + direction_u = 2 * flag[0] - 1 + self.iu2.append(self.iu2[0] + direction_u) + self.iv2.append(self.iv2[0]) + + if flag[1]!=-1: + direction_v = 2 * flag[1] - 1 + self.iu2.append(self.iu2[0]) + self.iv2.append(self.iv2[0]+ direction_v) + + if np.logical_and(flag[0]!=-1,flag[1]!=-1): + self.iu2.append(self.iu2[0]+ direction_u) + self.iv2.append(self.iv2[0]+ direction_v) + + + def add_patches(self, sum_idx,flag): + + n = self.iu1.__len__() + + direction = 2*flag -1 + + if sum_idx == 0: + for i in range(n): + self.iu1.append(self.iu1[0]+direction) + self.iv1.append(self.iv1[i]) + elif sum_idx == 1: + for i in range(n): + self.iu1.append(self.iu1[i]) + self.iv1.append(self.iv1[0]+direction) + +class IsecCurveSurf: + """ + Class which provides intersection of the given surface and curve + """ + + def __init__(self, surf, curv): + self.surf = surf + self.curv = curv + + + def _compute_jacobian_and_delta(self, uvt, iu, iv, it): + """ + Computes Jacobian matrix and delta vector of the function + :param uvt: vector of local coordinates [u,v,t] (array 3x1) + :param iu: index of the knot interval of the coordinate u + :param iv: index of the knot interval of the coordinate v + :param it: index of the knot interval of the coordinate t + :return: J: jacobian matrix (array 3x3) , deltaXYZ: vector of deltas in R^3 space (array 3x1) + """ + + surf = self.surf + curv = self.curv + + surf_poles = surf.poles[iu:iu + 3, iv:iv + 3, :] + t_poles = curv.poles[it:it + 3, :] + + + uf = surf.u_basis.eval_vector(iu, uvt[0, 0]) + vf = surf.v_basis.eval_vector(iv, uvt[1, 0]) + ufd = surf.u_basis.eval_diff_vector(iu, uvt[0, 0]) + vfd = surf.v_basis.eval_diff_vector(iv, uvt[1, 0]) + tf = curv.basis.eval_vector(it, uvt[2, 0]) + tfd = curv.basis.eval_diff_vector(it, uvt[2, 0]) + + dXYZt = np.tensordot(tfd, t_poles, axes=([0], [0])) + #print(dXYZt) + dXYZu1 = self._energetic_inner_product(ufd, vf, surf_poles) + dXYZv1 = self._energetic_inner_product(uf, vfd, surf_poles) + J = np.column_stack((dXYZu1, dXYZv1, -dXYZt)) + + XYZ1 = self._energetic_inner_product(uf, vf, surf_poles) + XYZ2 = np.tensordot(tf, t_poles, axes=([0], [0])) + #print(XYZ2) + #return + XYZ2 = XYZ2[:, None] + XYZ1 = XYZ1[:, None] + + deltaXYZ = XYZ1 - XYZ2 + + return J, deltaXYZ + + def get_intersection(self, uvt, iu, iv, it, max_it, rel_tol, abs_tol): + """ + Main solving method for solving + :param uvt: vector of initial guess of local coordinates [u,v,t] (array 3x1) + :param iu: index of the knot interval of the coordinate u + :param iv: index of the knot interval of the coordinate v + :param it: index of the knot interval of the coordinate t + :param max_it: maximum number of iteration + :param rel_tol: relative tolerance (absolute in parametric space) + :param abs_tol: absolute tolerance (in R3 space) + :return: uvt: vector of initial guess of local coordinates [u,v,t] (array 3x1), + conv as "0" if the methog does not achive desired accuracy + "1" if the methog achive desired accuracy + flag as intersection specification + """ + conv = 0 + #flag = -1 + flag = -np.ones([3],'int') + + ui = self._compute_bounds(self.surf.u_basis.knots, iu) + vi = self._compute_bounds(self.surf.v_basis.knots, iv) + ti = self._compute_bounds(self.curv.basis.knots, it) + + for i in range(max_it): + J, deltaXYZ = self._compute_jacobian_and_delta(uvt, iu, iv, it) + if la.norm(deltaXYZ) < abs_tol: + break + uvt = uvt - la.solve(J, deltaXYZ) + uvt = self._range_test(uvt, ui, vi, ti, 0.0) + + conv, XYZ = self._test_intesection_tolerance(uvt, iu, iv, it, abs_tol) + + if conv == 1: + flag[0] = self._patch_boundary_intersection(uvt[0, 0], ui, rel_tol) + flag[1] = self._patch_boundary_intersection(uvt[1, 0], vi, rel_tol) + flag[2] = self._curve_boundary_intersection(uvt[2, 0], ti, rel_tol) + + + return uvt, conv, flag, XYZ + + @staticmethod + def _curve_boundary_intersection(t,ti,rel_tol): + """ + + :param t: parameter value + :param it: interval array (2x1) + :return: + flag as "0" corresponds to lower bound of the curve + "1" corresponds to upper bound of the curve + "-1" corresponds to the interior points of the curve + """ + # interior boundary + if abs(ti[0] - t) < rel_tol: + flag = 0 + elif abs(ti[1] - t) < rel_tol: + flag = 1 + else: + flag = -1 + + return flag + + @staticmethod + def _patch_boundary_intersection(t,ti,rel_tol): + """ + + :param t: parameter value + :param it: interval array (2x1) + :return: + flag as "-1" corresponds to lower bound of the curve + "1" corresponds to upper bound of the curve + "0" corresponds to the boundary points of the curve + """ + # interior boundary + + flag = -1 + + if ti[0] != 0: + if abs(ti[0] - t) < rel_tol: + flag = 0 + elif ti[1] != 1: + if abs(ti[1] - t) < rel_tol: + flag = 1 + #else: + # flag = -1 + + return flag + + + @staticmethod + def _compute_bounds(knots, idx): + """ + Computes boundary of the given area (lower and upper) in parametric space + :param knots: knot vector + :param idx: index of the interval + :return: bounds as (array 2x1) (lower bound, upper bound) + """ + s = knots[idx + 2] + e = knots[idx + 3] + bounds = np.array([s, e]) + return bounds + + @staticmethod + def _range_test(uvt, ui, vi, ti, rel_tol): + """ + Test of the entries of uvt, lies in given intervals + with a given tolerance + :param uvt: vector of local coordinates [u,v,t] (array 3x1) + :param ui: knot interval of the coordinate u (array 2x1) + :param vi: knot interval of the coordinate v (array 2x1) + :param ti: knot interval of the coordinate t (array 2x1) + :param rel_tol: given relative tolerance (absolute in [u,v,t]) + :return: uvt: vector of local coordinates [u,v,t] (array 3x1) + """ + + test = 0 + + du = np.array([ ui[0] - uvt[0, 0],uvt[0, 0] - ui[1]]) + dv = np.array([ vi[0] - uvt[1, 0],uvt[1, 0] - vi[1]]) + dt = np.array([ ti[0] - uvt[2, 0],uvt[2, 0] - ti[1]]) + + for i in range(0, 2): + if (du[i] > rel_tol): + uvt[0, 0] = ui[i] + + for i in range(0, 2): + if (dv[i] > rel_tol): + uvt[1, 0] = vi[i] + + for i in range(0, 2): + if (dt[i] > rel_tol): + uvt[2, 0] = ti[i] + + #if np.logical_and(uvt[0, 0] >= ui[0], uvt[0, 0] <= ui[1]): + # if np.logical_and(uvt[1, 0] >= vi[0], uvt[1, 0] <= vi[1]): + # if np.logical_and(uvt[2, 0] >= ti[0], uvt[2, 0] <= ti[1]): + # test = 1 + + return uvt + + def _test_intesection_tolerance(self, uvt, iu, iv, it, abs_tol): + """ + Test of the tolerance of the intersections in R3 + :param uvt: vector of local coordinates [u,v,t] (array 3x1) + :param iu: index of the knot interval of the coordinate u + :param iv: index of the knot interval of the coordinate v + :param it: index of the knot interval of the coordinate t + :param abs_tol: given absolute tolerance in R^3 space + :return: conv as "0" if the methog does not achive desired accuracy + or "1" if the methog achive desired accuracy + """ + + conv = 0 + + surf_R3 = self.surf.eval_local(uvt[0, 0], uvt[1, 0], iu, iv) + curv_R3 = self.curv.eval_local(uvt[2, 0], it) + XYZ = (surf_R3 + curv_R3)/2 + dist = la.norm(surf_R3 - curv_R3) + #print('distance =', dist) + if dist <= abs_tol: + conv = 1 + + return conv, XYZ + + def get_initial_condition(self, iu, iv, it): + """ + Computes initial guess as mean value of the considered area + :param iu: index of the knot interval of the coordinate u + :param iv: index of the knot interval of the coordinate v + :param it: index of the knot interval of the coordinate t + :return: uvt as array (3x1) + """ + + uvt = np.zeros([3, 1]) + + uvt[0, 0] = self._get_mean_value(self.surf.u_basis.knots, iu) + uvt[1, 0] = self._get_mean_value(self.surf.v_basis.knots, iv) + uvt[2, 0] = self._get_mean_value(self.curv.basis.knots, it) + + return uvt + + @staticmethod + def _get_mean_value(knots, int): + """ + Computes mean value of the local coordinate in the given interval + :param knots: knot vector + :param idx: index of the patch + :return: mean + """ + mean = (knots[int + 2] + knots[int + 3])/2 + + return mean + + @staticmethod + def _energetic_inner_product(u, v, surf_poles): + """ + Computes energetic inner product u^T X v + :param u: vector of nonzero basis function in u + :param v: vector of nonzero basis function in v + :param X: tensor of poles in x,y,z + :return: xyz as array (3x1) + """ + uX = np.tensordot(u, surf_poles, axes=([0], [0])) + xyz = np.tensordot(uX, v, axes=([0], [0])) + return xyz + +class IsecSurfSurf: + def __init__(self, surf1, surf2, nt=2, max_it=10, rel_tol = 1e-16, abs_tol = 1e-14): + self.surf1 = surf1 + self.surf2 = surf2 + self.box1, self.tree1 = self.bounding_boxes(self.surf1) + self.box2, self.tree2 = self.bounding_boxes(self.surf2) + self.nt = nt + self.max_it = max_it + self.abs_tol = abs_tol + self.rel_tol = rel_tol + + self._ipoint_list = [] # append + # tolerance + + def bounding_boxes(self, surf): + """ + Compute bounding boxes and construct BIH tree for a given surface + :param surf: + :return: + """ + tree = bih.BIH() + n_patch = (surf.u_basis.n_intervals) * (surf.v_basis.n_intervals) + + patch_poles = np.zeros([9, 3, n_patch]) + i_patch = 0 + for iu in range(surf.u_basis.n_intervals): + for iv in range(surf.v_basis.n_intervals): + n_points = 0 + for i in range(0, 3): + for j in range(0, 3): + patch_poles[n_points, :, i_patch] = surf.poles[iu + i, iv + j, :] + n_points += 1 + #assert i_patch == (iu * surf.v_basis.n_intervals + iv) + assert i_patch == self._patch_pos2id(surf,iu,iv) + i_patch += 1 + + boxes = [bih.AABB(patch_poles[:, :, p].tolist()) for p in range(n_patch)] + # print(patch_poles[:, :, 0]) + tree.add_boxes(boxes) + tree.construct() + # print(boxes) + # for b in boxes: + # print(b.min()[0:2],b.max()[0:2]) + return boxes, tree + + def get_intersection(self): + """ + Main method to get intersection points + :return: + """ + + + point_list1 = self.get_intersections(self.surf1, self.surf2, self.tree2) # patches of surf 2 with respect threads of the surface 1 + point_list2 = self.get_intersections(self.surf2, self.surf1, self.tree1) # patches of surf 1 with respect threads of the surface 2 + + print('points') + for point in point_list1: + print(point.uv1) + print('points2') + for point in point_list2: + print(point.uv2) + print(point_list1.__len__(),point_list2.__len__()) + + assert point_list1.__len__() == 9 + assert point_list2.__len__() == 9 + + return point_list1, point_list2 + + connected_points = self._connect_points(point_list1,point_list2) + + + def _connect_points(self,point_list1,point_list2): + + patch_point_list1, boundary_points1 = self.make_patch_point_list(point_list1,point_list2) + patch_point_list2, boundary_points2 = self.make_patch_point_list(point_list2,point_list1) + + print(boundary_points2) + + patch_point_list = [[],[]] + boundary_points = [[],[]] + + patch_point_list[0] = patch_point_list1 # append + patch_point_list[1] = patch_point_list2 + boundary_points[0] = boundary_points1 + boundary_points[1] = boundary_points2 + + + #line = self._make_orderings(patch_point_list, boundary_points) + + + + + + return patch_point_list2 + + + def _make_orderings(self,patch_point_list, boundary_points): + + n_curves = 0 + + line = [] + + n_intersections = patch_point_list[0].__len__() + patch_point_list[1].__len__() + #get_start + #connect + #test end + +# while n_intersections != 0: +# for i in range(0,2): +# for patchpoint in boundary_points[i]: +# if patchpoint.connected == 0: +# patchpoint.connect() +# n_intersections -= 1 +# line[n_curves].append(patchpoint) +# while end_found == 0: +# #id = self._patch_pos2id(patchpoint.surf2,patchpoint.iu2[k],patchpoint.iv2[k]) +# pt = np.zeros([patchpoint.iu1.___len___]) +# for k in range(patchpoint.iu1.___len___): +# ids = self._patch_pos2id(patchpoint.surf1,patchpoint.iu1[k],patchpoint.iv1[k]) +# for pp in patch_point_list[1-i][ids]: +# if pp.connected == 0: +# pt[k] += 1 +# if np.sum(pt) > 1: +# print('problem') +# else: +# +# #found_start +# end_found = 0 +# while end_found == 0: +# #search next +# check_duplicities + + + + while n_intersections != 0: + for i in range(0,2): + for patchpoint in patch_point_list[i]: + if np.logical_and(patchpoint.connected == 0, patchpoint.boundary_flag == 1 ): + patchpoint.connect() + n_intersections -= 1 + line[n_curves].append(patchpoint) + while end_found == 0: + #id = self._patch_pos2id(patchpoint.surf2,patchpoint.iu2[k],patchpoint.iv2[k]) + pt = np.zeros([patchpoint.iu1.___len___]) + for k in range(patchpoint.iu1.___len___): + ids = self._patch_pos2id(patchpoint.surf1,patchpoint.iu1[k],patchpoint.iv1[k]) + for pp in patch_point_list[1-i][ids]: + if pp.connected == 0: + pt[k] += 1 + if np.sum(pt) > 1: + print('problem') + else: + + #found_start + end_found = 0 + while end_found == 0: + #search next + check_duplicities + + + + + n_curves += 1 + + + + return line + + #@staticmethod + def make_patch_point_list(self,point_list,point_list2): + + surf1 = point_list[0].surf1 + + list_len = surf1.u_basis.n_intervals * surf1.v_basis.n_intervals + patch_points = [] + boundary_points = [] + + boundary_points.append([[],[]]) + + + for i in range(list_len): + patch_points.append([[],[]]) + + idp = -1 + for point in point_list: + idp += 1 + for i_patch in range(point.iu1.__len__()): + id = self._patch_pos2id(surf1,point.iu1[i_patch],point.iv1[i_patch]) + patch_points[id][0].append(idp) + if point.boundary_flag == 1: + boundary_points[0].append(idp) + + idp = -1 + for point in point_list2: + idp += 1 + for i_patch in range(point.iu2.__len__()): + id = self._patch_pos2id(surf1, point.iu2[i_patch], point.iv2[i_patch]) + patch_points[id][1].append(idp) + if point.boundary_flag == 1: + boundary_points[1].append(idp) + + return patch_points, boundary_points + + @staticmethod + def _main_threads(surf,sum_idx): + """ + Construction of the main threads + :param surf: surface which is used to construction of the main threads + :param sum_idx: sum_idx == 0 --> u fixed, sum_idx == 1 --> v fixed + :return: curves as list of curves, w_val as list of value of the fixed local coordinates , patches as list of neighbour patches + """ + + poles = surf.poles + + if sum_idx == 0: + fix_basis = surf.u_basis + curv_basis = surf.v_basis + elif sum_idx == 1: + fix_basis = surf.v_basis + curv_basis = surf.u_basis + + curves = [] + w_val = [] + patches = [] + + if sum_idx == 0: + curv_pol = poles[0, :, :] + elif sum_idx == 1: + curv_pol = poles[:, 0, :] + + w_val.append(0.0) + patches.append([0]) + + curv = bs.Curve(curv_basis, curv_pol) + curves.append(curv) + + for iw in range(1,fix_basis.n_intervals): + w1f = fix_basis.eval_vector(iw, fix_basis.knots[iw + 2]) + if sum_idx == 0: + surf_pol = poles[iw:iw + 3, :, :] + curv_pol = np.tensordot(w1f, surf_pol, axes=([0], [sum_idx])) + elif sum_idx == 1: + surf_pol = poles[:,iw:iw + 3, :] + curv_pol = np.tensordot(w1f, surf_pol, axes=([0], [sum_idx])) + w_val.append(fix_basis.knots[iw + 2]) + patches.append([iw-1, iw]) + + curv = bs.Curve(curv_basis, curv_pol) + curves.append(curv) + + if sum_idx == 0: + curv_pol = poles[poles.shape[sum_idx]-1,:, :] + #print(poles.shape[sum_idx]-1,fix_basis.n_intervals - 1) + elif sum_idx == 1: + curv_pol = poles[:,poles.shape[sum_idx]-1 , :] # fix_basis.n_intervals - 1 + #print(poles.shape[sum_idx] - 1,fix_basis.n_intervals - 1) + + w_val.append(1.0) + patches.append([fix_basis.n_intervals-1]) + + curv = bs.Curve(curv_basis, curv_pol) + curves.append(curv) + + return curves, w_val, patches + + @staticmethod + def _patch_pos2id(surf,iu,iv): + + id = iu * surf.v_basis.n_intervals + iv + return id + + @staticmethod + def _patch_id2pos(surf,id): + + iu = int(np.floor(id / surf.v_basis.n_intervals)) + iv = int(id - (iu * surf.v_basis.n_intervals)) + + return iu, iv + + def get_intersections(self,surf1,surf2,tree2): + """ + Tries to compute intersection of the main curves from surface1 and patches of the surface2 which have a + non-empty intersection of corresponding bonding boxes + :param surf1: Surface used to construction of the main threads + :param surf2: Intersected surface + :param tree2: Bih tree of the patches of the surface 2 + :return: point_list as list of points of intersection + """ + + point_list = [] + crossing = np.zeros([surf1.u_basis.n_intervals+1,surf1.v_basis.n_intervals+1]) + + for sum_idx in range(2): + if sum_idx == 1: + print(crossing) + curves, w_val, patches = self._main_threads(surf1, sum_idx) + curve_id = -1 + for curve in curves: + curve_id += 1 + interval_id = -1 + for it in range(curve.basis.n_intervals): + interval_id += 1 + if self._already_found(crossing,interval_id,curve_id,sum_idx) == 1: + #print('continue') + #print(crossing) + #print(point.R3_coor) + continue + curv_surf_isec = IsecCurveSurf(surf2, curve) + boxes = bih.AABB(curve.poles[it:it+3, :].tolist()) + uv1 = np.zeros([2,1]) + if np.logical_or(it == 0, it == curve.basis.n_intervals - 1): + interval_list = [it] + else: + interval_list = [it, it] + uv1[sum_idx,0] = w_val[curve_id] + intersectioned_patches2 = tree2.find_box(boxes) + #print("curve_id=", curve_id) + #print("sum_idx=", sum_idx) + #print("intersectioned_patches2=",intersectioned_patches2) + for ipatch2 in intersectioned_patches2: + iu2, iv2 = self._patch_id2pos(surf2,ipatch2) + uvt = curv_surf_isec.get_initial_condition(iu2, iv2, it) + uvt, conv, flag, XYZ = curv_surf_isec.get_intersection(uvt, iu2, iv2, it, self.max_it, + self.rel_tol, self.abs_tol) + if conv == 1: + #check second surf + #else break + if np.logical_or(np.logical_and(flag[2] == 0, it == 0), + np.logical_and(flag[2] == 1, it == curve.basis.n_intervals - 1)): + boundary_flag = 1 + else: + boundary_flag = 0 + + uv1[1 - sum_idx, 0] = uvt[2, 0] + if sum_idx == 0: + point = IsecPoint(surf1, patches[curve_id], interval_list, uv1, boundary_flag, flag, + sum_idx, surf2, [iu2], [iv2], uvt[0:2, :], XYZ) + elif sum_idx == 1: + point = IsecPoint(surf1, interval_list, patches[curve_id], uv1, boundary_flag, flag, + sum_idx,surf2, [iu2], [iv2], uvt[0:2, :], XYZ) + + point_list.append(point) + + #if boundary_flag == 0: + # point_list[- 1].add_patches(sum_idx, flag[2]) # move to constructor + #print('flag') + #print(flag.shape) + #print(flag[2]) + if np.logical_or(flag[2] == 0,flag[2] == 1): + if sum_idx == 0: + crossing[curve_id, interval_id + flag[2]] = 1 + elif sum_idx == 1: + crossing[interval_id + flag[2],curve_id] = 1 # xxx + break + + + + return point_list + + @staticmethod + def _already_found(crossing,interval_id,curve_id,sum_idx): + + found = 0 + + + if sum_idx == 0: + for i in range(2): + if crossing[interval_id + i,curve_id] == 1: + found = 1 + elif sum_idx == 1: + for i in range(2): + if crossing[curve_id, interval_id + i] == 1: + found = 1 + + return found + + + +# np.logical_not(np.logical_and(crossing[interval_id+flag,curve_id] == 1,sum_idx == 1)) + ''' + Calculation and representation of intersection of two B-spline surfaces. + + Result is set of B-spline segments approximating the intersection branches. + Every segment consists of: a 3D curve and two 2D curves in parametric UV space of the surfaces. + ''' diff --git a/swaplines.m b/swaplines.m deleted file mode 100644 index b9ad811..0000000 --- a/swaplines.m +++ /dev/null @@ -1,8 +0,0 @@ -function [ mat ] = swaplines( mat,i,j ) - -temp = mat(j,:); -mat(j,:) = mat(i,:); -mat(i,:)= temp; - -end - diff --git a/tests/test_brep_writer.py b/tests/test_brep_writer.py new file mode 100644 index 0000000..598435f --- /dev/null +++ b/tests/test_brep_writer.py @@ -0,0 +1,157 @@ +import pytest +import unittest +import brep_writer as bw +import sys + + + +class TestConstructors(unittest.TestCase): + def test_Location(self): + print( "test locations") + la = bw.Location([[1, 0, 0, 4], [0, 1, 0, 8], [0, 0, 1, 12]]) + lb = bw.Location([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]) + lc = bw.ComposedLocation([ (la, 2), (lb, 1) ]) + ld = bw.ComposedLocation([ (lb, 2), (lc, 3) ]) + with self.assertRaises(bw.ParamError): + bw.Location([1,2,3]) + with self.assertRaises(bw.ParamError): + bw.Location([[1], [2], [3]]) + with self.assertRaises(bw.ParamError): + a = 1 + b = 'a' + lb = bw.Location([[a, b, a, b], [a, b, a, b], [a, b, a, b]]) + + def test_Shape(self): + + # check sub types + with self.assertRaises(bw.ParamError): + bw.Wire(['a', 'b']) + v=bw.Vertex([1,2,3]) + with self.assertRaises(bw.ParamError): + bw.Wire([v, v]) + + + def test_Vertex(self): + with self.assertRaises(bw.ParamError): + bw.Vertex(['a','b','c']) + + +class TestPlanarGeomeries(unittest.TestCase): + + def _cube(self): + # 0, 0; top, bottom + v1=bw.Vertex((0.0, 0.0, 1.0)) + v2=bw.Vertex((0.0, 0.0, 0.0)) + + v3=bw.Vertex((0.0, 1.0, 1.0)) + v4=bw.Vertex((0.0, 1.0, 0.0)) + + v5=bw.Vertex((1.0, 0.0, 1.0)) + v6=bw.Vertex((1.0, 0.0, 0.0)) + + # vertical edges + e1=bw.Edge([v1,v2]) + e2=bw.Edge([v3,v4]) + e3=bw.Edge([v5,v6]) + + # face v12 - v34 + # top + e4=bw.Edge([v1,v3]) + # bottom + e5=bw.Edge([v2,v4]) + f1 = bw.Face([e1.m(), e4, e2, e5.m()]) + + # face v34 - v56 + # top + e6=bw.Edge([v3, v5]) + # bottom + e7=bw.Edge([v4, v6]) + f2 = bw.Face([e2.m(), e6, e3, e7.m()]) + + # face v56 - v12 + # top + e8=bw.Edge([v5, v1]) + # bottom + e9=bw.Edge([v6, v2]) + f3 = bw.Face([e3.m(), e8, e1, e9.m()]) + + # top cup + f4 = bw.Face([e4, e6, e8]) + # bot cup + w5=bw.Wire([e5, e7, e9]) + f5 = bw.Face([w5.m()]) + + shell = bw.Shell([ f1, f2, f3, f4, f5.m() ]) + return shell + + def _permuted_cube(self): + # 0, 0; top, bottom + v1=bw.Vertex((0.0, 0.0, 1.0)) + v2=bw.Vertex((0.0, 0.0, 0.0)) + + v3=bw.Vertex((0.0, 1.0, 1.0)) + v4=bw.Vertex((0.0, 1.0, 0.0)) + + v5=bw.Vertex((1.0, 0.0, 1.0)) + v6=bw.Vertex((1.0, 0.0, 0.0)) + + # face v12 - v34 + # top + e4=bw.Edge([v1,v3]) + # top + e6=bw.Edge([v3, v5]) + # top + e8=bw.Edge([v5, v1]) + + # top cup + f4 = bw.Face([e4, e6, e8]) + + # bottom + e5=bw.Edge([v2,v4]) + # face v34 - v56 + # bottom + e7=bw.Edge([v4, v6]) + # face v56 - v12 + # bottom + e9=bw.Edge([v6, v2]) + + # bot cup + w5=bw.Wire([e5, e7, e9]) + f5 = bw.Face([w5.m()]) + + + # vertical edges + e1=bw.Edge([v1,v2]) + e2=bw.Edge([v3,v4]) + e3=bw.Edge([v5,v6]) + + f1 = bw.Face([e1.m(), e4, e2, e5.m()]) + + f2 = bw.Face([e2.m(), e6, e3, e7.m()]) + + f3 = bw.Face([e3.m(), e8, e1, e9.m()]) + + shell = bw.Shell([ f4, f5.m(), f1, f2, f3 ]) + return shell + + def test_cube(self): + + shell = self._cube() + #shell = self._permuted_cube() + + s1=bw.Solid([ shell ]) + + c1=bw.Compound([s1]) + + loc1=bw.Location([[0,0,1,0],[1,0,0,0],[0,1,0,0]]) + loc2=bw.Location([[0,0,1,0],[1,0,0,0],[0,1,0,0]]) #dej tam tu druhou z prikladu + cloc=bw.ComposedLocation([(loc1,1),(loc2,1)]) + + with open("_out_test_prism.brep", "w") as f: + bw.write_model(f, c1, cloc) + #bw.write_model(sys.stdout, c1, cloc) + print(c1) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/test_bs_approx.py b/tests/test_bs_approx.py new file mode 100644 index 0000000..655749a --- /dev/null +++ b/tests/test_bs_approx.py @@ -0,0 +1,262 @@ +import pytest +import numpy as np +import bspline as bs +import bspline_approx as bs_approx +import math +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +import bspline_plot as bs_plot +import time +import logging + + +def function_sin_cos(x): + return math.sin(x[0] * 4) * math.cos(x[1] * 4) + + +def gen_uv_grid(nu, nv): + # surface on unit square + U = np.linspace(0.0, 1.0, nu) + V = np.linspace(0.0, 1.0, nv) + V_grid, U_grid = np.meshgrid(V, U) + + return np.stack([U_grid.ravel(), V_grid.ravel()], axis=1) + +def eval_func_on_grid(func, xy_mat, xy_shift, z_mat, shape=(50, 50)): + nu, nv = shape + UV = gen_uv_grid(nu, nv) + XY = xy_mat.dot(UV.T).T + xy_shift + z_func_eval = np.array([z_mat[0] * func([u, v]) + z_mat[1] for u, v in UV]) + return np.concatenate((XY, z_func_eval[:, None]), axis=1).reshape(nu, nv, 3) + + +def eval_z_surface_on_grid(surface, xy_mat, xy_shift, shape=(50, 50)): + nu, nv = shape + UV = gen_uv_grid(nu, nv) + XY = xy_mat.dot(UV.T).T + xy_shift + Z = surface.z_eval_xy_array(XY) + return np.concatenate((XY, Z[:, None]), axis=1).reshape(nu, nv, 3) + + +def eval_surface_on_grid(surface, shape=(50, 50)): + nu, nv = shape + UV = gen_uv_grid(nu, nv) + return surface.eval_array(UV).reshape(nu, nv, 3) + +def plot_cmp(a_grid, b_grid): + plt = bs_plot.Plotting() + plt.plot_surface_3d(a_grid[:, :, 0], a_grid[:, :, 1], a_grid[:, :, 2]) + plt.plot_surface_3d(b_grid[:, :, 0], b_grid[:, :, 1], b_grid[:, :, 2]) + plt.show() + + diff = b_grid - a_grid + plt.plot_surface_3d(a_grid[:, :, 0], a_grid[:, :, 1], diff[:, :, 0]) + plt.plot_surface_3d(a_grid[:, :, 0], a_grid[:, :, 1], diff[:, :, 1]) + plt.plot_surface_3d(a_grid[:, :, 0], a_grid[:, :, 1], diff[:, :, 2]) + plt.show() + +def grid_cmp( a, b, tol): + a_z = a[:, :, 2].ravel() + b_z = b[:, :, 2].ravel() + eps = 0.0 + n_err=0 + for i, (za, zb) in enumerate(zip(a_z, b_z)): + diff = np.abs( za - zb) + eps = max(eps, diff) + if diff > tol: + n_err +=1 + if n_err < 10: + print(" {} =|a({}) - b({})| > tol({}), idx: {}".format(diff, za, zb, tol, i) ) + elif n_err == 10: + print("... skipping") + print("Max norm: ", eps, "Tol: ", tol) + if eps > tol: + plot_cmp(a, b) + +class TestSurfaceApprox: + # todo: numerical test + + + def test_approx_func(self): + logging.basicConfig(level=logging.DEBUG) + + xy_mat = np.array( [ [1.0, -1.0, 10 ], [1.0, 1.0, 20 ]]) # rotate left pi/4 and blow up 1.44 + z_mat = np.array( [2.0, -10] ) + xyz_func = eval_func_on_grid(function_sin_cos, xy_mat[:2,:2], xy_mat[:, 2], z_mat) + + # print("Compare: Func - GridSurface.transform") + # points = bs.make_function_grid(function_sin_cos, 200, 200) + # gs = bs.GridSurface(points.reshape(-1, 3)) + # gs.transform(xy_mat, z_mat) + # xyz_grid = eval_z_surface_on_grid(gs, xy_mat[:2,:2], xy_mat[:, 2]) + # grid_cmp(xyz_func, xyz_grid, 0.02) + # + # print("Compare: Func - GridSurface.transform.z_surf") + # xyz_grid = eval_z_surface_on_grid(gs.z_surface, xy_mat[:2, :2], xy_mat[:, 2]) + # xyz_grid[:,:,2] *= z_mat[0] + # xyz_grid[:, :, 2] += z_mat[1] + # grid_cmp(xyz_func, xyz_grid, 0.02) + + print("\nCompare: Func - GridSurface.approx") + points = bs.make_function_grid(function_sin_cos, 50, 50) + gs = bs.GridSurface(points.reshape(-1, 3)) + gs.transform(xy_mat, z_mat) + approx = bs_approx.SurfaceApprox.approx_from_grid_surface(gs) + surface = approx.compute_approximation() + xyz_grid = eval_z_surface_on_grid(surface, xy_mat[:2, :2], xy_mat[:, 2]) + print("Approx error: ", approx.error) + grid_cmp(xyz_func, xyz_grid, 0.02) + + print("\nCompare: Func - points.approx") + np.random.seed(seed=123) + uv = np.random.rand(1000,2) + xy = xy_mat[:2,:2].dot(uv.T).T + xy_mat[:, 2] + z = np.array( [function_sin_cos([u, v]) for u, v in uv] ) + xyz = np.concatenate((xy, z[:, None]), axis=1) + approx = bs_approx.SurfaceApprox(xyz) + quad = approx.compute_default_quad() + + nuv = approx.nuv + ref_quad = np.array([ [-1 , 1], [0,0], [1, 1], [0, 2] ]) + ref_quad += np.array([10, 20]) + assert np.allclose(ref_quad, quad, atol=1e-2) + + nuv = approx.compute_default_nuv() + assert np.allclose( np.array([8, 8]), nuv) + + surface = approx.compute_approximation() + surface.transform(xy_mat = None, z_mat=z_mat) + nu, nv = 50, 50 + uv_probe = gen_uv_grid(nu, nv) + uv_probe = (0.9*uv_probe + 0.05) + xy_probe = xy_mat[:2,:2].dot(uv_probe.T).T + xy_mat[:, 2] + z_func = np.array( [z_mat[0] * function_sin_cos([u, v]) + z_mat[1] for u, v in uv_probe] ) + xyz_func = np.concatenate((xy_probe, z_func[:, None]), axis=1).reshape(nu, nv, 3) + xyz_approx = surface.eval_xy_array(xy_probe).reshape(nu, nv, 3) + print("Approx error: ", approx.error) + grid_cmp(xyz_func, xyz_approx, 0.02) + + # approx = bs_approx.SurfaceApprox(xyz) + # surface = approx.compute_approximation() + # + # + # plt = bs_plot.Plotting() + # plt.plot_surface_3d(gs) + # plt.plot_surface_3d(surface) + # plt.scatter_3d(xyz[:,0], xyz[:,1], xyz[:,2]) + # plt.show() + + # xyz_grid = eval_z_surface_on_grid(z_surf, xy_mat[:2,:2], xy_mat[:, 2]) + # grid_cmp(xyz_func, xyz_grid, 0.02) + # + # print("Compare: Func - GridSurface.transform.Z_surf_approx.full_surface") + # xyz_grid = eval_surface_on_grid(z_surf.make_full_surface()) + # grid_cmp(xyz_func, xyz_grid, 0.01) + + def test_approx_real_surface(self): + # Test approximation of real surface grid using a random subset + # hard test of regularization. + logging.basicConfig(level=logging.DEBUG) + + xy_mat = np.array( [ [1.0, 0.0, 0 ], [0.0, 1.0, 0 ]]) # rotate left pi/4 and blow up 1.44 + z_mat = np.array( [1.0, 0] ) + xyz_func = eval_func_on_grid(function_sin_cos, xy_mat[:2,:2], xy_mat[:, 2], z_mat) + + print("Compare: Func - Randomized.approx") + points = bs.make_function_grid(function_sin_cos, 50, 50) + points = points.reshape( (-1, 3) ) + n_sample_points= 400 # this is near the limit number of points to keep desired precision + random_subset = np.random.random_integers(0, len(points)-1, n_sample_points) + points_random = points[random_subset, :] + + approx = bs_approx.SurfaceApprox(points_random) + approx.set_quad(None) # set unit square + surface = approx.compute_approximation() + xyz_grid = eval_z_surface_on_grid(surface, xy_mat[:2, :2], xy_mat[:, 2]) + print("Approx error: ", approx.error) + grid_cmp(xyz_func, xyz_grid, 0.02) + + + def test_transformed_quad(self): + xy_mat = np.array( [ [1.0, -1.0, 0 ], [1.0, 1.0, 0 ]]) # rotate left pi/4 and blow up 1.44 + np.random.seed(seed=123) + uv = np.random.rand(1000,2) + xy = xy_mat[:2,:2].dot(uv.T).T + xy_mat[:, 2] + z = np.array( [function_sin_cos([u, v]) for u, v in uv] ) + xyz = np.concatenate((xy, z[:, None]), axis=1) + approx = bs_approx.SurfaceApprox(xyz) + approx.quad = np.array([ [-1 , 1], [0,0], [1, 1], [0, 2] ]) + + xy_mat = np.array([[2, 1, -1],[1, 2, -2]]) + assert np.allclose( np.array([ [-2 , -1], [-1,-2], [2, 1], [1, 2] ]), approx.transformed_quad(xy_mat)) + + + def plot_approx_transformed_grid(self): + pass + + def plot_plane(self): + surf = bs_approx.plane_surface([ [0.0, 0, 0], [1.0, 0, 0], [0.0, 0, 1] ], overhang=0.1) + self.plot_surf(surf) + + def test_approx_speed(self): + logging.basicConfig(level=logging.DEBUG) + + print("Performance test for 100k points.") + np.random.seed(seed=123) + uv = np.random.rand(100000,2) + xy = uv + z = np.array( [function_sin_cos([u, v]) for u, v in uv] ) + xyz = np.concatenate((xy, z[:, None]), axis=1) + start_time = time.time() + approx = bs_approx.SurfaceApprox(xyz) + surface = approx.compute_approximation() + end_time = time.time() + print("\nApprox 100k points by 100x100 grid in: {} sec".format(end_time - start_time)) + assert end_time - start_time < 6 + # target is approximation of 1M points in one minute + # B matrix 3.6 sec, A matrix 0.7 sec, SVD + Z solve 0.7 sec + +class TestBoundingBox: + def test_hull_and_box(self): + points = np.random.randn(1000000,2) + + print() + start = time.perf_counter() + for i in range(1): + hull = bs_approx.convex_hull_2d(points) + end = time.perf_counter() + print("\nConvex hull of 1M points: {} s".format(end - start)) + + start = time.perf_counter() + for i in range(10): + quad = bs_approx.min_bounding_rect(hull) + end = time.perf_counter() + print("Min area bounding box: {} s".format(end - start)) + return + + plt = bs_plot.Plotting() + plt.scatter_2d(points[:,0], points[:,1]) + plt.plot_2d(hull[:, 0], hull[:, 1]) + box_lines = np.concatenate((quad, quad[0:1,:]), axis=0) + plt.plot_2d(box_lines[:, 0], box_lines[:, 1]) + plt.show() + +class TestCurveApprox: + + def plot_approx_2d(self): + x_vec = np.linspace(1.1, 3.0, 100) + y_vec = np.array([ np.sin(10*x) for x in x_vec ]) + points = np.stack( (x_vec, y_vec), axis=1) + curve = bs_approx.curve_from_grid(points) + + bs_plot.plot_curve_2d(curve, 1000) + bs_plot.plot_curve_poles_2d(curve) + + plt.show() + + #plt.plot(x_vec, y_vec, color='green') + #plt.show() + + def test_approx_2d(self): + # self.plot_approx_2d() + pass \ No newline at end of file diff --git a/tests/test_bspline.py b/tests/test_bspline.py new file mode 100644 index 0000000..d35c50b --- /dev/null +++ b/tests/test_bspline.py @@ -0,0 +1,481 @@ +import pytest +import bspline as bs +import numpy as np +import math +import bspline_plot as bp + +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.pyplot as plt + +plotting = bp.Plotting() +#plotting = bp.Plotting(bp.PlottingMatplot()) + + +class TestSplineBasis: + + def test_find_knot_interval(self): + """ + test methods: + - make_equidistant + - find_knot_interval + """ + eq_basis = bs.SplineBasis.make_equidistant(2, 100) + assert eq_basis.find_knot_interval(0.0) == 0 + assert eq_basis.find_knot_interval(0.001) == 0 + assert eq_basis.find_knot_interval(0.01) == 1 + assert eq_basis.find_knot_interval(0.011) == 1 + assert eq_basis.find_knot_interval(0.5001) == 50 + assert eq_basis.find_knot_interval(1.0 - 0.011) == 98 + assert eq_basis.find_knot_interval(1.0 - 0.01) == 99 + assert eq_basis.find_knot_interval(1.0 - 0.001) == 99 + assert eq_basis.find_knot_interval(1.0) == 99 + + knots = np.array([0, 0, 0, 0.1880192, 0.24545785, 0.51219762, 0.82239001, 1., 1. , 1.]) + basis = bs.SplineBasis(2, knots) + for interval in range(2, 7): + xx = np.linspace(knots[interval], knots[interval+1], 10) + for j, x in enumerate(xx[:-1]): + i_found = basis.find_knot_interval(x) + assert i_found == interval - 2, "i_found: {} i: {} j: {} x: {} ".format(i_found, interval-2, j, x) + + + def test_packed_knots(self): + """ + Test: + - make_from_packed_knots + - pack_knots + :return: + """ + packed = [(-0.1, 3), (0,1), (1,1), (1.1, 3)] + basis = bs.SplineBasis.make_from_packed_knots(2, packed) + assert packed == basis.pack_knots() + + + + + def plot_basis(self, eq_basis): + n_points = 401 + x_coord = np.linspace(eq_basis.domain[0], eq_basis.domain[1], n_points) + + for i_base in range(eq_basis.size): + y_coord = [ eq_basis.eval(i_base, x) for x in x_coord ] + plotting.plot_2d(x_coord, y_coord) + + plotting.show() + + + def test_eval(self): + #self.plot_basis(bs.SplineBasis.make_equidistant(0, 4)) + + knots = np.array([0, 0, 0, 0.1880192, 0.24545785, 0.51219762, 0.82239001, 1., 1. , 1.]) + basis = bs.SplineBasis(2, knots) + #self.plot_basis(basis) + + eq_basis = bs.SplineBasis.make_equidistant(0, 2) + assert eq_basis.eval(0, 0.0) == 1.0 + assert eq_basis.eval(1, 0.0) == 0.0 + assert eq_basis.eval(0, 0.5) == 0.0 + assert eq_basis.eval(1, 0.5) == 1.0 + assert eq_basis.eval(1, 1.0) == 1.0 + + eq_basis = bs.SplineBasis.make_equidistant(1, 4) + assert eq_basis.eval(0, 0.0) == 1.0 + assert eq_basis.eval(1, 0.0) == 0.0 + assert eq_basis.eval(2, 0.0) == 0.0 + assert eq_basis.eval(3, 0.0) == 0.0 + assert eq_basis.eval(4, 0.0) == 0.0 + + assert eq_basis.eval(0, 0.125) == 0.5 + assert eq_basis.eval(1, 0.125) == 0.5 + assert eq_basis.eval(2, 1.0) == 0.0 + + # check summation to one: + for deg in range(0, 10): + basis = bs.SplineBasis.make_equidistant(deg, 2) + for x in np.linspace(basis.domain[0], basis.domain[1], 10): + s = sum([ basis.eval(i, x) for i in range(basis.size) ]) + assert np.isclose(s, 1.0) + + def fn_supp(self): + basis = bs.SplineBasis.make_equidistant(2, 4) + for i in range(basis.size): + supp = basis.fn_supp(i) + for x in np.linspace(supp[0] - 0.1, supp[0], 10): + assert basis.eval(i, x) == 0.0 + for x in np.linspace(supp[0] + 0.001, supp[1] - 0.001, 10): + assert basis.eval(i, x) > 0.0 + for x in np.linspace(supp[1], supp[1] + 0.1): + assert basis.eval(i, x) == 0.0 + + def test_linear_poles(self): + eq_basis = bs.SplineBasis.make_equidistant(2, 4) + poles = eq_basis.make_linear_poles() + + t_vec = np.linspace(0.0, 1.0, 21) + for t in t_vec: + b_vals = np.array([ eq_basis.eval(i, t) for i in range(eq_basis.size) ]) + x = np.dot(b_vals, poles) + assert np.abs( x - t ) < 1e-15 + + def check_eval_vec(self, basis, i, t): + vec = basis.eval_vector(i, t) + for j in range(basis.degree + 1): + assert vec[j] == basis.eval(i + j, t) + + def plot_basis_vec(self, basis): + n_points = 401 + x_coord = np.linspace(basis.domain[0], basis.domain[1], n_points) + + y_coords = np.zeros( (basis.size, x_coord.shape[0]) ) + for i, x in enumerate(x_coord): + idx = basis.find_knot_interval(x) + y_coords[idx : idx + basis.degree + 1, i] = basis.eval_vector(idx, x) + + for i_base in range(basis.size): + plotting.plot_2d(x_coord, y_coords[i_base, :]) + + plotting.show() + + + def test_eval_vec(self): + basis = bs.SplineBasis.make_equidistant(2, 4) + #self.plot_basis_vec(basis) + self.check_eval_vec(basis, 0, 0.1) + self.check_eval_vec(basis, 1, 0.3) + self.check_eval_vec(basis, 2, 0.6) + self.check_eval_vec(basis, 3, 0.8) + self.check_eval_vec(basis, 3, 1.0) + + basis = bs.SplineBasis.make_equidistant(3, 4) + # self.plot_basis_vec(basis) + self.check_eval_vec(basis, 0, 0.1) + self.check_eval_vec(basis, 1, 0.3) + self.check_eval_vec(basis, 2, 0.6) + self.check_eval_vec(basis, 3, 0.8) + self.check_eval_vec(basis, 3, 1.0) + + + + def check_diff_vec(self, basis, i, t): + vec = basis.eval_diff_vector(i, t) + for j in range(basis.degree + 1): + assert np.abs(vec[j] - basis.eval_diff(i + j, t)) < 1e-15 + + def plot_basis_diff(self, basis): + n_points = 401 + x_coord = np.linspace(basis.domain[0], basis.domain[1], n_points) + + y_coords = np.zeros( (basis.size, x_coord.shape[0]) ) + for i, x in enumerate(x_coord): + idx = basis.find_knot_interval(x) + y_coords[idx : idx + basis.degree + 1, i] = basis.eval_diff_vector(idx, x) + + for i_base in range(basis.size): + plotting.plot_2d(x_coord, y_coords[i_base, :]) + + plotting.show() + + + def test_eval_diff_base_vec(self): + basis = bs.SplineBasis.make_equidistant(2, 4) + #self.plot_basis_diff(basis) + self.check_diff_vec(basis, 0, 0.1) + self.check_diff_vec(basis, 1, 0.3) + self.check_diff_vec(basis, 2, 0.6) + self.check_diff_vec(basis, 3, 0.8) + self.check_diff_vec(basis, 3, 1.0) + + basis = bs.SplineBasis.make_equidistant(3, 4) + # self.plot_basis_diff(basis) + self.check_diff_vec(basis, 0, 0.1) + self.check_diff_vec(basis, 1, 0.3) + self.check_diff_vec(basis, 2, 0.6) + self.check_diff_vec(basis, 3, 0.8) + self.check_diff_vec(basis, 3, 1.0) + + + +class TestCurve: + + def plot_4p(self): + degree = 2 + poles = [ [0., 0.], [1.0, 0.5], [2., -2.], [3., 1.] ] + basis = bs.SplineBasis.make_equidistant(degree, 2) + curve = bs.Curve(basis, poles) + + plotting.plot_curve_2d(curve, poles=True) + b00, b11 = curve.aabb() + b01 = [b00[0], b11[1]] + b10 = [b11[0], b00[1]] + bb = np.array([b00, b10, b11, b01, b00]) + + plotting.plot_2d( bb[:, 0], bb[:, 1]) + plotting.show() + + def test_evaluate(self): + #self.plot_4p() + # TODO: make numerical tests with explicitely computed values + # TODO: test rational curves, e.g. circle + + pass + + def test_aabb(self): + poles = [ [0., 0.], [1.0, 0.5], [2., -2.], [3., 1.] ] + basis = bs.SplineBasis.make_equidistant(2, 2) + curve = bs.Curve(basis, poles) + box = curve.aabb() + assert np.allclose( box, np.array([ [0,-2], [3, 1]]) ) + + + + + + + + + +class TestSurface: + + def plot_extrude(self): + + # curve extruded to surface + poles_yz = [[0., 0.], [1.0, 0.5], [2., -2.], [3., 1.]] + poles_x = [0, 1, 2] + poles = [ [ [x] + yz for yz in poles_yz ] for x in poles_x ] + u_basis = bs.SplineBasis.make_equidistant(2, 1) + v_basis = bs.SplineBasis.make_equidistant(2, 2) + surface_extrude = bs.Surface( (u_basis, v_basis), poles) + plotting.plot_surface_3d(surface_extrude, poles = True) + plotting.show() + + def plot_function(self): + fig = plt.figure() + ax = fig.gca(projection='3d') + + # function surface + def function(x): + return math.sin(x[0]) * math.cos(x[1]) + + poles = bs.make_function_grid(function, 4, 5) + u_basis = bs.SplineBasis.make_equidistant(2, 2) + v_basis = bs.SplineBasis.make_equidistant(2, 3) + surface_func = bs.Surface( (u_basis, v_basis), poles) + plotting.plot_surface_3d(surface_func, poles = True) + plotting.show() + + def test_evaluate(self): + #self.plot_extrude() + #self.plot_function() + # TODO: test rational surfaces, e.g. sphere + pass + + + def test_aabb(self): + # function surface + def function(x): + return x[0] * (x[1] + 1.0) + 3.0 + + poles = bs.make_function_grid(function, 4, 5) + u_basis = bs.SplineBasis.make_equidistant(2, 2) + v_basis = bs.SplineBasis.make_equidistant(2, 3) + surface_func = bs.Surface( (u_basis, v_basis), np.array(poles)) + box = surface_func.aabb() + assert np.allclose( box, np.array([ [0,0, 3], [1, 1, 5]]) ) + + + +class TestZ_Surface: + + # TODO: Compute max norm of the difference of two surfaces and assert that it is cose to zero. + + def make_z_surf(self, func, quad): + poles = bs.make_function_grid(func, 4, 5) + u_basis = bs.SplineBasis.make_equidistant(2, 2) + v_basis = bs.SplineBasis.make_equidistant(2, 3) + surface_func = bs.Surface( (u_basis, v_basis), poles[:,:, [2] ]) + return bs.Z_Surface(quad, surface_func) + + def plot_function_uv(self): + # function surface + def function(x): + return math.sin(x[0]*4) * math.cos(x[1]*4) + + quad = np.array([[0, 0], [1, 0], [1, 1], [0, 1]]) + z_surf = self.make_z_surf(function, quad) + + z_surf.transform(None, np.array([2.0, 0]) ) + full_surf = z_surf.make_full_surface() + z_surf.transform(None, np.array([1.0, 0.1]) ) + plotting.plot_surface_3d(z_surf) + plotting.plot_surface_3d(full_surf) + + plotting.show() + + def test_eval_uv(self): + #self.plot_function_uv() + pass + + def test_aabb(self): + # function surface + def function(x): + return x[0] * (x[1] + 1.0) + 3.0 + + quad = np.array( [ [0, 0], [0, 0.5], [1, 0.1], [1.1, 1.1] ] ) + z_surf = self.make_z_surf(function, quad) + + box = z_surf.aabb() + assert np.allclose(box, np.array([[0, 0, 3], [1.1, 1.1, 5]])) + + z_surf.transform(None, np.array([2.0, 1.0])) + box = z_surf.aabb() + assert np.allclose(box, np.array([[0, 0, 7], [1.1, 1.1, 11]])) + + def test_reset_transform(self): + def function(x): + return math.sin(x[0]*4) * math.cos(x[1]*4) + + quad = np.array([[1., 3.5], [1., 2.], [2., 2.2], [2, 3.7]]) + z_surf = self.make_z_surf(function, quad) + xy_mat = np.array([[2, 1, -1],[1, 2, -2]]) + z_mat = np.array([2, 1]) + z_surf.transform(xy_mat, z_mat) + assert np.allclose( z_surf.quad[0], np.array([4.5, 6]) ) + z_surf.reset_transform() + assert np.allclose(z_surf.quad, quad) + + def test_get_copy(self): + def function(x): + return math.sin(x[0]*4) * math.cos(x[1]*4) + + quad = np.array([[1., 3.5], [1., 2.], [2., 2.2], [2, 3.7]]) + a_surf = self.make_z_surf(function, quad) + b_surf = a_surf.get_copy() + b_surf.transform(None, np.array([2, 1])) + assert np.allclose(a_surf.center()[0:1], b_surf.center()[0:1]) + assert np.isclose(2 * a_surf.center()[2] + 1, b_surf.center()[2]) + za = a_surf.center()[2] + zb = b_surf.center()[2] + + b_surf.apply_z_transform() + assert np.allclose(a_surf.center()[0:1], b_surf.center()[0:1]) + assert np.isclose(zb, b_surf.center()[2]) + assert np.isclose(za, a_surf.center()[2]) + + + +class TestPointGrid: + + @staticmethod + def function(x): + return math.sin(x[0]) * math.cos(x[1]) + + + def make_point_grid(self): + nu, nv = 5,6 + grid = bs.make_function_grid(TestPointGrid.function, 5, 6).reshape(nu*nv, 3) + surf = bs.GridSurface(grid) + return surf + + + def plot_check_surface(self, XYZ_grid_eval, XYZ_surf_eval, XYZ_func_eval): + plotting.plot_surface(XYZ_grid_eval[:, :, 0], XYZ_grid_eval[:, :, 1], XYZ_grid_eval[:, :, 2], color='blue') + plotting.plot_surface(XYZ_surf_eval[:, :, 0], XYZ_surf_eval[:, :, 1], XYZ_surf_eval[:, :, 2], color='red') + plotting.plot_surface(XYZ_func_eval[:, :, 0], XYZ_func_eval[:, :, 1], XYZ_func_eval[:, :, 2], color='green') + plotting.show() + + def grid_cmp(self, a, b, tol): + a_z = a[:, :, 2].ravel() + b_z = b[:, :, 2].ravel() + eps = 0.0 + for i, (za, zb) in enumerate(zip(a_z, b_z)): + diff = np.abs( za - zb) + eps = max(eps, diff) + assert diff < tol, " |a({}) - b({})| > tol({}), idx: {}".format(za, zb, tol, i) + print("Max norm: ", eps, "Tol: ", tol) + + def check_surface(self, surf, xy_mat, xy_shift, z_mat): + """ + TODO: Make this a general function - evaluate a surface on a grid, use it also in other tests + to compare evaluation on the grid to the original function. Can be done after we have approximations. + """ + + nu, nv = 30, 40 + # surface on unit square + U = np.linspace(0.0, 1.0, nu) + V = np.linspace(0.0, 1.0, nv) + V_grid, U_grid = np.meshgrid(V,U) + + UV = np.stack( [U_grid.ravel(), V_grid.ravel()], axis = 1 ) + XY = xy_mat.dot(UV.T).T + xy_shift + Z = surf.z_eval_xy_array(XY) + XYZ_grid_eval = np.concatenate( (XY, Z[:, None]) , axis = 1).reshape(nu, nv, 3) + + XYZ_surf_eval = surf.eval_array(UV).reshape(nu, nv, 3) + + z_func_eval = np.array([ z_mat[0]*TestPointGrid.function([u,v]) + z_mat[1] for u, v in UV ]) + XYZ_func_eval = np.concatenate( (XY, z_func_eval[:, None]), axis =1 ).reshape(nu, nv, 3) + + #self.plot_check_surface(XYZ_grid_eval, XYZ_surf_eval, XYZ_func_eval) + + eps = 0.0 + hx = 1.0 / surf.shape[0] + hy = 1.0 / surf.shape[1] + tol = 0.5* ( hx*hx + 2*hx*hy + hy*hy) + + self.grid_cmp(XYZ_func_eval, XYZ_grid_eval, tol) + self.grid_cmp(XYZ_func_eval, XYZ_surf_eval, tol) + + + def test_grid_surface(self): + xy_mat = np.array([ [1.0, 0.0], [0.0, 1.0] ]) + xy_shift = np.array([0.0, 0.0 ]) + z_shift = np.array([1.0, 0.0]) + surface = self.make_point_grid() + + # fig = plt.figure() + # ax = fig.gca(projection='3d') + # bs_plot.plot_grid_surface_3d(surface, ax) + # plt.show() + # self.check_surface(surface, xy_mat, xy_shift, z_shift) + + # transformed surface + xy_mat = np.array([ [3.0, -3.0], [2.0, 2.0] ]) / math.sqrt(2) + xy_shift = np.array([[-2.0, 5.0 ]]) + z_shift = np.array([1.0, 1.3]) + new_quad = np.array([ [0, 1.0], [0,0], [1, 0], [1, 1]]) + new_quad = new_quad.dot(xy_mat[:2,:2].T) + xy_shift + + surface = self.make_point_grid() + surface.transform(np.concatenate((xy_mat, xy_shift.T), axis=1), z_shift) + assert np.all(surface.quad == new_quad), "surf: {} ref: {}".format(surface.quad, new_quad) + self.check_surface(surface, xy_mat, xy_shift, z_shift) + + surf_center = surface.center() + ref_center = np.array( [-2.0, math.sqrt(2.0)+5.0, (math.sin(0.5)*math.cos(0.5) + 1.3) ] ) + #print(surf_center) + #print(ref_center) + assert np.allclose( ref_center, surf_center, rtol=0.01) + + v_min, v_max = surface.aabb() + assert np.allclose(v_min, np.array([-3.0/math.sqrt(2) -2, 0.0 + 5, 1.3])) + assert np.allclose(v_max, np.array([3.0 / math.sqrt(2) - 2, 4.0 / math.sqrt(2) + 5, math.sin(1.0) + 1.3])) + + def test_grid_surface_transform(self): + surface = self.make_point_grid() + xy_mat = np.array([ [3.0, 0.0], [0.0, 2.0] ]) + xy_shift = np.array([[-2.0, 5.0 ]]) + surface.transform(np.concatenate((xy_mat, xy_shift.T), axis=1), None) + quad = surface.quad + #print(quad) + assert quad[0][0] == -2 + assert quad[0][1] == 7 + assert quad[2][0] == 1 + assert quad[2][1] == 5 + + surface.reset_transform() + quad = surface.quad + assert quad[0][0] == 0 + assert quad[0][1] == 1 + assert quad[2][0] == 1 + assert quad[2][1] == 0 + diff --git a/tests/test_bspline_plot.py b/tests/test_bspline_plot.py new file mode 100644 index 0000000..98caf36 --- /dev/null +++ b/tests/test_bspline_plot.py @@ -0,0 +1,62 @@ +import pytest +import bspline as bs +import numpy as np +import math +import bspline_plot as bp + + +class TestPlottingMatplot: + # plotting = bp.Plotting(bp.PlottingMatplot()) + + def plotting_2d(self, plotting): + # simple 2d graphs + n_points = 201 + eq_basis = bs.SplineBasis.make_equidistant(2, 4) + x_coord = np.linspace(eq_basis.domain[0], eq_basis.domain[1], n_points) + + for i_base in range(eq_basis.size): + y_coord = [ eq_basis.eval(i_base, x) for x in x_coord ] + plotting.plot_2d(x_coord, y_coord) + plotting.show() + + # 2d curves + degree = 2 + poles = [ [0., 0.], [1.0, 0.5], [2., -2.], [3., 1.] ] + basis = bs.SplineBasis.make_equidistant(degree, 2) + curve = bs.Curve(basis, poles) + plotting.plot_curve_2d(curve, poles=True) + + poles = [[0., 0.], [-1.0, 0.5], [-2., -2.], [3., 1.]] + basis = bs.SplineBasis.make_equidistant(degree, 2) + curve = bs.Curve(basis, poles) + plotting.plot_curve_2d(curve, poles=True) + + plotting.show() + + def test_plot_2d(self): + self.plotting_2d(bp.Plotting(bp.PlottingMatplot())) + self.plotting_2d(bp.Plotting((bp.PlottingPlotly()))) + + def plotting_3d(self, plotting): + # plotting 3d surfaces + def function(x): + return math.sin(x[0]*4) * math.cos(x[1]*4) + + poles = bs.make_function_grid(function, 4, 5) + u_basis = bs.SplineBasis.make_equidistant(2, 2) + v_basis = bs.SplineBasis.make_equidistant(2, 3) + surface_func = bs.Surface( (u_basis, v_basis), poles[:,:, [2] ]) + + #quad = np.array( [ [0, 0], [0, 0.5], [1, 0.1], [1.1, 1.1] ] ) + quad = np.array([[0, 0], [1, 0], [1, 1], [0, 1]]) + z_surf = bs.Z_Surface(quad, surface_func) + full_surf = z_surf.make_full_surface() + z_surf.transform(np.array([[1., 0, 0], [0, 1, 0]]), np.array([2.0, 0]) ) + plotting.plot_surface_3d(z_surf) + plotting.plot_surface_3d(full_surf) + + plotting.show() + + def test_plot_3d(self): + self.plotting_3d(bp.Plotting(bp.PlottingMatplot())) + self.plotting_3d(bp.Plotting(bp.PlottingPlotly())) \ No newline at end of file diff --git a/tests/test_isec.py b/tests/test_isec.py new file mode 100644 index 0000000..bc5de20 --- /dev/null +++ b/tests/test_isec.py @@ -0,0 +1,114 @@ +import pytest +import isec_surf_surf as iss +import bspline as bs +import numpy as np +import math +import bspline_plot as bp + + +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.pyplot as plt + +class TestSurface: + + def plot_extrude(self): + #fig1 = plt.figure() + + #ax1 = fig1.gca(projection='3d') + + + + def function(x): + return math.sin(x[0]*4) * math.cos(x[1]*4) + + def function2(x): + return math.cos(x[0]*4) * math.sin(x[1]*4) + + u1_int = 4 + v1_int = 4 + u2_int = 4 + v2_int = 4 + + u_basis = bs.SplineBasis.make_equidistant(2, u1_int) #10 + v_basis = bs.SplineBasis.make_equidistant(2, v1_int) #15 + u2_basis = bs.SplineBasis.make_equidistant(2, u2_int) #10 + v2_basis = bs.SplineBasis.make_equidistant(2, v2_int) #15 + poles = bs.make_function_grid(function, u1_int + 2, v1_int + 2) #12, 17 + surface_extrude = bs.Surface((u_basis, v_basis), poles) + + myplot = bp.Plotting((bp.PlottingPlotly())) + myplot.plot_surface_3d(surface_extrude, poles = False) + poles2 = bs.make_function_grid(function2, u2_int + 2, v2_int + 2) #12, 17 + surface_extrude2 = bs.Surface((u2_basis, v2_basis), poles2) + myplot.plot_surface_3d(surface_extrude2, poles=False) + + + + return surface_extrude, surface_extrude2, myplot + + + + #def plot_function(self): + # fig = plt.figure() + # ax = fig.gca(projection='3d') + + # function surface + # def function(x): + # return math.sin(x[0]) * math.cos(x[1]) + + #poles = bs.make_function_grid(function, 4, 5) + #u_basis = bs.SplineBasis.make_equidistant(2, 2) + #v_basis = bs.SplineBasis.make_equidistant(2, 3) + #surface_func = bs.Surface( (u_basis, v_basis), poles) + #bs_plot.plot_surface_3d(surface_func, ax) + #bs_plot.plot_surface_poles_3d(surface_func, ax) + + #plt.show() + + + #def boudingbox(self): + # SurfSurf = IsecSurfSurf.bounding_boxes() + + def test_isec(self): + + + surf1, surf2, myplot = self.plot_extrude() + isec = iss.IsecSurfSurf(surf1, surf2) + #box1, tree1 = isec.bounding_boxes(surf1) + #box2, tree2 = isec.bounding_boxes(surf2) + #isec.get_intersection(surf1,surf2,tree1,tree2,box1,box2,isec.nt,isec.nit) # surf1,surf2,tree1,tree2 + point_list1, point_list2 = isec.get_intersection() + + m = point_list1.__len__() + point_list2.__len__() + + X= np.zeros([m]) + Y = np.zeros([m]) + Z = np.zeros([m]) + + i = -1 + + for point in point_list1: + i += 1 + X[i] = point.R3_coor[0] + Y[i]= point.R3_coor[1] + Z[i]= point.R3_coor[2] + + for point in point_list2: + i += 1 + X[i] = point.R3_coor[0] + Y[i]= point.R3_coor[1] + Z[i]= point.R3_coor[2] + + + + myplot.scatter_3d(X, Y, Z) + + + myplot.show() # view + + #print(tree1.find_box(boxes2[0])) + #print(surf1.poles[:,:,1]) + #print(surf1.u_basis.n_intervals) + #print(surf1.u_basis.knots) + + diff --git a/transform.m b/transform.m deleted file mode 100644 index 932957f..0000000 --- a/transform.m +++ /dev/null @@ -1,55 +0,0 @@ -function [ Xp X] = transform( X,P0,P1,P2,P3 ) - -[np k] = size(X); - -Pi = [P0 P1 P2 P3]; - -%%% Compute normalized normals - -Ni = zeros(2,4); - - -Ni(:,1) = Pi(:,1) - Pi(:,4); -nt = Ni(1,1); -Ni(1,1) = -Ni(2,1); -Ni(2,1) = nt; -Ni(:,1) = Ni(:,1)/norm(Ni(:,1)); - - - -for i =2:4 - Ni(:,i) = Pi(:,i) - Pi(:,i-1); - nt = Ni(1,i); - Ni(1,i) = -Ni(2,i); - Ni(2,i) = nt; - Ni(:,i) = Ni(:,i)/norm(Ni(:,i)); -end - -%%% Compute local coordinates and drop all - -h = 0; -for j=1:np - - P = [X(j,1); X(j,2)]; - - u = (P - Pi(:,1))' * Ni(:,1) / ( (P - Pi(:,1))' * Ni(:,1) + (P - Pi(:,3))' * Ni(:,3) ) ; - v = (P - Pi(:,1))' * Ni(:,2) / ( (P - Pi(:,1))' * Ni(:,2) + (P - Pi(:,4))' * Ni(:,4) ) ; -% P - Pi(:,1) -% P - Pi(:,4) -% Ni(:,2) -% Ni(:,4) -% (P - Pi(:,1))' * Ni(:,2) -% (P - Pi(:,4))' * Ni(:,4) -% pause - - if (u >= 0) && (u <= 1) && (v >= 0) && (v <= 1) - h = h+1; - Xp(h,1) = u; - Xp(h,2) = v; - Xp(h,3:4) = X(j,3:4); - Xp(h,5:6) = X(j,1:2); - end -end - -end - diff --git a/vec2mat.m b/vec2mat.m deleted file mode 100644 index 7f0d202..0000000 --- a/vec2mat.m +++ /dev/null @@ -1,11 +0,0 @@ -function [ Z_coor ] = vec2mat( zs,u_n_basf,v_n_basf ) -%UNTITLED3 Summary of this function goes here -% Detailed explanation goes here -Z_coor = zeros(u_n_basf,v_n_basf); - -for i =1:v_n_basf - Z_coor(:,i) = zs(1+(i-1)*(u_n_basf):i*(u_n_basf)); -end - -end -