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..3888434 --- /dev/null +++ b/src/brep_writer.py @@ -0,0 +1,1061 @@ +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 + return curve_dim(curve.poles, curve.basis.pack_knots(), curve.rational, curve.basis.degree) + +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 _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 _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: + """ + return Surface(surf.poles, (surf.u_basis.pack_knots(), surf.v_basis.pack_knots()), + surf.rational, (surf.u_basis.degree, surf.v_basis.degree) ) + + +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 _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 + 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 + 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 + """ + 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 + """ + 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 + """ + 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..9ad2166 --- /dev/null +++ b/src/bspline.py @@ -0,0 +1,971 @@ +""" +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 matplotlib.pyplot as plt +import numpy as np +import math +import numpy.linalg as la + + +__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(self, t): + """ + Evaluate a B-spline curve for paramater 't'. + :param t: Evaluation point. + :return: D-dimensional pnumpy array. D - is dimension given by dimension of poles. + TODO: + - use basis.eval_vector + - test evaluation for rational curves + """ + + it = self.basis.find_knot_interval(t) + dt = self.basis.degree + 1 + t_base_vec = np.array([self.basis.eval(jt, t) for jt in range(it, it + dt)]) + + 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_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.dim = len(poles[0][0]) - rational + # Surface dimension, D. + + self.poles=np.array(poles, dtype=float) + # 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(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. + TODO: + - use basis.eval_vector + - 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) + du = self.u_basis.degree + 1 + dv = self.v_basis.degree + 1 + u_base_vec = np.array([self.u_basis.eval(ju, u) for ju in range(iu, iu + du)]) + v_base_vec = np.array([self.v_basis.eval(jv, v) for jv in range(iv, iv + dv)]) + + 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_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. + + # Build envelope quadrilateral polygon in XY plane. + 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 + + + #TODO: remove this, after fixing GridSurface z_eval + #self.z_scale =1.0 + #self.z_shift = 1.0 + + 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) + poles = np.concatenate( (xy_poles, self.z_surface.poles), axis = 2 ) + + return Surface(basis, poles) + + + def transform(self, xy_mat, z_mat=np.array( [1.0, 0.0] ) ): + """ + 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: [ z_scale, z_shift] + :return: None + """ + assert xy_mat.shape == (2, 3) + assert z_mat.shape == (2, ) + + 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) + + # apply z-transfrom directly to the poles + self.z_surface.poles *= z_mat[0] + self.z_surface.poles += z_mat[1] + + #TODO: remove this, after fixing GridSurface z_eval + #self.z_scale *= z_mat[0] + #self.z_shift = self.z_shift*z_mat[0] + z_mat[1] + + + + """ + 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) + 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) + 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) + z_points = self.z_surface.eval_array(uv_points) + return np.concatenate( (xy_points, z_points), axis = 1) + + + 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) + 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) + z_points = self.z_surface.eval_array(uv_points) + return z_points.reshape(-1) + + + 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_surface.aabb()[:,0] + 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 ) + + # 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) ) + 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 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 + 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 + 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): + return self.z_surface.aabb() + + +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..2197212 --- /dev/null +++ b/src/bspline_approx.py @@ -0,0 +1,416 @@ +""" +Collection of functions to produce Bapline curves and surfaces as approximation of various analytical curves and surfaces. +""" + +import bspline as bs +import numpy as np +import math +import time +#import numpy.linalg +import scipy.sparse +import scipy.sparse.linalg +import scipy.interpolate + + +""" +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, **kwargs): + """ + 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(grid_surface, nuv, **kwargs) + return approx.get_approximation() + + +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 + + + + + + +class _SurfaceApprox: + """ + TODO: + - Check efficiency of scipy methods, compare it to our approach assuming theoretical number of operations. + - Optimize construction of A regularization matrix. + - 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 + + """ + + + + def __init__(self, grid_surface, nuv=None, **kwargs): + self.degree = np.array((2, 2)) + self.grid_surf = grid_surface + self.nuv = nuv + self.regularization_weight = kwargs.get('reg_weight', 0.001) + + def get_approximation(self): + if not hasattr(self, 'z_surf'): + self.z_surf = self.approx_chol() + return self.z_surf + + + def approx_chol(self): + """ + This function tries to approximate terrain data with B-Spline surface patches + using Cholesky decomposition + :param terrain_data: matrix of 3D terrain data + :param quad: points defining quadrangle area (array) + :param u_knots: array of u knots + :param v_knots: array of v knots + :param sparse: if sparse matrix is used + :param filter_thresh: threshold of filter + :return: B-Spline patch + """ + + + print('Transforming points to parametric space ...') + start_time = time.time() + + points_uvz = self.grid_surf.grid_uvz.reshape(-1, 3) + points_uv = points_uvz[:, 0:2] + + # 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] + self.points_z = points_uvz[in_idx, 2][:,None] + print("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.points_uv = np.minimum(points_uv, np.array([1.0, 1.0])) + + # determine number of knots + assert len(self.points_uv) == len(self.points_z) + n_points = len(self.points_uv) + + grid_shape = np.asarray(self.grid_surf.shape, dtype=int) + if self.nuv == None: + nuv = np.floor_divide( grid_shape / 3 ) + else: + nuv = np.floor(np.array(self.nuv)) + nuv = np.minimum( nuv, grid_shape - self.degree - 2) + size_u, size_v = nuv + self.degree + + # try to make number of unknowns less then number of remaining points + # +1 to improve determination + if (size_u + 1) * (size_v + 1) > n_points: + sv = np.floor(np.sqrt( n_points * size_v / size_u )) + su = np.floor(sv * size_u / size_v) + nuv = np.array( [su, sv] ) - self.degree + nuv = nuv.astype(int) + if nuv[0] < 1 or nuv[1] < 1: + raise Exception("Two few points, {}, to make approximation, degree: {}".format(n_points, self.degree)) + + print("Using {} x {} B-spline approximation.".format(nuv[0], nuv[1])) + self.u_basis = bs.SplineBasis.make_equidistant(2, nuv[0]) + self.v_basis = bs.SplineBasis.make_equidistant(2, nuv[1]) + + end_time = time.time() + print('Computed in {0} seconds.'.format(end_time - start_time)) + + + # Own computation of approximation + print('Creating B matrix ...') + start_time = time.time() + b_mat, interval = self.build_ls_matrix() + end_time = time.time() + print('Computed in {0} seconds.'.format(end_time - start_time)) + + print('Creating A matrix ...') + start_time = time.time() + a_mat = self.build_sparse_reg_matrix() + end_time = time.time() + print('Computed in {0} seconds.'.format(end_time - start_time)) + + print('Computing B^T B matrix ...') + start_time = time.time() + bb_mat = b_mat.transpose().dot(b_mat) + + end_time = time.time() + print('Computed in {0} seconds.'.format(end_time - start_time)) + + + print('Computing A and B svds approximation ...') + start_time = time.time() + + bb_norm = scipy.sparse.linalg.svds(bb_mat, k=1, ncv=10, tol=1e-4, which='LM', v0=None, + maxiter=300, return_singular_vectors=False) + a_norm = scipy.sparse.linalg.svds(a_mat, k=1, ncv=10, tol=1e-4, which='LM', v0=None, + maxiter=300, return_singular_vectors=False) + c_mat = bb_mat + self.regularization_weight * (bb_norm[0] / a_norm[0]) * a_mat + + g_vec = self.points_z[:, 0] + b_vec = b_mat.transpose().dot( g_vec ) + end_time = time.time() + print('Computed in {0} seconds.'.format(end_time - start_time)) + + + print('Computing 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." + + print(type(z_vec)) + end_time = time.time() + print('Computed in {0} seconds.'.format(end_time - start_time)) + + + print('Computing differences ...') + start_time = time.time() + diff = b_mat.dot(z_vec) - g_vec + max_diff = np.max(diff) + print("Approximation error (max norm): {}".format(max_diff) ) + end_time = time.time() + print('Computed in {0} seconds.'.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]) + z_surf = bs.Z_Surface(self.grid_surf.quad[0:3], surface_z) + + return z_surf + + + + 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.points_uv.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.points_uv[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._weights[j] * np.outer(u_base_vec,u_base_vec) + d_point_val_outer[:, :, i] += self._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 / math.sqrt(20)), (0.5 + 1 / math.sqrt(20)), 1] + self._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) + + 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) + + print("vnint: {} unint: {} nqp: {} prod: {}".format(v_n_inter, u_n_inter, nq_points, v_n_inter* u_n_inter* nq_points*nq_points)) + for i in range(v_n_inter): + for l in range(u_n_inter): + jac = 1.0 / u_n_inter / v_n_inter + idx_range = n_uv_loc_nz * n_uv_loc_nz + v_val_outer_loc = v_val_outer[:, :, i] + dv_val_outer_loc = v_diff_val_outer[:, : , i] + u_val_outer_loc = u_val_outer[:, :, i] + du_val_outer_loc = u_diff_val_outer[:, : , i] + data_m[nnz_a:nnz_a + idx_range] = jac * ( np.kron(v_val_outer_loc, du_val_outer_loc) + + np.kron(dv_val_outer_loc, u_val_outer_loc) ).ravel() + colv = np.repeat((i + linsp) * u_n_basf,self.u_basis.degree+1) + llinsp + 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..e24b7cf --- /dev/null +++ b/src/bspline_plot.py @@ -0,0 +1,137 @@ +""" +Functions to plot Bspline curves and surfaces. +""" + +import matplotlib.pyplot as plt +from matplotlib import cm +import numpy as np + + + +def plot_curve_2d(curve, n_points=100, **kwargs): + """ + Plot 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. + :return: Plot object. + """ + + 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) + + poles = kwargs.pop('poles', False) + img = plt.plot(x_coord, y_coord, **kwargs) + if poles: + plot_curve_poles_2d(curve) + + return img + + +def plot_curve_poles_2d(curve, **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 = curve.poles.T[0:2, :] # remove weights + return plt.plot(x_poles, y_poles, 'bo', color='red', **kwargs) + + + +def plot_surface_3d(surface, fig_ax, n_points=(100, 100), **kwargs): + """ + Plot a surface in 3d. + Usage: + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.pyplot as plt + fig = plt.figure() + ax = fig.gca(projection='3d') + plot_surface_3d(surface, ax) + plt.show() + + :param surface: Parametric surface in 3d. + :param fig_ax: Axes object: + :param n_points: (nu, nv), nu*nv - number of evaluation point + :param kwargs: surface_plot additional options + :return: The plot object. + """ + + 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. + poles = kwargs.pop('poles', False) + + img = fig_ax.plot_surface(X, Y, Z, **kwargs) + if poles: + plot_curve_poles_2d(surface) + return img + + + +def plot_grid_surface_3d(surface, fig_ax, n_points=(100, 100), **kwargs): + """ + Plot a surface in 3d. + Usage: + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.pyplot as plt + fig = plt.figure() + ax = fig.gca(projection='3d') + plot_surface_3d(surface, ax) + plt.show() + + :param surface: Parametric surface in 3d. + :param fig_ax: Axes object: + :param n_points: (nu, nv), nu*nv - number of evaluation point + :param kwargs: surface_plot additional options + :return: The plot object. + """ + + 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 + + X = X.reshape(U.shape) + Y = Y.reshape(U.shape) + Z = Z.reshape(U.shape) + + #Z = surface.z_eval_array(points).reshape(U.shape) + # Plot the surface. + img = fig_ax.plot_surface(U, V, Z, **kwargs) + return img + + +def plot_surface_poles_3d(surface, fig_ax, **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 fig_ax.scatter(x_poles, y_poles, z_poles, color='red', **kwargs) + + + 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..4ff7684 --- /dev/null +++ b/src/isec_surf_surf.py @@ -0,0 +1,8 @@ +class IsecSurfSurf: + ''' + 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. + ''' + \ No newline at end of file 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..2db43d3 --- /dev/null +++ b/tests/test_bs_approx.py @@ -0,0 +1,168 @@ +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 + +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): + fig = plt.figure() + ax = fig.gca(projection='3d') + ax.plot_surface(a_grid[:, :, 0], a_grid[:, :, 1], a_grid[:, :, 2], color='green') + ax.plot_surface(b_grid[:, :, 0], b_grid[:, :, 1], b_grid[:, :, 2], color='red') + plt.show() + + diff = b_grid - a_grid + fig = plt.figure() + ax = fig.gca(projection='3d') + ax.plot_surface(a_grid[:, :, 0], a_grid[:, :, 1], diff[:, :, 0], color='red') + ax.plot_surface(a_grid[:, :, 0], a_grid[:, :, 1], diff[:, :, 1], color='green') + ax.plot_surface(a_grid[:, :, 0], a_grid[:, :, 1], diff[:, :, 2], color='blue') + 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 plot_surf(self, surf): + fig = plt.figure() + ax = fig.gca(projection='3d') + bs_plot.plot_surface_3d(surf, ax) + plt.show() + + + def test_approx_func(self): + 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("Compare: Func - GridSurface.transform.Z_surf_approx") + #z_surf = bs_approx.surface_from_grid(gs, (16, 24)) + z_surf = bs_approx.surface_from_grid(gs, (100, 100)) + #z_surf.transform(xy_mat, z_mat) + 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) + + + + + + + #points = bs.make_function_grid(function_sin_cos, 20, 30) + #points = np.dot(points, mat.T) + np.array([10, 20, -5.0]) + #gs = bs.GridSurface(points.reshape(-1, 3)) + #z_surf = bs_approx.surface_from_grid(gs, (3,4) ) + #self.plot_surf(z_surf) + + + + + + + + + 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_surface_approx(self): + # self.plot_approx_grid() + # self.plot_approx_transformed_grid() + # self.plot_plane() + pass + + + + +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..762f8cc --- /dev/null +++ b/tests/test_bspline.py @@ -0,0 +1,428 @@ +import pytest +import bspline as bs +import numpy as np +import math +import bspline_plot as bs_plot + +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.pyplot as plt + +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 ] + plt.plot(x_coord, y_coord) + + plt.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): + plt.plot(x_coord, y_coords[i_base, :]) + + plt.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): + plt.plot(x_coord, y_coords[i_base, :]) + + plt.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) + + bs_plot.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]) + + plt.plot( bb[:, 0], bb[:, 1], color='green') + plt.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): + fig = plt.figure() + ax = fig.gca(projection='3d') + + # 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) + bs_plot.plot_surface_3d(surface_extrude, ax, poles = True) + plt.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) + bs_plot.plot_surface_3d(surface_func, ax) + bs_plot.plot_surface_poles_3d(surface_func, ax) + + plt.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 plot_function_uv(self): + # function surface + def function(x): + return math.sin(x[0]*4) * math.cos(x[1]*4) + + fig = plt.figure() + ax = fig.gca(projection='3d') + + 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], [0, 1], [1, 0], [1, 1]]) + z_surf = bs.Z_Surface(quad, surface_func) + full_surf = z_surf.make_full_surface() + + bs_plot.plot_surface_3d(z_surf, ax) + bs_plot.plot_surface_3d(full_surf, ax, color='red') + #bs_plot.plot_surface_poles_3d(surface_func, ax) + + plt.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 + + 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] ] ) + z_surf = bs.Z_Surface(quad, surface_func) + box = z_surf.aabb() + assert np.allclose(box, np.array([[0, 0, 3], [1.1, 1.1, 5]])) + + +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): + fig = plt.figure() + ax = fig.gca(projection='3d') + ax.plot_surface(XYZ_grid_eval[:, :, 0], XYZ_grid_eval[:, :, 1], XYZ_grid_eval[:, :, 2], color='blue') + ax.plot_surface(XYZ_surf_eval[:, :, 0], XYZ_surf_eval[:, :, 1], XYZ_surf_eval[:, :, 2], color='red') + ax.plot_surface(XYZ_func_eval[:, :, 0], XYZ_func_eval[:, :, 1], XYZ_func_eval[:, :, 2], color='green') + plt.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) 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 -