diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8941317 --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +# Upstream docs / package archives are large; keep Matlab sources + Python engine tracked. +reference/package/*.tar.gz +reference/package/*.zip +__pycache__/ +*.py[cod] +.venv/ +dist/ +*.egg-info/ +.pytest_cache/ diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..94a9ed0 --- /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. + + + Copyright (C) + + 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: + + Copyright (C) + 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 index 61d2b19..445cd6d 100644 --- a/README.md +++ b/README.md @@ -1,43 +1,37 @@ -Fractal-Mosaics -=============== +# Fractal Mosaics (Tower fork) -Fractal Mosaics is a rotation, scale, and translation invariant photomosaic algorithm (i.e., the images can be placed at arbitrary locations, can be scaled to any size, and can be rotated by any angle). +Fork of [s-ben/Fractal-Mosaics](https://github.com/s-ben/Fractal-Mosaics) adapted for the Tower **mosaik** drop-inbox. -The figure below shows a Fractal Mosaic and traditional photomosaic made using the same target image and image set. +Upstream implements a rotation / scale / translation invariant photomosaic algorithm (log-polar registration + CBIR prefilter + large-first placement). Original code is MATLAB (GPL-3.0). This repository keeps the unpacked MATLAB sources under `reference/upstream-matlab/` and ships a **Python CLI** for headless use on Arch/Linux. -![Photomosaic New Algorithm Fractal Eye](https://raw.github.com/s-ben/Fractal-Mosaics/gh-pages/images/fractal_trad_compare.png) +## Credit -For a detailed description of the algorithm, see my adademic-style "[paper](http://s-ben.github.io/Fractal-Mosaics/)" on it and documentation below. +- Original author: Seth Benton — https://github.com/s-ben/Fractal-Mosaics +- Algorithm write-up: https://s-ben.github.io/Fractal-Mosaics/ +- License: GNU GPL v3 (see `LICENSE` and `reference/upstream-matlab/gpl_license.txt`) -# Prerequisites +## Python CLI (Tower) -Fractal Mosaics is implemented in Matlab. You’ll need a copy of Matlab with the Signal Processing toolbox - most versions should have the Signal Processing toolbox. +```bash +python3 -m venv .venv && source .venv/bin/activate +pip install -e '.[cli]' +fractal-mosaic \ + --master path/to/target.jpg \ + --tiles path/to/tile_pool/ \ + --output out.png \ + --profile mittel +``` -# Getting Started +Profiles: `wenig` (fast, larger tiles), `mittel`, `dicht` (more placements, slower). -1. Clone the Fractal-Mosaics repository (or, if you don't have git, just download the [zip file](https://github.com/s-ben/Fractal-Mosaics/archive/gh-pages.zip)). -2. Navigate to the /code folder and unpack the 'FractalMosaics_package' zip or tar.gz file. This contains the demo example. -3. Download the [example image set](https://docs.google.com/file/d/0B_2ApIVBvXm1Y0VMVTVqNUhCaVk/edit?usp=sharing) (9.1 GB database of greyscale images of graffiti). -4. Unpack the zip file containing the example image set and copy its contents into the /library_images folder in the demo example. -5. Navigate to the /m-files directory in the demo example. -6. Run ‘fractal_mosaic.m’. +**Constraints for Tower v1:** geometric transforms + edge alpha only — **no tile recoloring**. -Your mosaic is now rendering. +## Layout -WARNING: it takes a long time to render a mosaic (~10.5 hrs on my macbook pro). If you want to optimize it, please do ;) - -The demo example creates a black and white mosaic. This is because color takes ~3X as long. To create a color mosaic, download the color image set and set Fractal Mosaics to color mode (see the [Variables](https://github.com/s-ben/Fractal-Mosaics/wiki/Variables) page for how to do this). Due to file size limits in Google Drive, the color image set is divided into two zip files ([image set 1](https://docs.google.com/file/d/0B3M2Bu_2k-QSSnhSVzNWRERtenc/edit?usp=sharing), [image set 2](https://docs.google.com/file/d/0B3M2Bu_2k-QSTkxxbEVGay1Gcjg/edit?usp=sharing)). - -# Documentation - -* [Fractal Mosaic “paper”](http://s-ben.github.io/Fractal-Mosaics/): academic-style paper explaining how the algorithm works. -* [Wiki](https://github.com/s-ben/Fractal-Mosaics/wiki): details on how to create your own mosaic. -* [Gallery](http://www.flickr.com/photos/travelingseth/sets/72157623789223762/): Flickr page showing finished mosaics. - -# Licencing - -Fractal Mosaics is licensed under the GNU General Public License ([GPL](http://www.gnu.org/licenses/gpl.txt)). Some functions used in Fractal Mosaics are from the Matlab File Exchange and use the [BSD license](http://opensource.org/licenses/bsd-license.php). This license is "permissive" (can be used for any purpose, commercial or non). If you wish to use Fractal Mosaics for proprietary software, let's talk! (seth.benton@gmail.com) - - -[![githalytics.com alpha](https://cruel-carlota.pagodabox.com/80d211511c492cf9dce9dd9841acf603 "githalytics.com")](http://githalytics.com/s-ben/Fractal-Mosaics) +| Path | Role | +|------|------| +| `reference/upstream-matlab/` | Unpacked `.m` files from upstream package | +| `reference/package/` | Original `FractalMosaics_package.tar.gz` | +| `src/fractal_mosaic/` | Python engine + CLI | +| docs site (`index.html`, …) | Upstream gh-pages documentation (preserved) | diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..24e8d61 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,26 @@ +[build-system] +requires = ["setuptools>=68", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "fractal-mosaic" +version = "0.1.0" +description = "RST collage / fractal photomosaic engine (fork of s-ben/Fractal-Mosaics) for Tower inbox" +readme = "README.md" +license = { text = "GPL-3.0-or-later" } +requires-python = ">=3.10" +authors = [{ name = "Tower / copilot" }] +dependencies = [ + "numpy>=1.24", + "Pillow>=10.0", + "scipy>=1.10", +] + +[project.optional-dependencies] +cli = [] + +[project.scripts] +fractal-mosaic = "fractal_mosaic.cli:main" + +[tool.setuptools.packages.find] +where = ["src"] diff --git a/reference/upstream-matlab/DC_calc.m b/reference/upstream-matlab/DC_calc.m new file mode 100644 index 0000000..4a1d3a5 --- /dev/null +++ b/reference/upstream-matlab/DC_calc.m @@ -0,0 +1,68 @@ + +clear all +close all + +color = 0; +% ----------------------------- get jpeg file paths -------------------- + +'getting image paths...' + + if color ==0 + path = ['../library_images/graffiti_256_black_white/']; + else + path = ['C:\Americano\library_images\graffiti_256_color\']; + end + + + % files = dir('C:\Image_registration\flickr_spider\graffiti_library\*.jpg'); % get file names + directories = dir(path); % get directory names + + num_files = 0; + file_count = 0; + + dir_temp = directories(3).name; % get first directory name + files = dir([path,dir_temp,'/*.jpg']); % get file names for first directory + file_count = length(files); + + % add directory path to jpeg files name + for j=1:length(files) + files(j).name = [path,dir_temp,'/',files(j).name]; + end + + for i=2:(length(directories)-2) % offset by 2 because first two entries are '.' and '..' + + + dir_temp = directories(i+2).name; % get next directory name + + files_temp = dir([path,dir_temp,'/*.jpg']); % get file names for given directory + + % add directory path to jpeg files name + for j=1:length(files_temp) + files_temp(j).name = [path,dir_temp,'/',files_temp(j).name]; + end + + + + % concatenate file paths from next directory to existing + files((file_count+1):(file_count+length(files_temp))) = files_temp; + + + file_count = file_count+length(files_temp); + % pause + end + +DC_v = zeros(length(files),1); + +for z=1:length(files) + +% z + Im_lib = imread([files(z).name]); % read in library image + + DC_component = sum(sum(Im_lib))/(256*256); + + DC_v(z) = DC_component; + + +end + +save DC_vals_graffiti DC_v diff --git a/reference/upstream-matlab/fractal_mosaic.m b/reference/upstream-matlab/fractal_mosaic.m new file mode 100644 index 0000000..5228997 --- /dev/null +++ b/reference/upstream-matlab/fractal_mosaic.m @@ -0,0 +1,1127 @@ +% Copyright (c) 2013 Seth Benton (seth.benton@gmail.com) +% +% This program (Fractal Mosaics) 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 . + +% The latest source code and documentation for Fractal Mosaics can be +% found at: https://github.com/s-ben/Fractal-Mosaics + + + +clear all +close all +load DC_vals_graffiti + +tic + +% ---------- User Parameters ---------------- + +% Use these parameters to control the behavior of the Fractal Mosaic +% algorithm + +color = 0; % if color = '0', black and white images are used for matching. color = '1' uses color images +threshold_layer_on = 0; % '1' if a "transparency layer" (thresh_layer) is used to control the acceptable error (rms_bias_v) on a pixel-by-pixel basis. + % if threshold_layer_on = '0', the same error + + +target_path = ['rose_black_bg__corner_2048.jpg']; % path to target image +threshold_path = ['eye_pic_clean_thresh_layer.jpg']; % path to (optional) transparency image + +% path = ['../library_images/san_francisco_256_color/']; % path to library images +path = ['../library_images/graffiti_256_black_white/']; % path to library images + + +bias_floor = 24; % bias floor +rms_bias_range = 6; % range of scaling according to transparency layer (0-1 greyscale image) + +max_repeats_v = [1 2 4 4 20 50 100 100]; % maximum number of repeated images in the mosaic (level dependent) +keeper_thresh_v = [10 10 10 11 14 15 15 15]; % for each of the eight "levels" [1 2 3 4 5 6 7 8] (scales of analysis), set an rms threshold. + +num_alt_images = 3; % number of alternate images to store per pixel + + +perc_max_overlap_v = [.20 .20 .20 .20 .30 .35 .40 .50]; % max overlap of candidate images on placed image + +num_lib_test_v = [1 1 1 1 1 1 1 1]; % number of candidate matches to stop at at each level (i.e., if '1' for a given level, the algorithm won't continue to look +% "underneath" a placed pixel for more candidate images. You probably +% won't touch this parameter + +num_lib_test_pix_v = [3000 2500 2000 1000 750 1000 500 500]; % number of library images to test for each pixel, per level + +num_lib_test_pix_v = num_lib_test_pix_v/2; % REMOVE + + +% ------- resolution levels ------------ +res_v = [8 16 32 64 128 256]; +level_v = [5 5 5 5 5 5 4 3]; % for each level, the res_v() to go up to + + +% Create matrix for calculation of color histograms +range = 0.0:0.1:1.0; +rangeNew = 0.0:0.05:1.0; +[x1,y1,z1] = meshgrid(range); +[x2,y2,z2] = meshgrid(rangeNew); + + +% --------- rendering parameters ------------------- +rendering_level = 5; % indexes res_v (e.g., 1 -> 8x8, 2 -> 16x16, ....) +scaling_constant_v = [7 8 9 10 11 12 13 14]; % vector of values used to scale candidate pixels to their respective analysis levels (j) +lowest_res_image = res_v(rendering_level);%128; % used below to calculate size of "analysis tiles" for each Level + +% --------- calculate size of analysis tile for each level --------------------- +for i=1:length(level_v) + lib_im_v(i) = lowest_res_image*2^(length(level_v)-2-i); +end + + +% Define vector for storing rms error values +rms_err_v = zeros(length(max(num_lib_test_pix_v)),1); + + +% Calculate rms thresholds for each Level (i) and layer (j) +rms_thresh_v = [(bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+1) ; % i=1 + (bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+1); % i=2 + (bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+1); % i=3 + (bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+3); % i=4 + (bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+4); % i=5 + (bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+6); % i=6 + (bias_floor+4) (bias_floor+4) (bias_floor+4) (bias_floor+7) (bias_floor+7); % i=7 + (bias_floor+7) (bias_floor+7) (bias_floor+10) (bias_floor+10) (bias_floor+10)]; % i=8 + + +% number of candidate pixels (via level 'i' and resolution 'j') that are passed through to the next level +pass_through_v_abs = [1 4 4 1 1 1 ; 1 4 4 1 1 1 ; 1 4 4 1 1 1; 1 4 4 1 1 1; 1 4 1 1 1 1; 1 1 1 1 1 1; 1 1 1 1 1 1; 1 1 1 1 1 1]; + + +% ------------- Load target image ------------------ + +if color == 0 + + I1 = imread(target_path); + + I1_size = size(I1); + + % If target an RGB image, converty to greyscale + if (I1_size(3) == 3) + I1 = rgb2gray(I1); + end + + if threshold_layer_on == 1 + thresh_layer = imread(threshold_path); + thresh_layer = rgb2gray(thresh_layer); + end + + + +else +% I1 = imread('old_car_blue_2048_clean.jpg'); +% thresh_layer = imread('old_car_blue_2048_clean_transparency.jpg'); + + I1 = imread(target_path); + I1_size = size(I1); + + + if threshold_layer_on == 1 + thresh_layer = imread(threshold_path); + end + +end + + + +Im_target = I1; +clear I1; + +% ----------------------------- get jpeg file paths -------------------- + +'getting image paths...' + + + directories = dir(path); % get directory names + + num_files = 0; + file_count = 0; + + dir_temp = directories(3).name; % get first directory name + files = dir([path,dir_temp,'/*.jpg']); % get file names for first directory + file_count = length(files); + + % add directory path to jpeg files name + for j=1:length(files) + files(j).name = [path,dir_temp,'/',files(j).name]; + end + + for i=2:(length(directories)-2) % offset by 2 because first two entries are '.' and '..' + + + dir_temp = directories(i+2).name; % get next directory name + + files_temp = dir([path,dir_temp,'/*.jpg']); % get file names for given directory + + % add directory path to jpeg files name + for j=1:length(files_temp) + files_temp(j).name = [path,dir_temp,'/',files_temp(j).name]; + end + + % concatenate file paths from next directory to existing + files((file_count+1):(file_count+length(files_temp))) = files_temp; + + file_count = file_count+length(files_temp); + end + + +% --- create varying resolution target images for matching ------ +if color == 0 + size_target = size(Im_target); +else + size_target = size(Im_target(:,:,1)); +end + +partitions = I1_size(1)/1024; + +res_scale = (res_v(1)*partitions)/size_target(1); % height of lowest res target image for matching / orig target size +Im_target_lr = imresize(Im_target,res_scale); % lowest resolution target image +size_target_low_res = size(Im_target_lr); + + +for i=1:10 % for each resolution needed (hard-coded for now, need to generalize) + + % ------- resize target to res ------- + + scale_factor = 2^(i-1); + Im_target_tmp(i).lr = imresize(Im_target,[size_target_low_res(1)*scale_factor size_target_low_res(2)*scale_factor]); + +end + +% ---------- create final 'Render_target' variable ------ +Render_target = Im_target_tmp(10).lr; +size_render_target = size(Render_target); + +% Render_target_valid is a logical bitmap specifying whether an images has been "placed" on the target image. +% declar all '0s' logical image +Render_target_valid = im2bw(Render_target,1); +% --------------------------------------------------------------- + +if color ==0 + + % -------- load histograms of library images ------------ + load('hist_matrix_graffiti_128'); +% load('hist_matrix_SF_128'); + % ------------- Load pdfs for hist_dist, gfd, std_diff ------------------ + load('pdfs'); + + % -------- load GFDs of library images ------------ + load('gfd_matrix_graffiti_128'); +% load('gfd_matrix_SF_128'); + + load('std_matrix_graffiti_128'); +% load('std_matrix_SF_128'); + +else + load san_francisco_color_hists; + %load graffiti_color_hists; +end + + +% -------- resize threshold layer (assuming i=8 is highest) ---------------- +Im_size_lr = size(Im_target_tmp(8).lr); +if color == 0 + if threshold_layer_on == 1 + thresh_layer = imresize(thresh_layer, Im_size_lr); + end + +else + if threshold_layer_on == 1 + thresh_layer = imresize(thresh_layer, [Im_size_lr(1) Im_size_lr(2)]); + end +end + +% ---------- calculate distance windows ----------- + +for i=1:length(level_v) + + perc_sloped = .05; + perc_sloped_fixed = .01; % perc sloped for fixed window (one that's used to cut off clipping) + + fd = zeros(lib_im_v(i)); + fd_fixed = zeros(lib_im_v(i)); + + y_d = ones(1,lib_im_v(i)/2); + y_d_fixed = ones(1,lib_im_v(i)/2); + + x_ind = linspace(0,lib_im_v(i)/2,lib_im_v(i)/2); + for u=1:(perc_sloped*(lib_im_v(i)/2)) + y_d(u) = (1/(perc_sloped*(lib_im_v(i)/2)))*x_ind(u); + end + for u=1:(perc_sloped_fixed*(lib_im_v(i)/2)) + y_d_fixed(u) = (1/(perc_sloped_fixed*(lib_im_v(i)/2)))*x_ind(u); + end + + fd(1:lib_im_v(i)/2) = y_d; + fd(lib_im_v(i)/2+1:lib_im_v(i)) = fliplr(y_d); + fd_2 = fd*fd'; + window(i).win = fd_2; + + fd_fixed(1:lib_im_v(i)/2) = y_d_fixed; + fd_fixed(lib_im_v(i)/2+1:lib_im_v(i)) = fliplr(y_d_fixed); + fd_2_fixed = fd_fixed*fd_fixed'; + window(i).win_fixed = fd_2_fixed; + +end + + + +% ----- allocate array to store candidate pixels ------- +cand_pixels = zeros(2,4^(length(res_v)-1)); +cand_pixels_next = zeros(2,4^(length(res_v)-1)); +cand_pixels_rms = zeros(2,4^(length(res_v)-1)); + + + +'rendering mosaic' + +% set counters to zero +num_cand_pixels_final = 0; +num_cand_pixels_alt = 0; +cand_pixels_final_abs = zeros(10000,3); +dist = zeros(1,30000); + +% calculate histogram bins +bins = 16; +hist_x = [(255/bins)/2:255/bins:(255-(255/bins)/2)]; + +% declare final pdf (sum of pdfs) array +match_pdf = zeros(length(files),1); + +% declare sq_center for each level (i.e., declare pixel boundaries) +% the pixel boundaries (for now) manually squeeze in to avoid too many +% lower level (time consuming) matches on the borders +sq_center_i_v = [5 5 5 7 14 28 56 112]; + +% vector of scalings for rad_buff (establishes radius for palced matches) +buff_v = [3 3 3.5 3.8 3.9 4 4.2 4.5]; + + + + +for i=2:length(level_v) % for each scale (level) + + +% 'scale' +% i + + % ------- read lowest resolution target image for given scale ------- + Im_target_match = Im_target_tmp(i).lr; + + if color ==0 + size_Im_target_match = size(Im_target_match); + else + size_Im_target_match = size(Im_target_match(:,:,1)); + end + % ------- read in target image w/ 128 resolution for current scale ------ + Im_target_match_hist = Im_target_tmp((i-1)+level_v(i)).lr; + + + + % ----- for each pixel in valid area --------- + for y = sq_center_i_v(i):(size_Im_target_match(1)-sq_center_i_v(i)+1) % for each "translation" pixel in (lowest res) valid area (y_sample(v)*(2^(i-1))):(y_sample(v)*(2^(i-1))) + for x = sq_center_i_v(i):(size_Im_target_match(2)-sq_center_i_v(i)+1) %sq_center_i(2):(size_Im_target_match(2)-sq_center_i(2)+1) + + + img_tested = num_lib_test_v(i); % number of library images at this level to evaluate + num_img_in_rad = 0; % number of images within the radius of already selected candidate images + quota_reached = 0; % binary variable to flag if quota has been reached + cand_image_found = 0; + num_cand_image_found = 0; % number of candidate images found + + % ---------- scale cand pixel position according to scale and res --------------------- + x_tmper = round(x*2^(scaling_constant_v(rendering_level)-(i+1))); + y_tmper = round(y*2^(scaling_constant_v(rendering_level)-(i+1))); + + % ------------ find rms bias for given pixel ------------------------ + x_thresh = x*2^(8-i); + y_thresh = y*2^(8-i); + + + % ------------------------------------------------------------- + + % If threshold layer turned off, don't vary rms error threshold + % pixel-by-pixel + if (threshold_layer_on == 0) + rms_bias_v = 0; + else + rms_bias_v = thresh_layer(y_thresh,x_thresh); + rms_bias_v = rms_bias_range*(1-(double(rms_bias_v)/255)); + end + + + % ------------- calculate distance of cand pixel to placed images -------------------------------- + for w=1:num_cand_pixels_final + + % If a library image has already been placed which covers + % the candidate pixel, exit loop + if (Render_target_valid(y_tmper,x_tmper) == 1) + quota_reached = 1; + break + end + + % calculate the minimum distance from the center of a place + % image that a candidate image needs to be to be evaluated + min_dist = cand_pixels_final_abs(w,3) - perc_max_overlap_v(i)*lib_im_v(i) + lib_im_v(i)/2; + + % calculate distance from candidate pixel location to + % placed image location + dist(w) = sqrt((cand_pixels_final_abs(w,1) - x_tmper)^2 + (cand_pixels_final_abs(w,2) - y_tmper)^2); + + % if candidate pixel is too close to a placed image, exit + % loop + if (dist(w) < min_dist) + num_img_in_rad = num_img_in_rad+1; + if (num_img_in_rad >= num_lib_test_v(i)) + quota_reached = 1; + break + end + end + end + + % --------------- if quota reached for cand pixels in this area, skip loop (don't evaluate this cand pixel) ------------------------- + if (quota_reached == 1) +% 'quota reached, exiting loop' + continue +% z = 1000000; + else + img_tested = num_lib_test_v(i)-num_img_in_rad; + z = 1; + end + + + % ------- calculate 128x128 target_tile for current pixel -------- + x_tmp = x; + y_tmp = y; + + for z=2:(level_v(i)) + % ----update pixel location for higher res --------- + x_tmp = 2*(x_tmp-0.5); + y_tmp = 2*(y_tmp-0.5); + end + + xmin = x_tmp-(floor((res_v(level_v(i)))/2)+1) +1; + ymin = y_tmp-(floor((res_v(level_v(i)))/2)+1) +1; + width = res_v(level_v(i))-1; + height = res_v(level_v(i))-1; + + rect = [xmin ymin width height]; + Im_target_tile_hist = imcrop(Im_target_match_hist,rect); + + if (res_v(level_v(i))< res_v(rendering_level)) + Im_target_tile_hist = imresize(Im_target_tile_hist,[res_v(rendering_level) res_v(rendering_level)]); + end + % ------------------------------------------------------------- + + + % ------------- calculate hist, std, and GFD of 'target_tile' ----------------------- + + % Parameters for calculating General Fourier Descriptor (GFD) + radial_freq = 6; + angular_freq = 6; + size_v = [128 128]; + + + if color==0 + + % -------- calculate DC component -------------- + DC_component = sum(sum(Im_target_tile_hist))/(128*128); + + % ----- calculate histogram ---------- + Im_lib_array = reshape(Im_target_tile_hist,128*128,1); + hist_target_tile = hist(double(Im_lib_array),hist_x); + + if (max(max(Im_target_tile_hist)) == 0) + Im_target_tile_hist = Im_target_tile_hist + 0.01; + end + + % ------ calculate GFD ------------- + Im_lib_polar = polar_conversion(Im_target_tile_hist); + Im_lib_polar_FFT = abs(fftshift(fft2(Im_lib_polar))); + + DC = Im_lib_polar_FFT(size_v(1)/2+1,size_v(2)/2+1); + Im_lib_polar_FFT = Im_lib_polar_FFT./DC; + center_gfd = [size_v(1)/2+1 size_v(2)/2+1]; + + gfd = Im_lib_polar_FFT(center_gfd(1):(center_gfd(1)+angular_freq-1),center_gfd(2):(center_gfd(2)+radial_freq-1)); + + % ------ calculate std ------------- + std_tile = std2(Im_target_tile_hist); + + else + + % -------- Get color image histogram ---------------- + + RGB = Im_target_tile_hist; + RGBt = RGB; + RGB = rgb2hsv(RGB); + + % get image size: + [M,N,ttt] = size(RGB); + + range = 0.0:0.1:1.0; + + Hist = zeros(length(range),length(range),length(range)); + + for ii=1:M + for jj=1:N + + nn1 = round(RGB(ii,jj,1) * 10)+1; + nn2 = round(RGB(ii,jj,2) * 10)+1; + nn3 = round(RGB(ii,jj,3) * 10)+1; + + Hist(nn1, nn2, nn3) = Hist(nn1, nn2, nn3) + 1; + + end + end + + Hist = Hist / (M*N); + + % Calculate DC component + img_tmp = rgb2gray(Im_target_tile_hist); + DC_component = sum(sum(img_tmp))/(256*256); + + % ---------------------------------------------------------- + end + + % ------------ calculate statistical difference (hist, gfd, std) --------------- + if color ==0 + + % ------------ calculate histogram distance --------------- + neighbor_distances = kNearestNeighbors_no_sort(hist_matrix,hist_target_tile,length(files)); + + % ------------ calculate GFD distance --------------- + neighbor_distances_gfd = kNearestNeighbors_no_sort(gfd_matrix,gfd(:)',length(files)); + + % ------------ calculate std distance --------------- + neighbor_distances_std = abs(std_tile - std_matrix); + +% neighbor_distances_int = neighbor_distances; + + else + + Nfiles = length(Hists); + + % decision thresholds: + t = 0.010; + t2 = 0.8; + + Hist = interp3(x1,y1,z1,Hist,x2,y2,z2); + Similarity = zeros(Nfiles, 1); + + for (w=1:Nfiles) % for each file in database: + + % compute (normalized) eucledean distance for all hist bins: + HistT = Hists{w}; + % HistT = interp3(x,y,z,Hists{i},x2,y2,z2); + DIFF = abs(Hist-HistT) ./ Hist; + + % keep distance values for which the corresponding query image's values + % are larger than the predefined threshold: + DIFF = DIFF(Hist>t); + + % keep error values which are smaller than 1: + DIFF2 = DIFF(DIFF= max_repeats_v(i)) + z=z+1; + end + end + end + + + image_path = [files(nearest_neighbor_int_index(z)).name]; + for j=1:length(image_path) + if (image_path(j) == '\') + image_path(j) = '/'; + end + end + + + Im_lib = imread(image_path); % read in library image + + + perc_DC = 1; % percentage DC adjusted + DC_adjust = (DC_v(nearest_neighbor_int_index(z))-DC_component)*perc_DC; + Im_lib = Im_lib-DC_adjust; + + + for j=1:level_v(i) % for each resolution (up to the highest acceptable for each scale) + + + % ------- read in target image w/ correct resolution ------ + Im_target_match = Im_target_tmp((i-1)+j).lr; + + % -------resize library image------------- + Im_lib_tmp = imresize(Im_lib, [res_v(j) res_v(j)]); +% Im_lib_tmp = Im_lib; + + % ------ calculate valid pixel area ------- + sq_center = [floor(res_v(j)/2)+1 floor(res_v(j)/2)+1]; + + % ------ get size of current res target image ------- + if color == 0 + size_target = size(Im_target_match); + else + size_target = size(Im_target_match(:,:,1)); + end + + % calculate the number of pixels to evaluate at this + % resolution + if j==1 + num_cand_pixels = 1; + else + [rms_sort rms_index] = sort(cand_pixels_rms(1:num_cand_pixels_new)); + num_cand_pixels = min(num_cand_pixels_new, pass_through_v_abs(i,j)); + end + + + + if num_cand_pixels == 0 + break + end + + + % ---------- assign candidates from last pass for anaylsis ------------ + if j>1 + for b=1:num_cand_pixels + % cand_pixels(1,b) = cand_pixels_next(1,b); + % cand_pixels(2,b) = cand_pixels_next(2,b); + + cand_pixels(1,b) = cand_pixels_next(1,rms_index(b)); + cand_pixels(2,b) = cand_pixels_next(2,rms_index(b)); + + end + end + + +% % ------------- calculate log-polar of target tile ---------------------------- + + + % ---- log polar transform library image ------ + if color==0 + L2 = transformImage_log(Im_lib_tmp, res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + else + L2_r = transformImage_log(Im_lib_tmp(:,:,1), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + L2_g = transformImage_log(Im_lib_tmp(:,:,2), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + L2_b = transformImage_log(Im_lib_tmp(:,:,3), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + end + + num_cand_pixels_new = 0; + + for b=1:num_cand_pixels + + % ---------- if desired number of candidate images have been reached, break from search ----------------------- + + + % ----- update pixel position at higher resolutions (if necessary) ----- + if j==1 + x_c = x; + y_c = y; + else + + %--- read in candidate pixel location--- + y_c = cand_pixels(1,b); + x_c = cand_pixels(2,b); + % ----update pixel location for higher res --------- + x_c = 2*(x_c-0.5); + y_c = 2*(y_c-0.5); + end + + + % ------ determine indexing parameters for higher resolution (if necessary) ------ +% res_scale_up = res_v(j)/res_v(1); + if j==1 + x_start = x_c; + x_end = x_c; + y_start = y_c; + y_end = y_c; + pix_dist = 1; + else + pix_dist = 2; + x_start = x_c; + x_end = x_start + pix_dist-1; + y_start = y_c; + y_end = y_start + pix_dist-1; + end + + + % ----- for each sub-pixel (if not lowest res) -------------- + for p = y_start:y_end + for q = x_start:x_end + + + if color==0 + Im_target_tile_masked = zeros([res_v(j) res_v(j)]); + else + Im_target_tile_masked = zeros([res_v(j) res_v(j) 3]); + end + valid_pix = 0; + mask = zeros([res_v(j) res_v(j)]); + + xmin = q-sq_center(2) +1; + ymin = p-sq_center(1) +1; + width = res_v(j)-1; + height = res_v(j)-1; + + rect = [xmin ymin width height]; + Im_target_tile = imcrop(Im_target_match,rect); + + + if color==0 + L1 = transformImage_log(Im_target_tile, res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + else + L1_r = transformImage_log(Im_target_tile(:,:,1), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + L1_g = transformImage_log(Im_target_tile(:,:,2), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + L1_b = transformImage_log(Im_target_tile(:,:,3), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid'); + end + + if color==0 + cc=ifft2(fft2(L1).*conj(fft2(L2))); + cc=flipud(cc); + cc = circshift(cc,[(1) (0)]); + else + cc=ifft2(fft2(L1_r).*conj(fft2(L2_r)))+ifft2(fft2(L1_g).*conj(fft2(L2_g)))+ifft2(fft2(L1_b).*conj(fft2(L2_b))); + cc=flipud(cc); + cc = circshift(cc,[(1) (0)]); + end + + % ----------- set up scale finding function ---------- + minScale = 1; + maxScale = res_v(j)/2-1; + rad_lin = linspace(log(minScale), log(maxScale),res_v(j))'; + rho = exp(rad_lin); + + [min_difference, array_position] = min(abs(rho - 2)); + + + % -------- find max of cross correlation --------------------- + + [y_max, x_ind] = max(cc(1:array_position,:)); % returns max values per row along with x-index value + [x_max, y_ind] = max(y_max); + y_max_index = x_ind(y_ind); + x_max_index = y_ind; + + % --------- calculate rotation degree ------------------------- + Theta = 360*((x_max_index -1)/res_v(j)); + + + + % --------- calculate scale value ---------------------------- + scale_est = rho(y_max_index); + + I_sc_tmp = imresize(Im_lib_tmp, 1/scale_est); + I_sc_tmp_size = size(I_sc_tmp); + + tmp_height = I_sc_tmp_size(1); + tmp_width = I_sc_tmp_size(2); + + mask_tmp = ones([I_sc_tmp_size(1) I_sc_tmp_size(2)]); + + if color==0 + I_sc = zeros([res_v(j) res_v(j)]); + I_comb = zeros([res_v(j) res_v(j)]); + else + I_sc = zeros([res_v(j) res_v(j) 3]); + I_comb = zeros([res_v(j) res_v(j) 3]); + end + + + if color==0 + I_sc(round(res_v(j)/2 - tmp_height/2)+1:round(res_v(j)/2 - tmp_height/2)+tmp_height,round(res_v(j)/2 - tmp_width/2)+1:round(res_v(j)/2 - tmp_width/2)+tmp_width) = I_sc_tmp; + else + I_sc(round(res_v(j)/2 - tmp_height/2)+1:round(res_v(j)/2 - tmp_height/2)+tmp_height,round(res_v(j)/2 - tmp_width/2)+1:round(res_v(j)/2 - tmp_width/2)+tmp_width,:) = I_sc_tmp; + end + mask(round(res_v(j)/2 - tmp_height/2)+1:round(res_v(j)/2 - tmp_height/2)+tmp_height,round(res_v(j)/2 - tmp_width/2)+1:round(res_v(j)/2 - tmp_width/2)+tmp_width) = mask_tmp; + + % -------Rotate image back by theta estimate -------- + + + R1 = imrotate(I_sc, -Theta, 'bilinear', 'crop'); + mask = imrotate(mask, -Theta, 'bilinear','crop'); + + alpha = 0.5; + I_comb = (alpha)*double(Im_target_tile) + (1-alpha)*R1; + + + + if color==0 + Im_target_tile_masked = double(Im_target_tile).*mask; + else + Im_target_tile_masked(:,:,1) = double(Im_target_tile(:,:,1)).*mask; + Im_target_tile_masked(:,:,2) = double(Im_target_tile(:,:,2)).*mask; + Im_target_tile_masked(:,:,3) = double(Im_target_tile(:,:,3)).*mask; + end + + rms_err = sqrt( (sum(sum( (abs(R1-Im_target_tile_masked)).^2 )))/(sum(sum(mask)))); + + + + if color==1 + rms_err = sum(rms_err)/3; + end + + rms_err_v(z,j) = rms_err; + + + + if rms_err < (rms_thresh_v(i,j)-rms_bias_v) + + num_cand_pixels_new = num_cand_pixels_new+1; + + + % -------store candidate pixels for next pass --------- + cand_pixels_next(1,num_cand_pixels_new) = p; + cand_pixels_next(2,num_cand_pixels_new) = q; + + cand_pixels_rms(num_cand_pixels_new) = rms_err; + + + + % ---------- If candidate pixel is in highest res for given scale, save ------------ + if (j == level_v(i)) + + + + % ---------- scale cand pixel position according to scale and res --------------------- + x_tmper = round(q*2^(scaling_constant_v(rendering_level)-(i+j))); + y_tmper = round(p*2^(scaling_constant_v(rendering_level)-(i+j))); + + radius = ((lib_im_v(i)/2)/scale_est); + + + if rms_err < max_rms_final + + % ----------- if >1 candidates, store alternatives-------------- + + if ((num_cand_image_found > 1)&& (num_cand_image_found <= num_alt_images)) + num_cand_pixels_alt = num_cand_pixels_alt+1; + + cand_pixels_final_alt(num_cand_pixels_alt).DC = cand_pixels_final(num_cand_pixels_final+1).DC; + cand_pixels_final_alt(num_cand_pixels_alt).rms = cand_pixels_final(num_cand_pixels_final+1).rms; + cand_pixels_final_alt(num_cand_pixels_alt).x_cand = cand_pixels_final(num_cand_pixels_final+1).x_cand; + cand_pixels_final_alt(num_cand_pixels_alt).y_cand = cand_pixels_final(num_cand_pixels_final+1).y_cand; + cand_pixels_final_alt(num_cand_pixels_alt).res = cand_pixels_final(num_cand_pixels_final+1).res; + cand_pixels_final_alt(num_cand_pixels_alt).scale = cand_pixels_final(num_cand_pixels_final+1).scale; + cand_pixels_final_alt(num_cand_pixels_alt).rot = cand_pixels_final(num_cand_pixels_final+1).rot; + cand_pixels_final_alt(num_cand_pixels_alt).lib_im = cand_pixels_final(num_cand_pixels_final+1).lib_im; + cand_pixels_final_alt(num_cand_pixels_alt).level = cand_pixels_final(num_cand_pixels_final+1).level; + cand_pixels_final_alt(num_cand_pixels_alt).valid_pix = cand_pixels_final(num_cand_pixels_final+1).valid_pix; + cand_pixels_final_alt(num_cand_pixels_alt).file_path = cand_pixels_final(num_cand_pixels_final+1).file_path; + + end + + cand_pixels_final(num_cand_pixels_final+1).DC = DC_adjust; + + cand_pixels_final(num_cand_pixels_final+1).rms = rms_err; + cand_pixels_final(num_cand_pixels_final+1).x_cand = q; + cand_pixels_final(num_cand_pixels_final+1).y_cand = p; + cand_pixels_final(num_cand_pixels_final+1).res = j; + cand_pixels_final(num_cand_pixels_final+1).scale = scale_est; + cand_pixels_final(num_cand_pixels_final+1).rot = Theta; + cand_pixels_final(num_cand_pixels_final+1).lib_im = nearest_neighbor_int_index(z); + cand_pixels_final(num_cand_pixels_final+1).level = i; + cand_pixels_final(num_cand_pixels_final+1).valid_pix = valid_pix; + cand_path = files(nearest_neighbor_int_index(z)).name; + cand_pixels_final(num_cand_pixels_final+1).file_path = cand_path(length(path):end); + + + + max_rms_final = rms_err; + + + cand_image_found = 1; + num_cand_image_found = num_cand_image_found+1; % count number of images below rms threshold for given pixel + + cand_pixels_final_abs(num_cand_pixels_final+1,1) = x_tmper; + cand_pixels_final_abs(num_cand_pixels_final+1,2) = y_tmper; + cand_pixels_final_abs(num_cand_pixels_final+1,3) = radius; + + + end + end + end + +% figure(1),subplot(2,2,1),imshow(Im_lib_tmp), title('Im_lib') +% subplot(2,2,2),imshow(Im_target_tile), title('Target Tile') +% subplot(2,2,3),imshow(uint8(R1)), title('transformed Image') +% subplot(2,2,4),imshow(uint8(I_comb)), title('Combined Images') +% pause +% +% rms_err + + end + end + end + end + + + z = z+1; + + end + + + % ------------ if highest res math found, increment final pixel count, Render match ------------- + if (cand_image_found ==1) + + num_cand_pixels_final = num_cand_pixels_final+1; + + % ------------ Render match ---------------- + DC_adjust = cand_pixels_final(num_cand_pixels_final).DC; + i_r = cand_pixels_final(num_cand_pixels_final).level; + j_r = cand_pixels_final(num_cand_pixels_final).res; + x_r = cand_pixels_final(num_cand_pixels_final).x_cand; + y_r = cand_pixels_final(num_cand_pixels_final).y_cand; + scale_est_r = cand_pixels_final(num_cand_pixels_final).scale; + Theta_r = cand_pixels_final(num_cand_pixels_final).rot; + lib_im_num_r = cand_pixels_final(num_cand_pixels_final).lib_im; + + + % ---------- scale cand pixel position according to scale and res --------------------- + x_r = round(x_r*2^(scaling_constant_v(rendering_level)-(i_r+j_r))); + y_r = round(y_r*2^(scaling_constant_v(rendering_level)-(i_r+j_r))); + + % ---------- read in candidate library image ------------------ + Im_lib = imread([files(lib_im_num_r).name]); % read in library image + + perc_DC = 1; % percentage DC adjusted +% DC_adjust = (DC_v(nearest_neighbor_int_index(z))-DC_component)*perc_DC; + Im_lib = Im_lib-DC_adjust; + + % -------resize library image for level ------------- + Im_lib_tmp = imresize(Im_lib, [lib_im_v(i_r) lib_im_v(i_r)]); + + + % ------- crop analysis window ------------------------- + xmin = x_r - lib_im_v(i_r)/2; + ymin = y_r - lib_im_v(i_r)/2; + width = lib_im_v(i_r)-1; + height = lib_im_v(i_r)-1; + + rect = [xmin ymin width height]; + Im_target_tile = imcrop(Render_target,rect); + + Im_target_tile_masked = zeros([lib_im_v(i_r) lib_im_v(i_r)]); + mask = zeros([lib_im_v(i_r) lib_im_v(i_r)]); + + % ---------- matrix for backing up Render target pixels ----------------- + Render_target_backup = zeros([lib_im_v(i_r) lib_im_v(i_r)]); + + % --------- scale image ---------------- + + I_sc_tmp = imresize(Im_lib_tmp, 1/scale_est_r); + I_sc_tmp_size = size(I_sc_tmp); + + tmp_height = I_sc_tmp_size(1); + tmp_width = I_sc_tmp_size(2); + + radius = ((lib_im_v(i_r)/2)/scale_est_r)*(1-perc_sloped); + + fd_2 = window(i).win; + + + if color == 0 + fd_2 = imresize(fd_2,I_sc_tmp_size); + else + fd_2 = imresize(fd_2,[I_sc_tmp_size(1) I_sc_tmp_size(2)]); + end + + + if color == 0 + mask_tmp = ones(I_sc_tmp_size); + else + mask_tmp = ones([I_sc_tmp_size(1) I_sc_tmp_size(2)]); + end + + + if color == 0 + I_sc = zeros([lib_im_v(i_r) lib_im_v(i_r)]); + else + I_sc = zeros([lib_im_v(i_r) lib_im_v(i_r) 3]); + end + + + I_comb = zeros([lib_im_v(i_r) lib_im_v(i_r)]); + + if color == 0 + I_sc(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width) = I_sc_tmp; + else + I_sc(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width, :) = I_sc_tmp; + end + mask(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width) = fd_2; + + % -------Rotate image back by theta estimate -------- + + R1 = imrotate(I_sc, -Theta_r, 'bilinear', 'crop'); + mask = imrotate(mask, -Theta_r, 'bilinear', 'crop'); + + fd_2_tmp = window(i).win_fixed; + mask = mask.*fd_2_tmp; + + alpha = 0.5; + I_comb = (alpha)*double(Im_target_tile) + (1-alpha)*R1; + + % ----------------- Render match --------------------- + for u= 1:lib_im_v(i_r) + for v = 1:lib_im_v(i_r) + + Render_target_backup(u,v) = Render_target(ymin+u-1,xmin+v-1); + if (mask(u,v) > .01) + + if color == 0 + Render_target(ymin+u-1,xmin+v-1) = uint8(R1(u,v)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1)); + else + Render_target(ymin+u-1,xmin+v-1,1) = uint8(R1(u,v,1)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,1)); + Render_target(ymin+u-1,xmin+v-1,2) = uint8(R1(u,v,2)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,2)); + Render_target(ymin+u-1,xmin+v-1,3) = uint8(R1(u,v,3)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,3)); + end + + % update bitmap of valid placed pixels + Render_target_valid(ymin+u-1,xmin+v-1) = 1; + + end + + end + end + + end + + end + end + + level_txt = 1;%num2str(i); + + if (num_cand_pixels_final>0) + save(['Fractal_Mosaic_' level_txt]) + end + +end +% end + +imwrite(Render_target,'fractal_mosaic','jpg','quality',100); + + +time = toc + diff --git a/reference/upstream-matlab/fractal_mosaic_render.m b/reference/upstream-matlab/fractal_mosaic_render.m new file mode 100644 index 0000000..7ab4b65 --- /dev/null +++ b/reference/upstream-matlab/fractal_mosaic_render.m @@ -0,0 +1,223 @@ + +clear all + +load Fractal_Mosaic % load Matlab workspace with mosaic parameters + +black_bg = 0; % determines if a black background is used or the original target image. '1' = use target image as background + + +color_render = 1; % If color_render = '1', use color images, even if matching was done in black and white mode. If color ='0', use black and white images + + +perc_sloped = 0.1; + +I1 = imread(target_path); + +size_r = size(Render_target); + +I1 = imresize(I1, [size_r(1) size_r(2)]); +% I1 = imresize(I1, [8192 8192]); +Render_target = I1; +clear I1 + + +pat = '\'; +regexprep(path, pat, '/') + +if black_bg == 1 + 'overlaying on target image' +else + Render_target(:,:,:) = 0; +end + +%---------- calculate distance windows ----------- + + +for i=1:length(level_v) +% i + + perc_sloped_fixed = 0.01;%.1; % perc sloped for fixed window (one that's used to cut off clipping) + + fd = zeros(lib_im_v(i)); + fd_fixed = zeros(lib_im_v(i)); + + y_d = ones(1,lib_im_v(i)/2); + y_d_fixed = ones(1,lib_im_v(i)/2); + + x_ind = linspace(0,lib_im_v(i)/2,lib_im_v(i)/2); + for u=1:(perc_sloped*(lib_im_v(i)/2)) + y_d(u) = (1/(perc_sloped*(lib_im_v(i)/2)))*x_ind(u); + end + for u=1:(perc_sloped_fixed*(lib_im_v(i)/2)) + y_d_fixed(u) = (1/(perc_sloped_fixed*(lib_im_v(i)/2)))*x_ind(u); + end + % + % x_ind = linspace(0,lib_im_v(i)/2,lib_im_v(i)/2); %triangular window. + % y_d = (-1/(lib_im_v(i)/2))*x_ind + 1; + + fd(1:lib_im_v(i)/2) = y_d; + fd(lib_im_v(i)/2+1:lib_im_v(i)) = fliplr(y_d); + fd_2 = fd*fd'; + window(i).win = fd_2; + + fd_fixed(1:lib_im_v(i)/2) = y_d_fixed; + fd_fixed(lib_im_v(i)/2+1:lib_im_v(i)) = fliplr(y_d_fixed); + fd_2_fixed = fd_fixed*fd_fixed'; + window(i).win_fixed = fd_2_fixed; + + +% figure(2),mesh(fd_2_fixed) +% pause +end + + +for i=1:num_cand_pixels_final + + + if color_render == 0 + + color_path = '../library_images/graffiti_256_black_white'; + + else + color_path = '../library_images/graffiti_256_color'; + + end + + + image_path = cand_pixels_final(i).file_path; % read in file path + + % if file path uses '\' in path (default for PC), replace + % with '/' (works on both PC and mac) + + for j=1:length(image_path) + if (image_path(j) == '\') + image_path(j) = '/'; + end + end + + path = [color_path image_path]; + +% % -------- (for mac) convert '\' to '/' ------ +% pat = '\'; +% path = regexprep(path, pat, '/'); + + + % ---------- read in candidate library image ------------------ + + Im_lib = imread(path); % read in library image + + + DC_adjust = cand_pixels_final(i).DC; + i_r = cand_pixels_final(i).level; + j_r = cand_pixels_final(i).res; + x_r = cand_pixels_final(i).x_cand; + y_r = cand_pixels_final(i).y_cand; + scale_est_r = cand_pixels_final(i).scale; + Theta_r = cand_pixels_final(i).rot; + lib_im_num_r = cand_pixels_final(i).lib_im; + cand_pixels_final(i).valid_pix; + +% x_admin = cand_pixels_final(i).x +% y_admin = cand_pixels_final(i).y + + perc_DC = 1; % percentage DC adjusted +% DC_adjust = (DC_v(nearest_neighbor_int_index(z))-DC_component)*perc_DC; + Im_lib = Im_lib-DC_adjust; + + + % ---------- scale cand pixel position according to scale and res --------------------- + x_r = round(x_r*2^(scaling_constant_v(rendering_level)-(i_r+j_r))); + y_r = round(y_r*2^(scaling_constant_v(rendering_level)-(i_r+j_r))); + +% radius_cand = +% ((lib_im_v(i)/2)/scale_est)*(1-perc_sloped); + % ---------- read in candidate library image ------------------ +% Im_lib = imread([files(lib_im_num_r).name]); % read in library image + + % -------resize library image for level ------------- + Im_lib_tmp = imresize(Im_lib, [lib_im_v(i_r) lib_im_v(i_r)]); + + + + % ------- crop analysis window ------------------------- + xmin = x_r - lib_im_v(i_r)/2; + ymin = y_r - lib_im_v(i_r)/2; + width = lib_im_v(i_r)-1; + height = lib_im_v(i_r)-1; + + rect = [xmin ymin width height]; + Im_target_tile = imcrop(Render_target,rect); + +% Im_target_tile_masked = zeros([lib_im_v(i_r) lib_im_v(i_r)]); + mask = zeros([lib_im_v(i_r) lib_im_v(i_r)]); + + % ---------- matrix for backing up Render target pixels ----------------- + Render_target_backup = zeros([lib_im_v(i_r) lib_im_v(i_r) 3]); + + % --------- scale image ---------------- + + I_sc_tmp = imresize(Im_lib_tmp, 1/scale_est_r); + I_sc_tmp_size = size(I_sc_tmp); + + tmp_height = I_sc_tmp_size(1); + tmp_width = I_sc_tmp_size(2); + + radius = ((lib_im_v(i_r)/2)/scale_est_r)*(1-perc_sloped); + + fd_2 = window(i_r).win; + fd_2 = imresize(fd_2,[I_sc_tmp_size(1) I_sc_tmp_size(2)]); + mask_tmp = ones([I_sc_tmp_size(1) I_sc_tmp_size(2)]); + + I_sc = zeros([lib_im_v(i_r) lib_im_v(i_r) 3]); +% I_comb = zeros([lib_im_v(i_r) lib_im_v(i_r) 3]); + + I_sc(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width, :) = I_sc_tmp; + + % mask(round(lib_im_v(i)/2 - tmp_height/2)+1:round(lib_im_v(i)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i)/2 - tmp_width/2)+1:round(lib_im_v(i)/2 - tmp_width/2)+tmp_width) = mask_tmp; + mask(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width) = fd_2; + + % -------Rotate image back by theta estimate -------- + + R1 = imrotate(I_sc, -Theta_r, 'bilinear', 'crop'); + mask = imrotate(mask, -Theta_r, 'bilinear', 'crop'); + + fd_2_tmp = window(i_r).win_fixed; + mask = mask.*fd_2_tmp; + + % pause + % window(i).win_fixed + % mask_t = mask*window(i).win; +% +% alpha = 0.5; +% I_comb = (alpha)*double(Im_target_tile) + (1-alpha)*R1; +% + + % ----------------- Render match --------------------- + for u= 1:lib_im_v(i_r)%ymin:(ymin+height) + for v = 1:lib_im_v(i_r)%xmin:(xmin+width) + + Render_target_backup(u,v) = Render_target(ymin+u-1,xmin+v-1); + if (mask(u,v) > .01) + Render_target(ymin+u-1,xmin+v-1,1) = uint8(R1(u,v,1)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,1)); + Render_target(ymin+u-1,xmin+v-1,2) = uint8(R1(u,v,2)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,2)); + Render_target(ymin+u-1,xmin+v-1,3) = uint8(R1(u,v,3)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,3)); + + end + + end + end + +% figure(2),subplot(1,3,1),imshow(uint8(Render_target_backup)) +% subplot(1,3,2),imshow(Render_target(max(y-sq_center_disp(1),1):min(y+sq_center_disp(1),size_render_target(1)), max(x-sq_center_disp(2),1):min(x+sq_center_disp(2),size_render_target(2)))); +% subplot(1,3,3),imshow(Render_target) +% figure(2),imshow(Render_target) +% + figure(2),imshow(Render_target) + drawnow +% pause + +end + +imwrite(Render_target,'fractal_mosaic_graffiti_color','jpg','quality',100); + + diff --git a/reference/upstream-matlab/gfd_calc.m b/reference/upstream-matlab/gfd_calc.m new file mode 100644 index 0000000..d73a979 --- /dev/null +++ b/reference/upstream-matlab/gfd_calc.m @@ -0,0 +1,175 @@ + +clear all +close all + +color = 0; +% ----------------------------- get jpeg file paths -------------------- + +'getting image paths...' + + if color ==0 + path = ['../library_images/graffiti_256_black_white/']; + else + path = ['C:\Americano\library_images\graffiti_256_color\']; + end + + + % files = dir('C:\Image_registration\flickr_spider\graffiti_library\*.jpg'); % get file names + directories = dir(path); % get directory names + + num_files = 0; + file_count = 0; + + dir_temp = directories(3).name; % get first directory name + files = dir([path,dir_temp,'/*.jpg']); % get file names for first directory + file_count = length(files); + + % add directory path to jpeg files name + for j=1:length(files) + files(j).name = [path,dir_temp,'/',files(j).name]; + end + + for i=2:(length(directories)-2) % offset by 2 because first two entries are '.' and '..' + + + dir_temp = directories(i+2).name; % get next directory name + + files_temp = dir([path,dir_temp,'/*.jpg']); % get file names for given directory + + % add directory path to jpeg files name + for j=1:length(files_temp) + files_temp(j).name = [path,dir_temp,'/',files_temp(j).name]; + end + + + + % concatenate file paths from next directory to existing + files((file_count+1):(file_count+length(files_temp))) = files_temp; + + + file_count = file_count+length(files_temp); + % pause + end + + +% bins = 16; +% hist_x = [(255/bins)/2:255/bins:(255-(255/bins)/2)]; +% hist_matrix = zeros(length(files),bins); + +radial_freq = 6; +angular_freq = 6; +vector_len = radial_freq*angular_freq; +gfd_matrix = zeros(length(files),vector_len); +% gfd_matrix = zeros(8000,vector_len); + +size_v = [128 128]; + +avg_gfd = zeros(radial_freq,angular_freq); + +for i=1:length(files) + + i + + + Im_lib = imread([files(i).name]); % read in library image + + Im_lib = imresize(Im_lib,size_v); % resize image to 128 + +% size_Im_lib = size(Im_lib); +% +% % ------ reshape Im_lib into 1D array ------ +% Im_lib_array = reshape(Im_lib,size_Im_lib(1)*size_Im_lib(2),1); +% +% histo = hist(double(Im_lib_array),hist_x); + + Im_lib_polar = polar_conversion(Im_lib); + Im_lib_polar_FFT = abs(fftshift(fft2(Im_lib_polar))); + + DC = Im_lib_polar_FFT(size_v(1)/2+1,size_v(2)/2+1); +% energy_norm = sum(sum(Im_lib_polar_FFT)); + + % --------------- normalize coefs ------------------------ +% Im_lib_polar_FFT(size_v(1)/2+1,size_v(2)/2+1) = Im_lib_polar_FFT(size_v(1)/2+1,size_v(2)/2+1)/(size_v(1)*size_v(2)); + + Im_lib_polar_FFT = Im_lib_polar_FFT./DC; +% Im_lib_polar_FFT = Im_lib_polar_FFT./energy_norm; + + +% Im_lib_polar_FFT(size_v(1)/2+1,size_v(2)/2+1) = DC; + + + + center_gfd = [size_v(1)/2+1 size_v(2)/2+1]; + + gfd = Im_lib_polar_FFT(center_gfd(1):(center_gfd(1)+angular_freq-1),center_gfd(2):(center_gfd(2)+radial_freq-1)); + + avg_gfd = avg_gfd + gfd; + + +% figure(1),subplot(1,2,1),imshow(uint8(Im_lib)) +% subplot(1,2,2),mesh(gfd) +% +% pause + +% I1_gfd = fourier_mellin(Im_lib); + +% hu_moment = humoments(Im_lib); + gfd_matrix(i,:) = gfd(:); + + +% figure(1),subplot(2,1,1),imshow(Im_lib) +% subplot(2,1,2),imshow(Im_lib) +% subplot(3,1,3),mesh(I1_fm) + + +% figure(1),subplot(1,2,1),imshow(Im_lib) +% subplot(1,2,2),mesh(I1_fm); +% pause + +% subplot(1,2,1),imshow(Im_lib) +% subplot(1,2,2),plot(histo) + + +end +% pause + +gfd_matrix = single(gfd_matrix); +save('gfd_matrix_graffiti_128', 'gfd_matrix') +% save('gfd_matrix_SF_128', 'gfd_matrix') +% save('avg_gfd_3045', 'avg_gfd') + +% avg_gfd = gfd_matrix(1,:); +% +% for i=2:4846 +% avg_gfd = avg_gfd+ gfd_matrix(i,:); +% +% end + + + + +% for i=1:1000%length(files) +% i +% +% fm = fm_matrix(i,:); +% +% % for j=1:1000 +% +% +% [nearest_neighbor_fm neighbor_distances_fm] = kNearestNeighbors(fm_matrix(1:1000,:), fm, 1000); +% min_fm_dist = neighbor_distances_fm(2); +% max_fm_dist = neighbor_distances_fm(end); +% +% % end +% +% end + +% fm_matrix = single(fm_matrix); +% save('gfd_matrix_graffiti_128', 'gfd_matrix') +% load('hist_matrix_graffiti'); + + +% [index_vals,vector_vals,final_nodes] = kd_knn(tree,point,k,plot_stuff,node_number) + + + diff --git a/reference/upstream-matlab/gpl_license.txt b/reference/upstream-matlab/gpl_license.txt new file mode 100644 index 0000000..94a9ed0 --- /dev/null +++ b/reference/upstream-matlab/gpl_license.txt @@ -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. + + + Copyright (C) + + 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: + + Copyright (C) + 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/reference/upstream-matlab/histogram_calc.m b/reference/upstream-matlab/histogram_calc.m new file mode 100644 index 0000000..cc80ffb --- /dev/null +++ b/reference/upstream-matlab/histogram_calc.m @@ -0,0 +1,85 @@ + +clear all +close all + +% ----------------------------- get jpeg file paths -------------------- + +'getting image paths...' + +path = ['../library_images/graffiti_256_black_white/']; % path to library images + +% files = dir('C:\Image_registration\flickr_spider\graffiti_library\*.jpg'); % get file names +directories = dir(path); % get directory names + +num_files = 0; +file_count = 0; + +dir_temp = directories(3).name; % get first directory name +files = dir([path,dir_temp,'/*.jpg']); % get file names for first directory +file_count = length(files); + +% add directory path to jpeg files name +for j=1:length(files) + files(j).name = [path,dir_temp,'/',files(j).name]; +end + +for i=2:(length(directories)-2) % offset by 2 because first two entries are '.' and '..' + + + dir_temp = directories(i+2).name; % get next directory name + + files_temp = dir([path,dir_temp,'/*.jpg']); % get file names for given directory + + % add directory path to jpeg files name + for j=1:length(files_temp) + files_temp(j).name = [path,dir_temp,'/',files_temp(j).name]; + end + + + + % concatenate file paths from next directory to existing + files((file_count+1):(file_count+length(files_temp))) = files_temp; + + + file_count = file_count+length(files_temp); +% pause +end + +bins = 16; +hist_x = [(255/bins)/2:255/bins:(255-(255/bins)/2)]; +hist_matrix = zeros(length(files),bins); + + + +for i=1:length(files) + + i + + + Im_lib = imread([files(i).name]); % read in library image + + Im_lib = imresize(Im_lib,[128 128]); % resize image to 128 + + size_Im_lib = size(Im_lib); + + % ------ reshape Im_lib into 1D array ------ + Im_lib_array = reshape(Im_lib,size_Im_lib(1)*size_Im_lib(2),1); + + histo = hist(double(Im_lib_array),hist_x); + + hist_matrix(i,:) = histo; + +% subplot(1,2,1),imshow(Im_lib) +% subplot(1,2,2),plot(histo) + + +end + +save('hist_matrix_graffiti_128', 'hist_matrix') +% load('hist_matrix_graffiti'); + + + + + + diff --git a/reference/upstream-matlab/kNearestNeighbors_no_sort.m b/reference/upstream-matlab/kNearestNeighbors_no_sort.m new file mode 100755 index 0000000..cd08797 --- /dev/null +++ b/reference/upstream-matlab/kNearestNeighbors_no_sort.m @@ -0,0 +1,41 @@ +function [neighborDistances] = kNearestNeighbors_no_sort(dataMatrix, queryMatrix, k) +%-------------------------------------------------------------------------- +% Program to find the k - nearest neighbors (kNN) within a set of points. +% Distance metric used: Euclidean distance +% +% Usage: +% [neighbors distances] = kNearestNeighbors(dataMatrix, queryMatrix, k); +% dataMatrix (N x D) - N vectors with dimensionality D (within which we search for the nearest neighbors) +% queryMatrix (M x D) - M query vectors with dimensionality D +% k (1 x 1) - Number of nearest neighbors desired +% +% Example: +% a = [1 1; 2 2; 3 2; 4 4; 5 6]; +% b = [1 1; 2 1; 6 2]; +% [neighbors distances] = kNearestNeighbors(a,b,2); +% +% Output: +% neighbors = +% 1 2 +% 1 2 +% 4 3 +% +% distances = +% 0 1.4142 +% 1.0000 1.0000 +% 2.8284 3.0000 +%-------------------------------------------------------------------------- + +% neighborIds = zeros(size(queryMatrix,1),k); +% neighborDistances = neighborIds; + +neighborDistances= zeros(size(queryMatrix,1),k); + +numDataVectors = size(dataMatrix,1); +numQueryVectors = size(queryMatrix,1); +for i=1:numQueryVectors, + neighborDistances = sum((repmat(queryMatrix(i,:),numDataVectors,1)-dataMatrix).^2,2); +% [sortval sortpos] = sort(dist,'ascend'); +% neighborIds(i,:) = sortpos(1:k); +% neighborDistances(i,:) = sqrt(sortval(1:k)); +end \ No newline at end of file diff --git a/reference/upstream-matlab/polar_conversion.m b/reference/upstream-matlab/polar_conversion.m new file mode 100755 index 0000000..561d74e --- /dev/null +++ b/reference/upstream-matlab/polar_conversion.m @@ -0,0 +1,27 @@ + +function polar = polar_conversion(A) + + +A_dim = size(A); % get dimensions of A +center = A_dim/2; % get center of image coordinates + + +% ---------- calculate polar coordinates ------------- +theta = linspace(0,2*pi,A_dim(1)+1); theta(end) = []; +% d = min([A_dim(1)-center(1) center(1)-1 A_dim(2)-center(2) center(2)-1]); +rho = linspace(0, center(1)-1,A_dim(2))'; + +% convert polar coordinates to cartesian coordinates and center +xx = rho*cos(theta) + center(1); +yy = rho*sin(theta) + center(2); + + +xx = double(xx); +yy = double(yy); + + +polar=interp2(double(A),xx,yy); + + + + diff --git a/reference/upstream-matlab/std_calc.m b/reference/upstream-matlab/std_calc.m new file mode 100644 index 0000000..297c0af --- /dev/null +++ b/reference/upstream-matlab/std_calc.m @@ -0,0 +1,102 @@ + +clear all +close all + +color = 0; +% ----------------------------- get jpeg file paths -------------------- + +'getting image paths...' + + if color ==0 + path = ['../library_images/graffiti_256_black_white/']; + else + path = ['C:\Americano\library_images\graffiti_256_color\']; + end + + + % files = dir('C:\Image_registration\flickr_spider\graffiti_library\*.jpg'); % get file names + directories = dir(path); % get directory names + + num_files = 0; + file_count = 0; + + dir_temp = directories(3).name; % get first directory name + files = dir([path,dir_temp,'/*.jpg']); % get file names for first directory + file_count = length(files); + + % add directory path to jpeg files name + for j=1:length(files) + files(j).name = [path,dir_temp,'/',files(j).name]; + end + + for i=2:(length(directories)-2) % offset by 2 because first two entries are '.' and '..' + + + dir_temp = directories(i+2).name; % get next directory name + + files_temp = dir([path,dir_temp,'/*.jpg']); % get file names for given directory + + % add directory path to jpeg files name + for j=1:length(files_temp) + files_temp(j).name = [path,dir_temp,'/',files_temp(j).name]; + end + + + + % concatenate file paths from next directory to existing + files((file_count+1):(file_count+length(files_temp))) = files_temp; + + + file_count = file_count+length(files_temp); + % pause + end + +% bins = 16; +% hist_x = [(255/bins)/2:255/bins:(255-(255/bins)/2)]; +% hist_matrix = zeros(length(files),bins); + +% radial_freq = 6; +% angular_freq = 6; +% vector_len = radial_freq*angular_freq; +std_matrix = zeros(length(files),1); +% gfd_matrix = zeros(8000,vector_len); + +size_v = [128 128]; + +% avg_gfd = zeros(radial_freq,angular_freq); + +for i=1:length(files) + + i + + + Im_lib = imread([files(i).name]); % read in library image + + Im_lib = imresize(Im_lib,size_v); % resize image to 128 + + std_lib = std2(Im_lib); + + std_matrix(i) = std_lib; + +% figure(1),subplot(2,1,1),imshow(Im_lib) +% subplot(2,1,2),imshow(Im_lib) +% subplot(3,1,3),mesh(I1_fm) + + +% figure(1),subplot(1,2,1),imshow(Im_lib) +% subplot(1,2,2),mesh(I1_fm); +% pause + +% subplot(1,2,1),imshow(Im_lib) +% subplot(1,2,2),plot(histo) + + +end +% pause + +% std_matrix = single(gfd_matrix); +save('std_matrix_graffiti_128', 'std_matrix') +% save('std_matrix_SF_128', 'std_matrix') + + + diff --git a/reference/upstream-matlab/transformImage_log.m b/reference/upstream-matlab/transformImage_log.m new file mode 100755 index 0000000..cf165d4 --- /dev/null +++ b/reference/upstream-matlab/transformImage_log.m @@ -0,0 +1,82 @@ +function [r,g,b] = transformImage_log(A, Ar, Ac, Nrho, Ntheta, Method, Center, Shape) + +% Inputs: A the input image +% Nrho the desired number of rows of transformed image +% Ntheta the desired number of columns of transformed image +% Method interpolation method (nearest,bilinear,bicubic) +% Center origin of input image +% Shape output size (full,valid) +% Class storage class of A + +global rho; + +theta = linspace(0,2*pi,Ntheta+1); theta(end) = []; + +switch Shape +case 'full' + corners = [1 1;Ar 1;Ar Ac;1 Ac]; + d = max(sqrt(sum((repmat(Center(:)',4,1)-corners).^2,2))); +case 'valid' + d = min([Ac-Center(1) Center(1)-1 Ar-Center(2) Center(2)-1]); +end + +% minScale = .5; +% rho = logspace(log10(minScale),log10(d),Nrho)'; % default 'base 10' logspace - play with d to change the scale of the log axis + +% rad_lin = linspace(logb(1,1.044), logb(256,1.044),256)'; +% rho = 0.5.*(1.044.^rad_lin); + + + +minScale = 1; +% maxScale = 127; + +% maxScale = 2; +maxScale = Ar/2-1; + + +% rad_lin = linspace(log(minScale), log(maxScale),256)'; + rad_lin = linspace(log(minScale), log(maxScale),Ar)'; +% rad_lin = linspace(logb(1,1.044), logb(127,1.044),256)'; + rho = exp(rad_lin); +% rho = (1.044.^rad_lin); + +% convert polar coordinates to cartesian coordinates and center +xx = rho*cos(theta) + Center(1); +yy = rho*sin(theta) + Center(2); + +if nargout==3 + if strcmp(Method,'nearest'), % Nearest neighbor interpolation + r=interp2(A(:,:,1),xx,yy,'nearest'); + g=interp2(A(:,:,2),xx,yy,'nearest'); + b=interp2(A(:,:,3),xx,yy,'nearest'); + elseif strcmp(Method,'bilinear'), % Linear interpolation + r=interp2(A(:,:,1),xx,yy,'linear'); + g=interp2(A(:,:,2),xx,yy,'linear'); + b=interp2(A(:,:,3),xx,yy,'linear'); + elseif strcmp(Method,'bicubic'), % Cubic interpolation + r=interp2(A(:,:,1),xx,yy,'cubic'); + g=interp2(A(:,:,2),xx,yy,'cubic'); + b=interp2(A(:,:,3),xx,yy,'cubic'); + else + error(['Unknown interpolation method: ',method]); + end + % any pixels outside , pad with black + mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1); + r(mask)=0; + g(mask)=0; + b(mask)=0; +else + if strcmp(Method,'nearest'), % Nearest neighbor interpolation + r=interp2(A,xx,yy,'nearest'); + elseif strcmp(Method,'bilinear'), % Linear interpolation + r=interp2(double(A),xx,yy,'linear'); + elseif strcmp(Method,'bicubic'), % Cubic interpolation + r=interp2(double(A),xx,yy,'cubic'); + else + error(['Unknown interpolation method: ',method]); + end + % any pixels outside warp, pad with black + mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1); + r(mask)=0; +end \ No newline at end of file diff --git a/src/fractal_mosaic/__init__.py b/src/fractal_mosaic/__init__.py new file mode 100644 index 0000000..755f114 --- /dev/null +++ b/src/fractal_mosaic/__init__.py @@ -0,0 +1,9 @@ +"""Fractal mosaic engine — RST placement inspired by s-ben/Fractal-Mosaics.""" + +from __future__ import annotations + +from fractal_mosaic.engine import generate_mosaic +from fractal_mosaic.profiles import PROFILES, get_profile + +__all__ = ["generate_mosaic", "PROFILES", "get_profile"] +__version__ = "0.1.0" diff --git a/src/fractal_mosaic/cli.py b/src/fractal_mosaic/cli.py new file mode 100644 index 0000000..7d83521 --- /dev/null +++ b/src/fractal_mosaic/cli.py @@ -0,0 +1,49 @@ +"""CLI: fractal-mosaic --master X --tiles DIR --output Y --profile mittel""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +from fractal_mosaic.engine import generate_mosaic +from fractal_mosaic.profiles import PROFILES, get_profile + + +def main(argv: list[str] | None = None) -> None: + p = argparse.ArgumentParser( + prog="fractal-mosaic", + description=( + "RST collage photomosaic (Fractal-Mosaics inspired). " + "Tiles are selected/placed with rotation, scale, overlap and edge alpha — " + "never recolored." + ), + ) + p.add_argument("--master", required=True, type=Path, help="Target / master image") + p.add_argument("--tiles", required=True, type=Path, help="Directory of tile images (pool)") + p.add_argument("--output", "-o", required=True, type=Path, help="Output image path") + p.add_argument( + "--profile", + default="mittel", + choices=sorted(PROFILES.keys()), + help="Density / quality profile", + ) + p.add_argument("--quiet", action="store_true", help="Less progress output") + args = p.parse_args(argv) + + if not args.master.is_file(): + raise SystemExit(f"master not found: {args.master}") + if not args.tiles.is_dir(): + raise SystemExit(f"tiles dir not found: {args.tiles}") + + profile = get_profile(args.profile) + generate_mosaic( + args.master, + args.tiles, + args.output, + profile, + progress=not args.quiet, + ) + + +if __name__ == "__main__": + main() diff --git a/src/fractal_mosaic/engine.py b/src/fractal_mosaic/engine.py new file mode 100644 index 0000000..617d072 --- /dev/null +++ b/src/fractal_mosaic/engine.py @@ -0,0 +1,314 @@ +"""Large-first collage placement with edge alpha — no tile recoloring.""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +from numpy.typing import NDArray +from PIL import Image + +from fractal_mosaic.logpolar import estimate_rotation_scale, to_gray_f32 +from fractal_mosaic.profiles import Profile +from fractal_mosaic.retrieval import TileRecord, load_tile_library, rank_tiles + + +def _resize_master(path: Path, max_edge: int) -> Image.Image: + im = Image.open(path).convert("RGB") + w, h = im.size + long_edge = max(w, h) + if long_edge > max_edge: + scale = max_edge / float(long_edge) + im = im.resize( + (max(1, int(w * scale)), max(1, int(h * scale))), + Image.Resampling.LANCZOS, + ) + return im + + +def _edge_alpha_mask(w: int, h: int, blend: float) -> NDArray[np.float32]: + """Soft rectangular falloff (upstream perc_sloped window idea).""" + blend = float(np.clip(blend, 0.02, 0.49)) + ys = np.linspace(-1.0, 1.0, h, dtype=np.float32) + xs = np.linspace(-1.0, 1.0, w, dtype=np.float32) + xv, yv = np.meshgrid(xs, ys) + # Distance from edge in each axis + ax = np.clip((1.0 - np.abs(xv)) / blend, 0.0, 1.0) + ay = np.clip((1.0 - np.abs(yv)) / blend, 0.0, 1.0) + return (ax * ay).astype(np.float32) + + +def _prepare_tile( + tile: TileRecord, + size: int, + rotation_deg: float, + scale: float, +) -> Image.Image: + with Image.open(tile.path) as im: + im = im.convert("RGBA") + # Center-crop to square then resize — crop only, no recolor + side = min(im.size) + left = (im.size[0] - side) // 2 + top = (im.size[1] - side) // 2 + im = im.crop((left, top, left + side, top + side)) + target = max(8, int(round(size * scale))) + im = im.resize((target, target), Image.Resampling.LANCZOS) + if abs(rotation_deg) > 0.5: + im = im.rotate(rotation_deg, resample=Image.Resampling.BICUBIC, expand=True) + return im + + +def _tile_geometry( + tile_rgba: Image.Image, + cx: int, + cy: int, + canvas_shape: tuple[int, int, int], + edge_blend: float, + tile_alpha: float, +) -> tuple[slice, slice, NDArray[np.float32], NDArray[np.float32]] | None: + """Return (ys, xs, rgb, alpha) clipped to canvas, or None if empty.""" + tw, th = tile_rgba.size + x0 = cx - tw // 2 + y0 = cy - th // 2 + x1 = x0 + tw + y1 = y0 + th + H, W, _ = canvas_shape + sx0 = max(0, -x0) + sy0 = max(0, -y0) + sx1 = tw - max(0, x1 - W) + sy1 = th - max(0, y1 - H) + if sx1 <= sx0 or sy1 <= sy0: + return None + dx0, dy0 = x0 + sx0, y0 + sy0 + dx1, dy1 = x0 + sx1, y0 + sy1 + patch = np.asarray(tile_rgba.crop((sx0, sy0, sx1, sy1)), dtype=np.float32) / 255.0 + ph, pw = patch.shape[:2] + mask = _edge_alpha_mask(pw, ph, edge_blend) + a = (patch[..., 3] * mask * tile_alpha).astype(np.float32) + return slice(dy0, dy1), slice(dx0, dx1), patch[..., :3], a + + +def _overlap_fraction( + coverage: NDArray[np.float32], + ys: slice, + xs: slice, + alpha: NDArray[np.float32], +) -> float: + solid = alpha > 0.05 + if not np.any(solid): + return 1.0 + return float(coverage[ys, xs][solid].mean()) + + +def _paste_rgba( + canvas: NDArray[np.float32], + coverage: NDArray[np.float32], + tile_rgba: Image.Image, + cx: int, + cy: int, + edge_blend: float, + tile_alpha: float, +) -> float: + """Paste tile centered at (cx,cy). Returns overlap fraction before paste.""" + geo = _tile_geometry(tile_rgba, cx, cy, canvas.shape, edge_blend, tile_alpha) + if geo is None: + return 1.0 + ys, xs, rgb, a = geo + overlap = _overlap_fraction(coverage, ys, xs, a) + dest = canvas[ys, xs] + aa = a[..., None] + dest[:] = dest * (1.0 - aa) + rgb * aa + coverage[ys, xs] = np.maximum(coverage[ys, xs], a) + return overlap + + +def _rms(patch: NDArray[np.uint8], tile_rgb: NDArray[np.uint8]) -> float: + a = patch.astype(np.float32) + b = np.asarray( + Image.fromarray(tile_rgb).resize((patch.shape[1], patch.shape[0]), Image.Resampling.BILINEAR) + ).astype(np.float32) + return float(np.sqrt(np.mean((a - b) ** 2))) + + +def generate_mosaic( + master_path: Path, + tiles_dir: Path, + output_path: Path, + profile: Profile, + *, + progress: bool = True, +) -> Path: + master_path = Path(master_path) + tiles_dir = Path(tiles_dir) + output_path = Path(output_path) + + library = load_tile_library(tiles_dir) + if len(library) < profile.min_tiles: + raise SystemExit( + f"tile pool has {len(library)} images; profile {profile.name!r} needs ≥ {profile.min_tiles}" + ) + + master_im = _resize_master(master_path, profile.master_max) + master = np.asarray(master_im, dtype=np.uint8) + H, W = master.shape[:2] + + if profile.use_master_bg: + canvas = master.astype(np.float32) / 255.0 + else: + canvas = np.zeros((H, W, 3), dtype=np.float32) + coverage = np.zeros((H, W), dtype=np.float32) + + rng = np.random.default_rng(42) + placed = 0 + reuse_counts: dict[int, int] = {} + + for tile_size in profile.tile_sizes: + if placed >= profile.max_placements: + break + stride = max(8, int(tile_size * profile.stride_frac)) + # Prefer uncovered regions: sample grid, sort by low coverage + centers: list[tuple[float, int, int]] = [] + for y in range(tile_size // 2, H - tile_size // 2, stride): + for x in range(tile_size // 2, W - tile_size // 2, stride): + y0, x0 = y - tile_size // 2, x - tile_size // 2 + y1, x1 = y0 + tile_size, x0 + tile_size + cov = float(coverage[y0:y1, x0:x1].mean()) + centers.append((cov, x, y)) + centers.sort(key=lambda t: t[0]) + # shuffle lightly among similar coverage + if len(centers) > 4: + head = centers[: max(1, len(centers) // 3)] + rng.shuffle(head) + centers = head + centers[len(head) :] + + for cov0, cx, cy in centers: + if placed >= profile.max_placements: + break + if cov0 > profile.max_overlap + 0.15: + continue + + y0, x0 = cy - tile_size // 2, cx - tile_size // 2 + y1, x1 = y0 + tile_size, x0 + tile_size + patch = master[y0:y1, x0:x1] + if patch.size == 0: + continue + + idxs = rank_tiles(patch, library, top_k=profile.candidates) + best: tuple[float, Image.Image, int] | None = None + + for ti in idxs: + # Soft reuse limit: denser profiles allow more repeats + limit = 3 if profile.name == "wenig" else (6 if profile.name == "mittel" else 12) + if reuse_counts.get(ti, 0) >= limit: + continue + tile = library[ti] + # Log-polar hint for rot/scale, then snap to discrete search + try: + rot_est, scale_est = estimate_rotation_scale( + to_gray_f32(patch), + to_gray_f32(tile.thumb), + n_rho=48, + n_theta=96, + ) + except Exception: + rot_est, scale_est = 0.0, 1.0 + + rot_list = [0.0] + step = float(profile.rotation_step) + for mul in (-2, -1, 1, 2): + rot_list.append(rot_est + mul * step) + rot_list.append(rot_est) + # Normalize to [-180,180) and unique + rot_list = sorted({((r + 180) % 360) - 180 for r in rot_list}) + + scales = list(profile.scales) + if 0.5 <= scale_est <= 2.0: + scales = sorted(set(scales + [float(np.clip(scale_est, 0.7, 1.3))])) + + for rot in rot_list: + for sc in scales: + prepared = _prepare_tile(tile, tile_size, rot, sc) + # Compare center crop of prepared RGB to patch (no recolor) + rgb = prepared.convert("RGB") + sample = np.asarray( + rgb.resize((patch.shape[1], patch.shape[0]), Image.Resampling.BILINEAR), + dtype=np.uint8, + ) + err = _rms(patch, sample) + if best is None or err < best[0]: + best = (err, prepared, ti) + + if best is None: + continue + err, prepared, ti = best + # Accept threshold scaled by tile size (coarser allowed for large tiles) + thresh = 55.0 + (220 - min(tile_size, 220)) * 0.08 + if err > thresh: + continue + + geo = _tile_geometry( + prepared, cx, cy, canvas.shape, profile.edge_blend, profile.tile_alpha + ) + if geo is None: + continue + ys, xs, _rgb, a = geo + overlap = _overlap_fraction(coverage, ys, xs, a) + if overlap > profile.max_overlap: + continue + + _paste_rgba( + canvas, + coverage, + prepared, + cx, + cy, + profile.edge_blend, + profile.tile_alpha, + ) + reuse_counts[ti] = reuse_counts.get(ti, 0) + 1 + placed += 1 + if progress and placed % 10 == 0: + print(f"placed {placed}/{profile.max_placements} (tile≈{tile_size}px)", flush=True) + + # Mild residual fill: boost uncovered with smaller opaque tiles + uncovered = coverage < 0.15 + if np.any(uncovered) and placed < profile.max_placements: + ys, xs = np.where(uncovered[::8, ::8]) + for yi, xi in zip(ys[:40], xs[:40], strict=False): + if placed >= profile.max_placements: + break + cy, cx = int(yi * 8), int(xi * 8) + if not (0 <= cy < H and 0 <= cx < W): + continue + ts = profile.tile_sizes[-1] + y0, x0 = max(0, cy - ts // 2), max(0, cx - ts // 2) + y1, x1 = min(H, y0 + ts), min(W, x0 + ts) + patch = master[y0:y1, x0:x1] + if patch.size < 16: + continue + idxs = rank_tiles(patch, library, top_k=3) + if not idxs: + continue + ti = idxs[0] + prepared = _prepare_tile(library[ti], ts, float(rng.uniform(-20, 20)), 1.0) + _paste_rgba( + canvas, + coverage, + prepared, + cx, + cy, + profile.edge_blend, + min(1.0, profile.tile_alpha + 0.05), + ) + placed += 1 + + out = (np.clip(canvas, 0.0, 1.0) * 255.0).astype(np.uint8) + output_path.parent.mkdir(parents=True, exist_ok=True) + Image.fromarray(out, mode="RGB").save(output_path, optimize=True) + if progress: + fill = float(coverage.mean()) + print( + f"done: {placed} placements, coverage≈{fill:.0%}, wrote {output_path}", + flush=True, + ) + return output_path diff --git a/src/fractal_mosaic/logpolar.py b/src/fractal_mosaic/logpolar.py new file mode 100644 index 0000000..ddc87f7 --- /dev/null +++ b/src/fractal_mosaic/logpolar.py @@ -0,0 +1,85 @@ +"""Log-polar transform + FFT phase correlation for rotation/scale (Fractal-Mosaics core idea).""" + +from __future__ import annotations + +import numpy as np +from numpy.typing import NDArray +from scipy import ndimage + + +def to_gray_f32(rgb: NDArray[np.floating] | NDArray[np.integer]) -> NDArray[np.float32]: + arr = np.asarray(rgb, dtype=np.float32) + if arr.ndim == 2: + return arr + if arr.shape[2] >= 3: + r, g, b = arr[..., 0], arr[..., 1], arr[..., 2] + return 0.2989 * r + 0.5870 * g + 0.1140 * b + return arr[..., 0] + + +def log_polar( + image: NDArray[np.floating], + *, + n_rho: int | None = None, + n_theta: int | None = None, +) -> NDArray[np.float32]: + """Map Cartesian image to log-polar (rows=rho, cols=theta), matching upstream intent.""" + gray = to_gray_f32(image) + h, w = gray.shape + cy, cx = (h - 1) / 2.0, (w - 1) / 2.0 + n_rho = n_rho or max(32, min(h, w) // 2) + n_theta = n_theta or max(64, min(h, w)) + max_rho = max(1.0, min(cx, cy, w - 1 - cx, h - 1 - cy)) + # log spacing like transformImage_log.m (exp of linspace log) + rad = np.exp(np.linspace(np.log(1.0), np.log(max_rho), n_rho, dtype=np.float64)) + theta = np.linspace(0.0, 2.0 * np.pi, n_theta, endpoint=False, dtype=np.float64) + tt, rr = np.meshgrid(theta, rad) + xx = rr * np.cos(tt) + cx + yy = rr * np.sin(tt) + cy + sampled = ndimage.map_coordinates( + gray, [yy, xx], order=1, mode="constant", cval=0.0 + ) + return sampled.astype(np.float32) + + +def estimate_rotation_scale( + reference: NDArray[np.floating], + moving: NDArray[np.floating], + *, + n_rho: int = 64, + n_theta: int = 128, +) -> tuple[float, float]: + """ + Estimate rotation (degrees) and scale (moving→reference) via log-polar phase correlation. + Returns (rotation_deg, scale). Scale ~1 means same size. + """ + a = log_polar(reference, n_rho=n_rho, n_theta=n_theta) + b = log_polar(moving, n_rho=n_rho, n_theta=n_theta) + # Normalize + a = a - a.mean() + b = b - b.mean() + fa = np.fft.fft2(a) + fb = np.fft.fft2(b) + cross = fa * np.conj(fb) + denom = np.abs(cross) + denom[denom < 1e-8] = 1e-8 + cross /= denom + corr = np.fft.ifft2(cross).real + peak = np.unravel_index(int(np.argmax(corr)), corr.shape) + d_rho, d_theta = int(peak[0]), int(peak[1]) + # Wrap shifts + if d_theta > n_theta // 2: + d_theta -= n_theta + if d_rho > n_rho // 2: + d_rho -= n_rho + rotation_deg = (d_theta / float(n_theta)) * 360.0 + # Scale from rho shift: log-polar row index ≈ log(rho) + # Approximation: scale = exp(d_rho * log(max/min) / n_rho) + # With rho from 1..max_rho: delta_log ≈ d_rho * (log(max)-log(1))/(n_rho-1) + max_rho = float(min(reference.shape[:2]) / 2.0 - 1.0) + max_rho = max(max_rho, 2.0) + delta_log = d_rho * (np.log(max_rho) - np.log(1.0)) / max(n_rho - 1, 1) + scale = float(np.exp(delta_log)) + # Clamp absurd scales from weak correlations + scale = float(np.clip(scale, 0.4, 2.5)) + return float(rotation_deg), scale diff --git a/src/fractal_mosaic/profiles.py b/src/fractal_mosaic/profiles.py new file mode 100644 index 0000000..d157789 --- /dev/null +++ b/src/fractal_mosaic/profiles.py @@ -0,0 +1,92 @@ +"""Named profiles mapped from upstream Fractal-Mosaics knobs (coarse for inbox).""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Mapping + + +@dataclass(frozen=True) +class Profile: + name: str + # Master long-edge cap (pixels) + master_max: int + # Tile base sizes tried (largest first) + tile_sizes: tuple[int, ...] + # How many placements to attempt (approx) + max_placements: int + # CBIR: how many library candidates per site + candidates: int + # Rotation search step in degrees (smaller = slower/better) + rotation_step: int + # Max fractional overlap of coverage map before rejecting a site + max_overlap: float + # Edge fade as fraction of half-width (upstream perc_sloped ~0.1–0.25) + edge_blend: float + # Global multiply on tile alpha (1.0 = full) + tile_alpha: float + # Minimum tiles required in pool + min_tiles: int + # Grid stride for candidate centers as fraction of tile size + stride_frac: float + # Whether to composite over master (True) or black (False) + use_master_bg: bool + # Scale factors relative to tile_size tried per candidate + scales: tuple[float, ...] + + +PROFILES: Mapping[str, Profile] = { + "wenig": Profile( + name="wenig", + master_max=960, + tile_sizes=(220, 160, 110), + max_placements=80, + candidates=8, + rotation_step=30, + max_overlap=0.45, + edge_blend=0.22, + tile_alpha=0.92, + min_tiles=8, + stride_frac=0.55, + use_master_bg=True, + scales=(1.0, 0.85), + ), + "mittel": Profile( + name="mittel", + master_max=1280, + tile_sizes=(200, 140, 100, 72), + max_placements=220, + candidates=12, + rotation_step=20, + max_overlap=0.35, + edge_blend=0.18, + tile_alpha=0.90, + min_tiles=20, + stride_frac=0.45, + use_master_bg=True, + scales=(1.1, 1.0, 0.85), + ), + "dicht": Profile( + name="dicht", + master_max=1600, + tile_sizes=(180, 130, 90, 64, 48), + max_placements=480, + candidates=16, + rotation_step=15, + max_overlap=0.28, + edge_blend=0.14, + tile_alpha=0.88, + min_tiles=40, + stride_frac=0.35, + use_master_bg=True, + scales=(1.15, 1.0, 0.85, 0.7), + ), +} + + +def get_profile(name: str) -> Profile: + key = name.strip().lower() + if key not in PROFILES: + known = ", ".join(sorted(PROFILES)) + raise SystemExit(f"unknown profile {name!r} (expected one of: {known})") + return PROFILES[key] diff --git a/src/fractal_mosaic/retrieval.py b/src/fractal_mosaic/retrieval.py new file mode 100644 index 0000000..953c2b9 --- /dev/null +++ b/src/fractal_mosaic/retrieval.py @@ -0,0 +1,93 @@ +"""CBIR-style prefilter (HSV / mean colour) — selection only, never recolors tiles.""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +from numpy.typing import NDArray +from PIL import Image + + +IMAGE_SUFFIXES = {".png", ".jpg", ".jpeg", ".webp", ".bmp", ".tif", ".tiff"} + + +@dataclass +class TileRecord: + path: Path + mean_rgb: NDArray[np.float32] + hist: NDArray[np.float32] + # cached thumb for matching (RGB uint8) + thumb: NDArray[np.uint8] + + +def _hist_hsv(rgb: NDArray[np.uint8], bins: int = 8) -> NDArray[np.float32]: + """Flattened 3D HSV histogram (normalized).""" + arr = rgb.astype(np.float32) / 255.0 + r, g, b = arr[..., 0], arr[..., 1], arr[..., 2] + mx = np.maximum(np.maximum(r, g), b) + mn = np.minimum(np.minimum(r, g), b) + diff = mx - mn + s = np.where(mx > 1e-6, diff / mx, 0.0) + v = mx + h = np.zeros_like(mx) + mask = diff > 1e-6 + # Hue in [0,1) + rc = ((g - b) / (diff + 1e-8)) % 6 + gc = (b - r) / (diff + 1e-8) + 2 + bc = (r - g) / (diff + 1e-8) + 4 + h = np.where(mask & (mx == r), rc, h) + h = np.where(mask & (mx == g), gc, h) + h = np.where(mask & (mx == b), bc, h) + h = (h / 6.0) % 1.0 + hist, _ = np.histogramdd( + np.stack([h.ravel(), s.ravel(), v.ravel()], axis=1), + bins=bins, + range=((0, 1), (0, 1), (0, 1)), + ) + hist = hist.astype(np.float32) + hist /= max(float(hist.sum()), 1e-6) + return hist.ravel() + + +def load_tile_library(tile_dir: Path, thumb_size: int = 64) -> list[TileRecord]: + paths = sorted( + p + for p in tile_dir.rglob("*") + if p.is_file() and p.suffix.lower() in IMAGE_SUFFIXES and not p.name.startswith(".") + ) + out: list[TileRecord] = [] + for p in paths: + try: + with Image.open(p) as im: + im = im.convert("RGB") + thumb_im = im.copy() + thumb_im.thumbnail((thumb_size, thumb_size), Image.Resampling.LANCZOS) + thumb = np.asarray(thumb_im, dtype=np.uint8) + except OSError: + continue + mean = thumb.reshape(-1, 3).mean(axis=0).astype(np.float32) + hist = _hist_hsv(thumb) + out.append(TileRecord(path=p, mean_rgb=mean, hist=hist, thumb=thumb)) + return out + + +def rank_tiles( + patch_rgb: NDArray[np.uint8], + library: list[TileRecord], + *, + top_k: int, +) -> list[int]: + """Return indices of best matching tiles (histogram L2 + mean colour).""" + if not library: + return [] + ph = _hist_hsv(patch_rgb) + pm = patch_rgb.reshape(-1, 3).mean(axis=0).astype(np.float32) + scores = [] + for i, t in enumerate(library): + d_h = float(np.linalg.norm(ph - t.hist)) + d_m = float(np.linalg.norm(pm - t.mean_rgb) / 255.0) + scores.append((d_h + 0.35 * d_m, i)) + scores.sort(key=lambda x: x[0]) + return [i for _, i in scores[:top_k]]