diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..29d2654 --- /dev/null +++ b/.flake8 @@ -0,0 +1,61 @@ +[flake8] +select = B,C,E,F,P,W,B9 +max-line-length = 100 +### DEFAULT IGNORES FOR 4-space INDENTED PROJECTS ### +# Main Explanation Docs: https://lintlyci.github.io/Flake8Rules/ +# +# E127, E128 are hard to silence in certain nested formatting situations. +# E203 doesn't work for slicing +# E265, E266 talk about comment formatting which is too opinionated. +# E402 warns on imports coming after statements. There are important use cases +# like demandimport (https://fburl.com/demandimport) that require statements +# before imports. +# E501 is not flexible enough, we're using B950 instead. +# E722 is a duplicate of B001. +# F811 looks for duplicate imports + noise for overload typing +# P207 is a duplicate of B003. +# P208 is a duplicate of C403. +# W503 talks about operator formatting which is too opinionated. +ignore = E127, E128, E203, E265, E266, E402, E501, E722, F811, P207, P208, W503 +### DEFAULT IGNORES FOR 2-space INDENTED PROJECTS (uncomment) ### +# ignore = E111, E114, E121, E127, E128, E265, E266, E402, E501, P207, P208, W503 +exclude = + .git, + .hg, + __pycache__, + _bin/*, + _build/*, + _ig_fbcode_wheel/*, + buck-out/*, + third-party-buck/*, + third-party2/* + +# Calculate max-complexity by changing the value below to 1, then surveying fbcode +# to see the distribution of complexity: +# find ./[a-z0-9]* -name 'buck-*' -prune -o -name 'third*party*' -prune -o \ +# -name '*.py' -print |\ +# parallel flake8 --config ./.flake8 |\ +# perl -ne 'if (/C901/) { s/.*\((\d+)\)/$1/; print; }' | stats +# NOTE: This will take a while to run (near an hour IME) so you probably want a +# second working dir to run it in. +# Pick a reasonable point from there (e.g. p95 or "95%") +# As of 2016-05-18 the rough distribution is: +# +# count: 134807 +# min: 2 +# max: 206 +# avg: 4.361 +# median: 3 +# sum: 587882 +# stddev: 4.317 +# variance: 18.635 +# +# percentiles: +# 75%: 5 +# 90%: 8 +# 95%: 11 +# 99%: 20 +# 99.9%: 48 +# 99.99%: 107 +# 99.999%: 160 +max-complexity = 12 diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md new file mode 100644 index 0000000..3232ed6 --- /dev/null +++ b/CODE_OF_CONDUCT.md @@ -0,0 +1,80 @@ +# Code of Conduct + +## Our Pledge + +In the interest of fostering an open and welcoming environment, we as +contributors and maintainers pledge to make participation in our project and +our community a harassment-free experience for everyone, regardless of age, body +size, disability, ethnicity, sex characteristics, gender identity and expression, +level of experience, education, socio-economic status, nationality, personal +appearance, race, religion, or sexual identity and orientation. + +## Our Standards + +Examples of behavior that contributes to creating a positive environment +include: + +* Using welcoming and inclusive language +* Being respectful of differing viewpoints and experiences +* Gracefully accepting constructive criticism +* Focusing on what is best for the community +* Showing empathy towards other community members + +Examples of unacceptable behavior by participants include: + +* The use of sexualized language or imagery and unwelcome sexual attention or +advances +* Trolling, insulting/derogatory comments, and personal or political attacks +* Public or private harassment +* Publishing others' private information, such as a physical or electronic +address, without explicit permission +* Other conduct which could reasonably be considered inappropriate in a +professional setting + +## Our Responsibilities + +Project maintainers are responsible for clarifying the standards of acceptable +behavior and are expected to take appropriate and fair corrective action in +response to any instances of unacceptable behavior. + +Project maintainers have the right and responsibility to remove, edit, or +reject comments, commits, code, wiki edits, issues, and other contributions +that are not aligned to this Code of Conduct, or to ban temporarily or +permanently any contributor for other behaviors that they deem inappropriate, +threatening, offensive, or harmful. + +## Scope + +This Code of Conduct applies within all project spaces, and it also applies when +an individual is representing the project or its community in public spaces. +Examples of representing a project or community include using an official +project e-mail address, posting via an official social media account, or acting +as an appointed representative at an online or offline event. Representation of +a project may be further defined and clarified by project maintainers. + +This Code of Conduct also applies outside the project spaces when there is a +reasonable belief that an individual's behavior may have a negative impact on +the project or its community. + +## Enforcement + +Instances of abusive, harassing, or otherwise unacceptable behavior may be +reported by contacting the project team at . All +complaints will be reviewed and investigated and will result in a response that +is deemed necessary and appropriate to the circumstances. The project team is +obligated to maintain confidentiality with regard to the reporter of an incident. +Further details of specific enforcement policies may be posted separately. + +Project maintainers who do not follow or enforce the Code of Conduct in good +faith may face temporary or permanent repercussions as determined by other +members of the project's leadership. + +## Attribution + +This Code of Conduct is adapted from the [Contributor Covenant][homepage], version 1.4, +available at https://www.contributor-covenant.org/version/1/4/code-of-conduct.html + +[homepage]: https://www.contributor-covenant.org + +For answers to common questions about this code of conduct, see +https://www.contributor-covenant.org/faq diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..ba48590 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,34 @@ +# Contributing to Meta Open Source Projects + +We want to make contributing to this project as easy and transparent as +possible. + +## Pull Requests +We actively welcome your pull requests. + +Note: pull requests are not imported into the GitHub directory in the usual way. There is an internal Meta repository that is the "source of truth" for the project. The GitHub repository is generated *from* the internal Meta repository. So we don't merge GitHub PRs directly to the GitHub repository -- they must first be imported into internal Meta repository. When Meta employees look at the GitHub PR, there is a special button visible only to them that executes that import. The changes are then automatically reflected from the internal Meta repository back to GitHub. This is why you won't see your PR having being directly merged, but you still see your changes in the repository once it reflects the imported changes. + +1. Fork the repo and create your branch from `main`. +2. If you've added code that should be tested, add tests. +3. If you've changed APIs, update the documentation. +4. Ensure the test suite passes. +5. Make sure your code lints. +6. If you haven't already, complete the Contributor License Agreement ("CLA"). + +## Contributor License Agreement ("CLA") +In order to accept your pull request, we need you to submit a CLA. You only need +to do this once to work on any of Meta's open source projects. + +Complete your CLA here: + +## Issues +We use GitHub issues to track public bugs. Please ensure your description is +clear and has sufficient instructions to be able to reproduce the issue. + +Meta has a [bounty program](https://www.facebook.com/whitehat/) for the safe +disclosure of security bugs. In those cases, please go through the process +outlined on that page and do not file a public issue. + +## License +By contributing to this project, you agree that your contributions will be licensed +under the LICENSE file in the root directory of this source tree. diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d1bbe80 --- /dev/null +++ b/LICENSE @@ -0,0 +1,400 @@ + +Attribution-NonCommercial 4.0 International + +======================================================================= + +Creative Commons Corporation ("Creative Commons") is not a law firm and +does not provide legal services or legal advice. Distribution of +Creative Commons public licenses does not create a lawyer-client or +other relationship. Creative Commons makes its licenses and related +information available on an "as-is" basis. Creative Commons gives no +warranties regarding its licenses, any material licensed under their +terms and conditions, or any related information. Creative Commons +disclaims all liability for damages resulting from their use to the +fullest extent possible. + +Using Creative Commons Public Licenses + +Creative Commons public licenses provide a standard set of terms and +conditions that creators and other rights holders may use to share +original works of authorship and other material subject to copyright +and certain other rights specified in the public license below. The +following considerations are for informational purposes only, are not +exhaustive, and do not form part of our licenses. + + Considerations for licensors: Our public licenses are + intended for use by those authorized to give the public + permission to use material in ways otherwise restricted by + copyright and certain other rights. Our licenses are + irrevocable. Licensors should read and understand the terms + and conditions of the license they choose before applying it. + Licensors should also secure all rights necessary before + applying our licenses so that the public can reuse the + material as expected. Licensors should clearly mark any + material not subject to the license. This includes other CC- + licensed material, or material used under an exception or + limitation to copyright. More considerations for licensors: + wiki.creativecommons.org/Considerations_for_licensors + + Considerations for the public: By using one of our public + licenses, a licensor grants the public permission to use the + licensed material under specified terms and conditions. If + the licensor's permission is not necessary for any reason--for + example, because of any applicable exception or limitation to + copyright--then that use is not regulated by the license. Our + licenses grant only permissions under copyright and certain + other rights that a licensor has authority to grant. Use of + the licensed material may still be restricted for other + reasons, including because others have copyright or other + rights in the material. A licensor may make special requests, + such as asking that all changes be marked or described. + Although not required by our licenses, you are encouraged to + respect those requests where reasonable. More_considerations + for the public: + wiki.creativecommons.org/Considerations_for_licensees + +======================================================================= + +Creative Commons Attribution-NonCommercial 4.0 International Public +License + +By exercising the Licensed Rights (defined below), You accept and agree +to be bound by the terms and conditions of this Creative Commons +Attribution-NonCommercial 4.0 International Public License ("Public +License"). To the extent this Public License may be interpreted as a +contract, You are granted the Licensed Rights in consideration of Your +acceptance of these terms and conditions, and the Licensor grants You +such rights in consideration of benefits the Licensor receives from +making the Licensed Material available under these terms and +conditions. + +Section 1 -- Definitions. + + a. Adapted Material means material subject to Copyright and Similar + Rights that is derived from or based upon the Licensed Material + and in which the Licensed Material is translated, altered, + arranged, transformed, or otherwise modified in a manner requiring + permission under the Copyright and Similar Rights held by the + Licensor. For purposes of this Public License, where the Licensed + Material is a musical work, performance, or sound recording, + Adapted Material is always produced where the Licensed Material is + synched in timed relation with a moving image. + + b. Adapter's License means the license You apply to Your Copyright + and Similar Rights in Your contributions to Adapted Material in + accordance with the terms and conditions of this Public License. + + c. Copyright and Similar Rights means copyright and/or similar rights + closely related to copyright including, without limitation, + performance, broadcast, sound recording, and Sui Generis Database + Rights, without regard to how the rights are labeled or + categorized. For purposes of this Public License, the rights + specified in Section 2(b)(1)-(2) are not Copyright and Similar + Rights. + d. Effective Technological Measures means those measures that, in the + absence of proper authority, may not be circumvented under laws + fulfilling obligations under Article 11 of the WIPO Copyright + Treaty adopted on December 20, 1996, and/or similar international + agreements. + + e. Exceptions and Limitations means fair use, fair dealing, and/or + any other exception or limitation to Copyright and Similar Rights + that applies to Your use of the Licensed Material. + + f. Licensed Material means the artistic or literary work, database, + or other material to which the Licensor applied this Public + License. + + g. Licensed Rights means the rights granted to You subject to the + terms and conditions of this Public License, which are limited to + all Copyright and Similar Rights that apply to Your use of the + Licensed Material and that the Licensor has authority to license. + + h. Licensor means the individual(s) or entity(ies) granting rights + under this Public License. + + i. NonCommercial means not primarily intended for or directed towards + commercial advantage or monetary compensation. For purposes of + this Public License, the exchange of the Licensed Material for + other material subject to Copyright and Similar Rights by digital + file-sharing or similar means is NonCommercial provided there is + no payment of monetary compensation in connection with the + exchange. + + j. Share means to provide material to the public by any means or + process that requires permission under the Licensed Rights, such + as reproduction, public display, public performance, distribution, + dissemination, communication, or importation, and to make material + available to the public including in ways that members of the + public may access the material from a place and at a time + individually chosen by them. + + k. Sui Generis Database Rights means rights other than copyright + resulting from Directive 96/9/EC of the European Parliament and of + the Council of 11 March 1996 on the legal protection of databases, + as amended and/or succeeded, as well as other essentially + equivalent rights anywhere in the world. + + l. You means the individual or entity exercising the Licensed Rights + under this Public License. Your has a corresponding meaning. + +Section 2 -- Scope. + + a. License grant. + + 1. Subject to the terms and conditions of this Public License, + the Licensor hereby grants You a worldwide, royalty-free, + non-sublicensable, non-exclusive, irrevocable license to + exercise the Licensed Rights in the Licensed Material to: + + a. reproduce and Share the Licensed Material, in whole or + in part, for NonCommercial purposes only; and + + b. produce, reproduce, and Share Adapted Material for + NonCommercial purposes only. + + 2. Exceptions and Limitations. For the avoidance of doubt, where + Exceptions and Limitations apply to Your use, this Public + License does not apply, and You do not need to comply with + its terms and conditions. + + 3. Term. The term of this Public License is specified in Section + 6(a). + + 4. Media and formats; technical modifications allowed. The + Licensor authorizes You to exercise the Licensed Rights in + all media and formats whether now known or hereafter created, + and to make technical modifications necessary to do so. The + Licensor waives and/or agrees not to assert any right or + authority to forbid You from making technical modifications + necessary to exercise the Licensed Rights, including + technical modifications necessary to circumvent Effective + Technological Measures. For purposes of this Public License, + simply making modifications authorized by this Section 2(a) + (4) never produces Adapted Material. + + 5. Downstream recipients. + + a. Offer from the Licensor -- Licensed Material. Every + recipient of the Licensed Material automatically + receives an offer from the Licensor to exercise the + Licensed Rights under the terms and conditions of this + Public License. + + b. No downstream restrictions. You may not offer or impose + any additional or different terms or conditions on, or + apply any Effective Technological Measures to, the + Licensed Material if doing so restricts exercise of the + Licensed Rights by any recipient of the Licensed + Material. + + 6. No endorsement. Nothing in this Public License constitutes or + may be construed as permission to assert or imply that You + are, or that Your use of the Licensed Material is, connected + with, or sponsored, endorsed, or granted official status by, + the Licensor or others designated to receive attribution as + provided in Section 3(a)(1)(A)(i). + + b. Other rights. + + 1. Moral rights, such as the right of integrity, are not + licensed under this Public License, nor are publicity, + privacy, and/or other similar personality rights; however, to + the extent possible, the Licensor waives and/or agrees not to + assert any such rights held by the Licensor to the limited + extent necessary to allow You to exercise the Licensed + Rights, but not otherwise. + + 2. Patent and trademark rights are not licensed under this + Public License. + + 3. To the extent possible, the Licensor waives any right to + collect royalties from You for the exercise of the Licensed + Rights, whether directly or through a collecting society + under any voluntary or waivable statutory or compulsory + licensing scheme. In all other cases the Licensor expressly + reserves any right to collect such royalties, including when + the Licensed Material is used other than for NonCommercial + purposes. + +Section 3 -- License Conditions. + +Your exercise of the Licensed Rights is expressly made subject to the +following conditions. + + a. Attribution. + + 1. If You Share the Licensed Material (including in modified + form), You must: + + a. retain the following if it is supplied by the Licensor + with the Licensed Material: + + i. identification of the creator(s) of the Licensed + Material and any others designated to receive + attribution, in any reasonable manner requested by + the Licensor (including by pseudonym if + designated); + + ii. a copyright notice; + + iii. a notice that refers to this Public License; + + iv. a notice that refers to the disclaimer of + warranties; + + v. a URI or hyperlink to the Licensed Material to the + extent reasonably practicable; + + b. indicate if You modified the Licensed Material and + retain an indication of any previous modifications; and + + c. indicate the Licensed Material is licensed under this + Public License, and include the text of, or the URI or + hyperlink to, this Public License. + + 2. You may satisfy the conditions in Section 3(a)(1) in any + reasonable manner based on the medium, means, and context in + which You Share the Licensed Material. For example, it may be + reasonable to satisfy the conditions by providing a URI or + hyperlink to a resource that includes the required + information. + + 3. If requested by the Licensor, You must remove any of the + information required by Section 3(a)(1)(A) to the extent + reasonably practicable. + + 4. If You Share Adapted Material You produce, the Adapter's + License You apply must not prevent recipients of the Adapted + Material from complying with this Public License. + +Section 4 -- Sui Generis Database Rights. + +Where the Licensed Rights include Sui Generis Database Rights that +apply to Your use of the Licensed Material: + + a. for the avoidance of doubt, Section 2(a)(1) grants You the right + to extract, reuse, reproduce, and Share all or a substantial + portion of the contents of the database for NonCommercial purposes + only; + + b. if You include all or a substantial portion of the database + contents in a database in which You have Sui Generis Database + Rights, then the database in which You have Sui Generis Database + Rights (but not its individual contents) is Adapted Material; and + + c. You must comply with the conditions in Section 3(a) if You Share + all or a substantial portion of the contents of the database. + +For the avoidance of doubt, this Section 4 supplements and does not +replace Your obligations under this Public License where the Licensed +Rights include other Copyright and Similar Rights. + +Section 5 -- Disclaimer of Warranties and Limitation of Liability. + + a. UNLESS OTHERWISE SEPARATELY UNDERTAKEN BY THE LICENSOR, TO THE + EXTENT POSSIBLE, THE LICENSOR OFFERS THE LICENSED MATERIAL AS-IS + AND AS-AVAILABLE, AND MAKES NO REPRESENTATIONS OR WARRANTIES OF + ANY KIND CONCERNING THE LICENSED MATERIAL, WHETHER EXPRESS, + IMPLIED, STATUTORY, OR OTHER. THIS INCLUDES, WITHOUT LIMITATION, + WARRANTIES OF TITLE, MERCHANTABILITY, FITNESS FOR A PARTICULAR + PURPOSE, NON-INFRINGEMENT, ABSENCE OF LATENT OR OTHER DEFECTS, + ACCURACY, OR THE PRESENCE OR ABSENCE OF ERRORS, WHETHER OR NOT + KNOWN OR DISCOVERABLE. WHERE DISCLAIMERS OF WARRANTIES ARE NOT + ALLOWED IN FULL OR IN PART, THIS DISCLAIMER MAY NOT APPLY TO YOU. + + b. TO THE EXTENT POSSIBLE, IN NO EVENT WILL THE LICENSOR BE LIABLE + TO YOU ON ANY LEGAL THEORY (INCLUDING, WITHOUT LIMITATION, + NEGLIGENCE) OR OTHERWISE FOR ANY DIRECT, SPECIAL, INDIRECT, + INCIDENTAL, CONSEQUENTIAL, PUNITIVE, EXEMPLARY, OR OTHER LOSSES, + COSTS, EXPENSES, OR DAMAGES ARISING OUT OF THIS PUBLIC LICENSE OR + USE OF THE LICENSED MATERIAL, EVEN IF THE LICENSOR HAS BEEN + ADVISED OF THE POSSIBILITY OF SUCH LOSSES, COSTS, EXPENSES, OR + DAMAGES. WHERE A LIMITATION OF LIABILITY IS NOT ALLOWED IN FULL OR + IN PART, THIS LIMITATION MAY NOT APPLY TO YOU. + + c. The disclaimer of warranties and limitation of liability provided + above shall be interpreted in a manner that, to the extent + possible, most closely approximates an absolute disclaimer and + waiver of all liability. + +Section 6 -- Term and Termination. + + a. This Public License applies for the term of the Copyright and + Similar Rights licensed here. However, if You fail to comply with + this Public License, then Your rights under this Public License + terminate automatically. + + b. Where Your right to use the Licensed Material has terminated under + Section 6(a), it reinstates: + + 1. automatically as of the date the violation is cured, provided + it is cured within 30 days of Your discovery of the + violation; or + + 2. upon express reinstatement by the Licensor. + + For the avoidance of doubt, this Section 6(b) does not affect any + right the Licensor may have to seek remedies for Your violations + of this Public License. + + c. For the avoidance of doubt, the Licensor may also offer the + Licensed Material under separate terms or conditions or stop + distributing the Licensed Material at any time; however, doing so + will not terminate this Public License. + + d. Sections 1, 5, 6, 7, and 8 survive termination of this Public + License. + +Section 7 -- Other Terms and Conditions. + + a. The Licensor shall not be bound by any additional or different + terms or conditions communicated by You unless expressly agreed. + + b. Any arrangements, understandings, or agreements regarding the + Licensed Material not stated herein are separate from and + independent of the terms and conditions of this Public License. + +Section 8 -- Interpretation. + + a. For the avoidance of doubt, this Public License does not, and + shall not be interpreted to, reduce, limit, restrict, or impose + conditions on any use of the Licensed Material that could lawfully + be made without permission under this Public License. + + b. To the extent possible, if any provision of this Public License is + deemed unenforceable, it shall be automatically reformed to the + minimum extent necessary to make it enforceable. If the provision + cannot be reformed, it shall be severed from this Public License + without affecting the enforceability of the remaining terms and + conditions. + + c. No term or condition of this Public License will be waived and no + failure to comply consented to unless expressly agreed to by the + Licensor. + + d. Nothing in this Public License constitutes or may be interpreted + as a limitation upon, or waiver of, any privileges and immunities + that apply to the Licensor or You, including from the legal + processes of any jurisdiction or authority. + +======================================================================= + +Creative Commons is not a party to its public +licenses. Notwithstanding, Creative Commons may elect to apply one of +its public licenses to material it publishes and in those instances +will be considered the “Licensor.” The text of the Creative Commons +public licenses is dedicated to the public domain under the CC0 Public +Domain Dedication. Except for the limited purpose of indicating that +material is shared under a Creative Commons public license or as +otherwise permitted by the Creative Commons policies published at +creativecommons.org/policies, Creative Commons does not authorize the +use of the trademark "Creative Commons" or any other trademark or logo +of Creative Commons without its prior written consent including, +without limitation, in connection with any unauthorized modifications +to any of its public licenses or any other arrangements, +understandings, or agreements concerning use of licensed material. For +the avoidance of doubt, this paragraph does not form part of the +public licenses. + +Creative Commons may be contacted at creativecommons.org. diff --git a/README.md b/README.md new file mode 100644 index 0000000..bf75ccd --- /dev/null +++ b/README.md @@ -0,0 +1,87 @@ +# DRTK - Differentiable Rendering Toolkit + +This package is a PyTorch library that provides functionality for differentiable rasterization. + +It consists of five main components + +* **transform** +* **rasterize** +* **render** +* **interpolate** +* **edge_grad** + +There are also optional components such as **msi** and **mipmap_grid_sampler**. New components may be added in the future. + +Typical flow looks like this: + +**transform** -> **rasterize** -> **render** -> **interpolate** -> **CUSTOM SHADING** -> **edge_grad** + +where: +- **transform**: projects the vertex positions from camera space to image space +- **rasterize**: performs rasterization, where pixels in the output image are associated with triangles +- **render**: computes depth and baricentric image +- **interpolate**: interpolates arbitrary vertex attributes +- **CUSTOM SHADING**: user implemented shading. +- **edge_grad**: special module that computes gradients for the **rasterize** step which is not differentiable on its own. For details please see [**Rasterized Edge Gradients: Handling Discontinuities Differentiably**](https://arxiv.org/abs/2405.02508) + +## Hellow Triangle +The "Hellow Triangle" with DRTK would look like this: +```python +import drtk +import torch as th +from torchvision.utils import save_image # to save images + +# create vertex buffer of shape [1 x n_vertices x 3], here for triangle `n_vertices` == 3 +v = th.as_tensor([[[0, 511, 1], [255, 0, 1], [511, 511, 1]]]).float().cuda() + +# create index buffer +vi = th.as_tensor([[0, 1, 2]]).int().cuda() + +# rasterize +index_img = drtk.rasterize(v, vi, height=512, width=512) + +# compute baricentrics +_, bary = drtk.render(v, vi, index_img) + +# we won't do shading, we'll just save the baricentrics and filter out the empty region +# which is marked with `-1` in `index_img` +img = bary * (index_img != -1) + +save_image(img, "render.png") +``` + +![hellow triangle](doc/hellow_triangle.png) + +## Dependencies +* PyTorch >= 2.1.0 + +## Building + +To build a wheel and install it: +``` +pip install wheel +python setup.py bdist_wheel +pip install dist/drtk-.whl +``` + +To build inplace, which is useful for package development: +``` +python setup.py build_ext --inplace -j 1 +``` + +## Contributing + +See the [CONTRIBUTING](CONTRIBUTING.md) file for how to help out. + +## License +DRTK is CC-BY-NC 4.0 licensed, as found in the LICENSE file. + +## Citation +``` +@article{pidhorskyi2024rasterized, + title={Rasterized Edge Gradients: Handling Discontinuities Differentiably}, + author={Pidhorskyi, Stanislav and Simon, Tomas and Schwartz, Gabriel and Wen, He and Sheikh, Yaser and Saragih, Jason}, + journal={arXiv preprint arXiv:2405.02508}, + year={2024} +} +``` diff --git a/doc/hellow_triangle.png b/doc/hellow_triangle.png new file mode 100644 index 0000000..b926276 Binary files /dev/null and b/doc/hellow_triangle.png differ diff --git a/drtk/__init__.py b/drtk/__init__.py new file mode 100644 index 0000000..08124f6 --- /dev/null +++ b/drtk/__init__.py @@ -0,0 +1,16 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from . import utils # noqa # noqa +from .edge_grad_estimator import edge_grad_estimator, edge_grad_estimator_ref # noqa +from .interpolate import interpolate, interpolate_ref # noqa +from .mipmap_grid_sample import mipmap_grid_sample, mipmap_grid_sample_ref # noqa +from .msi import msi # noqa +from .rasterize import rasterize # noqa +from .render import render, render_ref # noqa +from .transform import transform # noqa + +__version__ = "0.1.0" # noqa diff --git a/drtk/edge_grad_estimator.py b/drtk/edge_grad_estimator.py new file mode 100644 index 0000000..e0d9f9b --- /dev/null +++ b/drtk/edge_grad_estimator.py @@ -0,0 +1,346 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import Callable, Optional, Tuple + +import torch as th +import torch.nn.functional as thf +from drtk import edge_grad_ext +from drtk.interpolate import interpolate +from drtk.utils import index + + +th.ops.load_library(edge_grad_ext.__file__) + + +def edge_grad_estimator( + v_pix: th.Tensor, + vi: th.Tensor, + bary_img: th.Tensor, + img: th.Tensor, + index_img: th.Tensor, + v_pix_img_hook: Optional[Callable[[th.Tensor], None]] = None, +) -> th.Tensor: + """ + Args: + v_pix: Pixel-space vertex coordinates with preserved camera-space Z-values. + N x V x 3 + + vi: face vertex index list tensor + V x 3 + + bary_img: 3D barycentric coordinate image tensor + N x 3 x H x W + + img: The rendered image + N x C x H x W + + index_img: index image tensor + N x H x W + + v_pix_img_hook: a backward hook that will be registered to v_pix_img. Useful for examining + generated image space. Default None + + Returns: + returns the img argument unchanged. Optionally also returns computed + v_pix_img. Your loss should use the returned img, even though it is + unchanged. + + Note: + It is important to not spatially modify the rasterized image before passing it to edge_grad_estimator. + Any operation as long as it is differentiable is ok after the edge_grad_estimator. + + Examples of opeartions that can be done before edge_grad_estimator: + - Pixel-wise MLP + - Color mapping + - Color correction, gamma correction + If the operation is significantly non-linear, then it is recommended to do it before + edge_grad_estimator. All sorts of clipping and clamping (e.g. `x.clamp(min=0.0, max=1.0)`), must be + done before edge_grad_estimator. + + Examples of operations that are not allowed before edge_grad_estimator: + - Gaussian blur + - Warping, deformation + - Masking, cropping, making holes. + + Usage:: + + from drtk.renderlayer import edge_grad_estimator + + ... + + out = renderlayer(v, tex, campos, camrot, focal, princpt, + output_filters=["index_img", "render", "mask", "bary_img", "v_pix"]) + + img = out["render"] * out["mask"] + + img = edge_grad_estimator( + v_pix=out["v_pix"], + vi=rl.vi, + bary_img=out["bary_img"], + img=img, + index_img=out["index_img"] + ) + + optim.zero_grad() + image_loss = loss_func(img, img_gt) + image_loss.backward() + optim.step() + """ + + # Could use v_pix_img output from DRTK, but bary_img needs to be detached. + v_pix_img = interpolate(v_pix, vi, index_img, bary_img.detach()) + + img = th.ops.edge_grad_ext.edge_grad_estimator(v_pix, v_pix_img, vi, img, index_img) + + if v_pix_img_hook is not None: + v_pix_img.register_hook(v_pix_img_hook) + return img + + +def edge_grad_estimator_ref( + v_pix: th.Tensor, + vi: th.Tensor, + bary_img: th.Tensor, + img: th.Tensor, + index_img: th.Tensor, + v_pix_img_hook: Optional[Callable[[th.Tensor], None]] = None, +) -> th.Tensor: + """ + Python reference implementation for + :func:`drtk.edge_grad_estimator.edge_grad_estimator`. + """ + + # could use v_pix_img output from DRTK, but bary_img needs to be detached. + v_pix_img = interpolate(v_pix, vi, index_img, bary_img.detach()) + # pyre-fixme[16]: `EdgeGradEstimatorFunction` has no attribute `apply`. + img = EdgeGradEstimatorFunction.apply(v_pix, v_pix_img, vi, img, index_img) + + if v_pix_img_hook is not None: + v_pix_img.register_hook(v_pix_img_hook) + return img + + +class EdgeGradEstimatorFunction(th.autograd.Function): + @staticmethod + @th.cuda.amp.custom_fwd(cast_inputs=th.float32) + # pyre-fixme[14]: `forward` overrides method defined in `Function` inconsistently. + def forward( + ctx, + v_pix: th.Tensor, + v_pix_img: th.Tensor, + vi: th.Tensor, + img: th.Tensor, + index_img: th.Tensor, + ) -> th.Tensor: + ctx.save_for_backward(v_pix, img, index_img, vi) + return img + + @staticmethod + @th.cuda.amp.custom_bwd + # pyre-fixme[14]: `backward` overrides method defined in `Function` inconsistently. + def backward(ctx, grad_output: th.Tensor) -> Tuple[ + Optional[th.Tensor], + Optional[th.Tensor], + Optional[th.Tensor], + Optional[th.Tensor], + Optional[th.Tensor], + ]: + # early exit in case geometry is not optimized. + if not ctx.needs_input_grad[1]: + return None, None, None, grad_output, None + + v_pix, img, index_img, vi = ctx.saved_tensors + + x_grad = img[:, :, :, 1:] - img[:, :, :, :-1] + y_grad = img[:, :, 1:, :] - img[:, :, :-1, :] + + l_index = index_img[:, None, :, :-1] + r_index = index_img[:, None, :, 1:] + t_index = index_img[:, None, :-1, :] + b_index = index_img[:, None, 1:, :] + + x_mask = r_index != l_index + y_mask = b_index != t_index + + x_both_triangles = (r_index != -1) & (l_index != -1) + y_both_triangles = (b_index != -1) & (t_index != -1) + + iimg_clamped = index_img.clamp(min=0).long() + + # compute barycentric coordinates + b = v_pix.shape[0] + vi_img = index(vi, iimg_clamped, 0).long() + p0 = th.cat( + [index(v_pix[i], vi_img[i, ..., 0].data, 0)[None, ...] for i in range(b)], + dim=0, + ) + p1 = th.cat( + [index(v_pix[i], vi_img[i, ..., 1].data, 0)[None, ...] for i in range(b)], + dim=0, + ) + p2 = th.cat( + [index(v_pix[i], vi_img[i, ..., 2].data, 0)[None, ...] for i in range(b)], + dim=0, + ) + + v10 = p1 - p0 + v02 = p0 - p2 + n = th.cross(v02, v10) + + px, py = th.meshgrid( + th.arange(img.shape[-2], device=v_pix.device), + th.arange(img.shape[-1], device=v_pix.device), + ) + + def epsclamp(x: th.Tensor) -> th.Tensor: + return th.where(x < 0, x.clamp(max=-1e-8), x.clamp(min=1e-8)) + + # pyre-fixme[53]: Captured variable `n` is not annotated. + # pyre-fixme[53]: Captured variable `p0` is not annotated. + # pyre-fixme[53]: Captured variable `px` is not annotated. + # pyre-fixme[53]: Captured variable `py` is not annotated. + # pyre-fixme[53]: Captured variable `v02` is not annotated. + # pyre-fixme[53]: Captured variable `v10` is not annotated. + def check_if_point_inside_triangle(offset_x: int, offset_y: int) -> th.Tensor: + _px = px + offset_x + _py = py + offset_y + + vp0p = th.stack([p0[..., 0] - _px, p0[..., 1] - _py], dim=-1) / epsclamp( + n[..., 2:3] + ) + + bary_1 = v02[..., 0] * -vp0p[..., 1] + v02[..., 1] * vp0p[..., 0] + bary_2 = v10[..., 0] * -vp0p[..., 1] + v10[..., 1] * vp0p[..., 0] + + return ((bary_1 > 0) & (bary_2 > 0) & ((bary_1 + bary_2) < 1))[:, None] + + left_pnt_inside_right_triangle = ( + check_if_point_inside_triangle(-1, 0)[..., :, 1:] + & x_mask + & x_both_triangles + ) + right_pnt_inside_left_triangle = ( + check_if_point_inside_triangle(1, 0)[..., :, :-1] + & x_mask + & x_both_triangles + ) + down_pnt_inside_up_triangle = ( + check_if_point_inside_triangle(0, 1)[..., :-1, :] + & y_mask + & y_both_triangles + ) + up_pnt_inside_down_triangle = ( + check_if_point_inside_triangle(0, -1)[..., 1:, :] + & y_mask + & y_both_triangles + ) + + horizontal_intersection = ( + right_pnt_inside_left_triangle & left_pnt_inside_right_triangle + ) + vertical_intersection = ( + down_pnt_inside_up_triangle & up_pnt_inside_down_triangle + ) + + left_hangs_over_right = left_pnt_inside_right_triangle & ( + ~right_pnt_inside_left_triangle + ) + right_hangs_over_left = right_pnt_inside_left_triangle & ( + ~left_pnt_inside_right_triangle + ) + + up_hangs_over_down = up_pnt_inside_down_triangle & ( + ~down_pnt_inside_up_triangle + ) + down_hangs_over_up = down_pnt_inside_up_triangle & ( + ~up_pnt_inside_down_triangle + ) + + x_grad *= x_mask + y_grad *= y_mask + + grad_output_x = 0.5 * (grad_output[:, :, :, 1:] + grad_output[:, :, :, :-1]) + grad_output_y = 0.5 * (grad_output[:, :, 1:, :] + grad_output[:, :, :-1, :]) + + x_grad = (x_grad * grad_output_x).sum(dim=1) + y_grad = (y_grad * grad_output_y).sum(dim=1) + + x_grad_no_int = x_grad * (~horizontal_intersection[:, 0]) + y_grad_no_int = y_grad * (~vertical_intersection[:, 0]) + + x_grad_spread = th.zeros( + *x_grad_no_int.shape[:1], + x_grad_no_int.shape[1], + y_grad_no_int.shape[2], + dtype=x_grad_no_int.dtype, + device=x_grad_no_int.device, + ) + x_grad_spread[:, :, :-1] = x_grad_no_int * (~right_hangs_over_left[:, 0]) + x_grad_spread[:, :, 1:] += x_grad_no_int * (~left_hangs_over_right[:, 0]) + + y_grad_spread = th.zeros_like(x_grad_spread) + y_grad_spread[:, :-1, :] = y_grad_no_int * (~down_hangs_over_up[:, 0]) + y_grad_spread[:, 1:, :] += y_grad_no_int * (~up_hangs_over_down[:, 0]) + + # Intersections. Compute border sliding gradients + ################################################# + z_grad_spread = th.zeros_like(x_grad_spread) + x_grad_int = x_grad * horizontal_intersection[:, 0] + y_grad_int = y_grad * vertical_intersection[:, 0] + + n = thf.normalize(n, dim=-1) + n = n.permute(0, 3, 1, 2) + + n_left = n[..., :, :-1] + n_right = n[..., :, 1:] + n_up = n[..., :-1, :] + n_down = n[..., 1:, :] + + def get_dp_db(v_varying: th.Tensor, v_fixed: th.Tensor) -> th.Tensor: + """ + Computes derivative of the point position with respect to edge displacement + See drtk/src/edge_grad/edge_grad_kernel.cu + Please refer to the paper "Rasterized Edge Gradients: Handling Discontinuities Differentiably" + for details. + """ + + v_varying = thf.normalize(v_varying, dim=1) + v_fixed = thf.normalize(v_fixed, dim=1) + b = th.stack([-v_fixed[:, 1], v_fixed[:, 0]], dim=1) + b_dot_varying = (b * v_varying).sum(dim=1, keepdim=True) + return b[:, 0:1] / epsclamp(b_dot_varying) * v_varying + + # We compute partial derivatives by fixing one triangle and moving the + # other, and then vice versa. + + # Left triangle moves, right fixed + dp_dbx = get_dp_db(n_left[:, [0, 2]], -n_right[:, [0, 2]]) + x_grad_spread[:, :, :-1] += x_grad_int * dp_dbx[:, 0] + z_grad_spread[:, :, :-1] += x_grad_int * dp_dbx[:, 1] + + # Left triangle fixed, right moves + dp_dbx = get_dp_db(n_right[:, [0, 2]], n_left[:, [0, 2]]) + x_grad_spread[:, :, 1:] += x_grad_int * dp_dbx[:, 0] + z_grad_spread[:, :, 1:] += x_grad_int * dp_dbx[:, 1] + + # Upper triangle moves, lower fixed + dp_dby = get_dp_db(n_up[:, [1, 2]], -n_down[:, [1, 2]]) + y_grad_spread[:, :-1, :] += y_grad_int * dp_dby[:, 0] + z_grad_spread[:, :-1, :] += y_grad_int * dp_dby[:, 1] + + # Lower triangle moves, upper fixed + dp_dby = get_dp_db(n_down[:, [1, 2]], n_up[:, [1, 2]]) + y_grad_spread[:, 1:, :] += y_grad_int * dp_dby[:, 0] + z_grad_spread[:, 1:, :] += y_grad_int * dp_dby[:, 1] + + m = index_img == -1 + x_grad_spread[m] = 0.0 + y_grad_spread[m] = 0.0 + + grad_v_pix = -th.stack([x_grad_spread, y_grad_spread, z_grad_spread], dim=3) + + return None, grad_v_pix, None, grad_output, None diff --git a/drtk/edge_grad_ext.pyi b/drtk/edge_grad_ext.pyi new file mode 100644 index 0000000..6ce28de --- /dev/null +++ b/drtk/edge_grad_ext.pyi @@ -0,0 +1,15 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from torch import Tensor + +def edge_grad_estimator( + v_pix: Tensor, + v_pix_img: Tensor, + vi: Tensor, + img: Tensor, + index_img: Tensor, +) -> Tensor: ... diff --git a/drtk/interpolate.py b/drtk/interpolate.py new file mode 100644 index 0000000..de1c16b --- /dev/null +++ b/drtk/interpolate.py @@ -0,0 +1,106 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +import torch as th +from drtk import interpolate_ext + +th.ops.load_library(interpolate_ext.__file__) + + +def interpolate( + vert_attributes: th.Tensor, + vi: th.Tensor, + index_img: th.Tensor, + bary_img: th.Tensor, +) -> th.Tensor: + """ + Performs a linear interpolation of the vertex attributes given the barycentric coordinates + Args: + vert_attributes (th.Tensor): vertex attribute tensor + N x V x C + vi (th.Tensor): face vertex index list tensor + V x 3 + index_img (th.Tensor): index image tensor + N x H x W + bary_img (th.Tensor): 3D barycentric coordinate image tensor + N x 3 x H x W + Returns: + A tensor with interpolated vertex attributes with a shape [N, C, H, W] + Note: + 1. The default of `channels_last` is set to true to make this function backward compatible. + Please consider using the argument `channels_last` instead of permuting the result afterward. + 2. By default, the output is not contiguous. Make sure you cal .contiguous() if that is a requirement. + """ + return th.ops.interpolate_ext.interpolate(vert_attributes, vi, index_img, bary_img) + + +def interpolate_ref( + vert_attributes: th.Tensor, + vi: th.Tensor, + index_img: th.Tensor, + bary_img: th.Tensor, +) -> th.Tensor: + """ + A reference implementation for `interpolate`. See the doc string from `interpolate` + """ + + # Run reference implementation in double precision to get as good reference as possible + orig_dtype = vert_attributes.dtype + vert_attributes = vert_attributes.double() + bary_img = bary_img.double() + b = vert_attributes.shape[0] + iimg_clamped = index_img.clamp(min=0).long() + vi_img = vi[iimg_clamped].long() + + v_img = th.gather( + vert_attributes, + 1, + vi_img.view(b, -1, 1).expand(-1, -1, vert_attributes.shape[-1]), + ) + v_img = ( + v_img.view(*vi_img.shape[:3], 3, vert_attributes.shape[-1]) + .permute(0, 3, 1, 2, 4) + .contiguous() + ) + v_img = (v_img * bary_img[..., None]).sum(dim=1) + + # Do the sweep of value in the range -1..1 for the `index_img == -1` region, like + # in is done in the CUDA kernel. + undefined_region = th.stack( + [ + ( + th.arange(0, index_img.shape[-1], device=vert_attributes.device)[ + None, ... + ] + .repeat(index_img.shape[-2], 1) + .double() + * 2.0 + + 1.0 + ) + / index_img.shape[-1] + - 1.0, + ( + th.arange(0, index_img.shape[-2], device=vert_attributes.device)[ + ..., None + ] + .repeat(1, index_img.shape[-1]) + .double() + * 2.0 + + 1.0 + ) + / index_img.shape[-2] + - 1.0, + ], + dim=2, + ) + undefined_region = th.tile( + undefined_region[None], dims=[1, 1, 1, (vert_attributes.shape[-1] + 1) // 2] + )[:, :, :, : vert_attributes.shape[-1]] + v_img[index_img == -1] = undefined_region.expand(index_img.shape[0], -1, -1, -1)[ + index_img == -1, : + ] + + return v_img.permute(0, 3, 1, 2).to(orig_dtype) diff --git a/drtk/interpolate_ext.pyi b/drtk/interpolate_ext.pyi new file mode 100644 index 0000000..1824df0 --- /dev/null +++ b/drtk/interpolate_ext.pyi @@ -0,0 +1,14 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from torch import Tensor + +def interpolate( + vert_attributes: Tensor, + vi: Tensor, + index_img: Tensor, + bary_img: Tensor, +) -> Tensor: ... diff --git a/drtk/mipmap_grid_sample.py b/drtk/mipmap_grid_sample.py new file mode 100644 index 0000000..6ca2729 --- /dev/null +++ b/drtk/mipmap_grid_sample.py @@ -0,0 +1,291 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import List, Optional, Tuple + +import torch as th +import torch.nn.functional as thf +from drtk import mipmap_grid_sampler_ext + +th.ops.load_library(mipmap_grid_sampler_ext.__file__) + + +def mipmap_grid_sample( + input: List[th.Tensor], + grid: th.Tensor, + vt_dxdy_img: th.Tensor, + max_aniso: int, + mode: str = "bilinear", + padding_mode: str = "zeros", + align_corners: Optional[bool] = None, + force_max_aniso: Optional[bool] = False, + clip_grad: Optional[bool] = False, +) -> th.Tensor: + """ + Similar to torch.nn.functional.grid_sample, but supports mipmapping and anisotropic filtering which mimics + graphics hardware behaviour. + Currently, only spatial (4-D) inputs are supported. + No nearest filtering. + + Args: + input (List[th.Tensor]): A list of tensors which represents a mipmap pyramid of a texture. List should + contain highest resolution texture at index 0. + All subsequent elements of the list should contain miplayers that are twice smaller, but is not a + hard requirement. Also there is no hard requirement for all mip levels to be present. + List of tensors of shape [N x C x H_in x W_in], [N x C x H_in / 2 x W_in / 2] ... [N x C x 1 x 1] + + grid (th.Tensor): uv coordinates field according to which the inputs are sampled. + N x H_out x W_out x 2 + + vt_dxdy_img (th.Tensor): Jacobian of uv coordinates field with respect to pixel position. + N x H_out x W_out x 2 x 2 + + max_aniso (int): Maximum number of samples for anisotropic filtering. + + mode (str): Interpolation mode to calculate output values. Same as grid_sample, but without 'nearest'. + 'bilinear' | 'bicubic'. Default: 'bilinear' + + padding_mode (str): padding mode for outside grid values + 'zeros' | 'border' | 'reflection'. Default: 'zeros' + + align_corners (bool, optional): Same as in grid_sample. Default: False. + + force_max_aniso (bool, optional): Contols number of samples for anisotropic filtering. + When it is False, the extension will work similarly to the graphics hardware implementation, + e.i. it does only needed number of samples, which could be anything from 1 to max_aniso depending + on the ratio of max and min uv gradient. When force_max_aniso is True, the extension always + produces max_aniso number of samples. However this mode is only intended for + debugging/tesing/comparing with reference implementation and it is not intended for the real + usage. Default: False. + + clip_grad (bool, optional): Controls behaviour when mipmap layer is missing. + This mipmap implementation allows using not full mipmap pyramid, e.g. you can have only 2, 3 or 4 + layers for a texture of size 1024 instead of all 10. Hardware require all of the layers of the + pyramid to be present. Such relaxed requirement leads to ambiguity when it needs to sample from + the missing layer. The flag clip_grad only impacts the cases when needed layers from the mipmap + pyramid are missing. In this scenario: + - when False: it will sample from the last available layer. This will lead to aliasing and + sparsely placed taps. The downside is that it may sample from arbitrary far regions of + the texture. + - when True: it will sample from the last available layer, but it will adjust the step size + to match the sampling rate of the available layer. This will lead to aliasing but densely + placed taps. Which in turn forbids it to sample from arbitrary far regions of the texture. + + Returns: + output (Tensor): Result of sampling from inputs given the grid. + N x C x H_out x W_out + """ + + if mode != "bilinear" and mode != "bicubic": + raise ValueError( + "mipmap_grid_sample(): only 'bilinear' and 'bicubic' modes are supported " + "but got: '{}'".format(mode) + ) + if ( + padding_mode != "zeros" + and padding_mode != "border" + and padding_mode != "reflection" + ): + raise ValueError( + "mipmap_grid_sample(): expected padding_mode " + "to be 'zeros', 'border', or 'reflection', " + "but got: '{}'".format(padding_mode) + ) + if mode == "bilinear": + mode_enum = 0 + elif mode == "nearest": # not supported + mode_enum = 1 + else: # mode == 'bicubic' + mode_enum = 2 + + if padding_mode == "zeros": + padding_mode_enum = 0 + elif padding_mode == "border": + padding_mode_enum = 1 + else: # padding_mode == 'reflection' + padding_mode_enum = 2 + + if align_corners is None: + align_corners = False + + return th.ops.mipmap_grid_sampler_ext.mipmap_grid_sampler_2d( + input, + grid, + vt_dxdy_img, + max_aniso, + padding_mode_enum, + mode_enum, + align_corners, + force_max_aniso, + clip_grad, + ) + + +def mipmap_grid_sample_ref( + input: List[th.Tensor], + grid: th.Tensor, + vt_dxdy_img: th.Tensor, + max_aniso: int, + mode: str = "bilinear", + padding_mode: str = "border", + align_corners: Optional[bool] = False, + high_quality: bool = False, +) -> th.Tensor: + """ + A reference implementation for `mipmap_grid_sample`. See the doc string from `mipmap_grid_sample` + The CUDA version of `mipmap_grid_sample` should behave the same as this referense implementation when: + - `force_max_aniso` argument of `mipmap_grid_sample` is set to True + - `clip_grad` argument of `mipmap_grid_sample` is set to False + - `high_quality` argument of `mipmap_grid_sample_ref` is set to False + """ + + q = len(input) + base_level_size = list(input[0].shape[2:]) + + with th.no_grad(): + # vt_dxdy_img has assumes uv in range 0..1. + # For the comutations below we need to convert from normalized units to pixels + size = th.as_tensor( + [base_level_size[0], base_level_size[1]], + dtype=th.float32, + device=vt_dxdy_img.device, + ) + vt_dxdy_img_pixel = vt_dxdy_img * size[None, None, None, :] + + # x and y gradients magnitudes. We then need to find direction of maximum gradient and + # minimum gradients direction (principal axis) + px, py = _compute_grad_magnitude(vt_dxdy_img_pixel) + + if not high_quality: + # This is what hardware implements. + # We assume that maximum and minimum direction is either x or y. + # This assumption is a quite grude approximation + p_max = th.max(px, py) + p_min = th.min(px, py) if max_aniso != 1 else None + else: + # Instead, a more correct way would be to find principal axis using SVD + # Note this is not practical as it is very slow + u, s, v = th.linalg.svd(vt_dxdy_img_pixel) + p_max = s[..., 0] + p_min = s[..., 1] + + # Given the max and min gradients, select mipmap levels (assumes linear interpolation + # between mipmaps) + d1, a = _mipmap_selection(q, p_max, p_min, max_aniso) + + if max_aniso != 1: + if not high_quality: + uv_step_x = vt_dxdy_img[..., 0, :] + uv_step_y = vt_dxdy_img[..., 1, :] + + with th.enable_grad(): + uv_ext_x = th.cat( + [ + grid + uv_step_x * ((j + 1) / (max_aniso + 1) * 2.0 - 1.0) + for j in range(max_aniso) + ], + dim=0, + ) + uv_ext_y = th.cat( + [ + grid + uv_step_y * ((j + 1) / (max_aniso + 1) * 2.0 - 1.0) + for j in range(max_aniso) + ], + dim=0, + ) + uv_ext = th.where( + (px > py)[..., None].tile(max_aniso, 1, 1, 2), + uv_ext_x, + uv_ext_y, + ) + else: + # From SVD we have direction of the maximum gradient in the uv space. + # We ntegrate along this direction using `max_aniso` samples + uv_step = (v[..., 0, :] * s[..., 0:1]) / size[None, None, None, :] + with th.enable_grad(): + uv_ext = th.cat( + [ + grid + uv_step * ((j + 1) / (max_aniso + 1) * 2.0 - 1.0) + for j in range(max_aniso) + ], + dim=0, + ) + + result = [] + if max_aniso == 1: + for level in input: + r = thf.grid_sample( + level, + grid, + mode=mode, + padding_mode=padding_mode, + align_corners=align_corners, + ) + result.append(r) + else: + for level in input: + r = thf.grid_sample( + th.tile(level, (max_aniso, 1, 1, 1)), + uv_ext, + mode=mode, + padding_mode=padding_mode, + align_corners=align_corners, + ) + r = r.view(max_aniso, r.shape[0] // max_aniso, *r.shape[1:]).mean(dim=0) + result.append(r) + return _combine_sampled_mipmaps(result, d1, a) + + +def _mipmap_selection( + q: int, + p_max: th.Tensor, + p_min: Optional[th.Tensor], + max_aniso: int = 1, +) -> Tuple[th.Tensor, th.Tensor]: + if max_aniso != 1: + # See p.255 of OpenGL Core Profile + # N = min(ceil(Pmax/Pmin),maxAniso) + N = th.clamp(th.ceil(p_max / p_min), max=max_aniso) + N[th.isnan(N)] = 1 + + # Lambda' = log2(Pmax/N) + lambda_ = th.log2(p_max / N) + else: + lambda_ = th.log2(p_max) + + lambda_[th.isinf(lambda_)] = 0 + + # See eq. 8.15, 8.16 + # Substract small number (1e-6) so that `lambda_` is always < q - 1 + lambda_ = th.clamp(lambda_, min=0, max=q - 1 - 1e-6) + d1 = th.floor(lambda_).long() + a = lambda_ - d1.float() + return d1, a + + +def _compute_grad_magnitude(vt_dxdy_img: th.Tensor) -> Tuple[th.Tensor, th.Tensor]: + # See p.255 of OpenGL Core Profile + # Px = sqrt(dudx^2 + dvdx^2) + # Py = sqrt(dudy^2 + dvdy^2) + px = th.norm(vt_dxdy_img[..., 0, :], dim=-1) + py = th.norm(vt_dxdy_img[..., 1, :], dim=-1) + return px, py + + +def _combine_sampled_mipmaps( + sampled_mipmaps: List[th.Tensor], d1: th.Tensor, a: th.Tensor +) -> th.Tensor: + if len(sampled_mipmaps) == 1: + return sampled_mipmaps[0] + sampled_mipmaps = th.stack(sampled_mipmaps, dim=0) + indices = th.cat([d1[None, :, None], d1[None, :, None] + 1], dim=0) + samples = th.gather( + sampled_mipmaps, + dim=0, + index=indices.expand(-1, *sampled_mipmaps.shape[1:3], -1, -1), + ) + # Interpolate two nearest mipmaps. See p.266 + return th.lerp(samples[0], samples[1], a[:, None]) diff --git a/drtk/mipmap_grid_sample_ext.pyi b/drtk/mipmap_grid_sample_ext.pyi new file mode 100644 index 0000000..65100b8 --- /dev/null +++ b/drtk/mipmap_grid_sample_ext.pyi @@ -0,0 +1,21 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import List + +from torch import Tensor + +def mipmap_grid_sample_2d( + x: List[Tensor], + grid: Tensor, + vt_dxdy_img: Tensor, + max_aniso: int, + padding_mode: int, + interpolation_mode: int, + align_corners: bool, + force_max_ansio: bool, + clip_grad: bool, +) -> Tensor: ... diff --git a/drtk/msi.py b/drtk/msi.py new file mode 100644 index 0000000..95075cd --- /dev/null +++ b/drtk/msi.py @@ -0,0 +1,52 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +import torch as th +from drtk import msi_ext + +th.ops.load_library(msi_ext.__file__) + + +def msi( + ray_o: th.Tensor, + ray_d: th.Tensor, + texture: th.Tensor, + sub_step_count: int = 2, + min_inv_r: float = 1.0, + max_inv_r: float = 0.0, + stop_thresh: float = 1e-7, +) -> th.Tensor: + """ + Renders a Multi-Sphere Image which is similar to the one described in "NeRF++: Analyzing and Improving + Neural Radiance Fields" + The implementation performs bilinear sampling in the spatial dimensions of each layer and cubic between + the layers. + + Args: + ray_o (th.Tensor): Ray origins [N x 3] + + ray_d (th.Tensor): Ray directions [N x 3] + + texture (th.Tensor): The MSI texture [L x 4 x H x W], where L - number of layers. + The first 3 channels are the color channels, and the fourth one is the sigma (transmittance) + channel (negative log of alpha). + + sub_step_count (int, optional): Rate of the subsampling of the layers. Default is 2. + + min_inv_r (float, optional): Inverse of the minimum sphere radius. Default is 1 for unit radius. + + max_inv_r (float, optional): Inverse of the maximum sphere radius. Default is 0 for infinite radius. + + stop_thresh (bool, optional): The threshold for early ray termination when the accumulated + transmittance goes beyond the specified value. + + Returns: + output (Tensor): Result of the sampled MSI. First three channels are the color channels, and the 4th + one is sigma (transmittance). [N x 4] + """ + return th.ops.msi_ext.msi( + ray_o, ray_d, texture, sub_step_count, min_inv_r, max_inv_r, stop_thresh + ) diff --git a/drtk/msi_ext.pyi b/drtk/msi_ext.pyi new file mode 100644 index 0000000..d7f1163 --- /dev/null +++ b/drtk/msi_ext.pyi @@ -0,0 +1,17 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from torch import Tensor + +def msi( + ray_o: th.Tensor, + ray_d: th.Tensor, + texture: th.Tensor, + sub_step_count: int, + min_inv_r: float, + max_inv_r: float, + stop_thresh: float, +) -> th.Tensor: ... diff --git a/drtk/rasterize.py b/drtk/rasterize.py new file mode 100644 index 0000000..862a892 --- /dev/null +++ b/drtk/rasterize.py @@ -0,0 +1,91 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import Tuple + +import torch as th + +from drtk import rasterize_ext + +th.ops.load_library(rasterize_ext.__file__) + + +def rasterize( + v: th.Tensor, + vi: th.Tensor, + height: int, + width: int, + wireframe: bool = False, +) -> th.Tensor: + """ + Rasterizes a mesh defined by v and vi. + + Args: + v (th.Tensor): vertex positions. The first two components are the projected vertex's location (x, y) + on the image plane. The coordinates of the top left corner are (-0.5, -0.5), and the coordinates of + the bottom right corner are (width - 0.5, height - 0.5). The z component is expected to be in the + camera space coordinate frame (before projection). + N x V x 3 + + vi (th.Tensor): face vertex index list tensor. The most significant nibble of vi is reserved for + controlling visibility of the edges in wireframe mode. In non-wireframe mode, content of the most + significant nibble of vi will be ignored. + V x 3 + + height (int): height of the image in pixels. + + width (int): width of the image in pixels. + + wireframe (bool): If False (default), rasterizes triangles. If True, rasterizes lines, where the most + significant nibble of vi is reinterpreted as a bit field controlling the visibility of the edges. The + least significant bit controls the visibility of the first edge, the second bit controls the + visibility of the second edge, and the third bit controls the visibility of the third edge. This + limits the maximum number of vertices to 268435455. + + Returns: + The rasterized image of triangle indices which is represented with an index tensor of a shape + [N, H, W] of type int32 that stores a triangle ID for each pixel. If a triangle covers a pixel and is + the closest triangle to the camera, then the pixel will have the ID of that triangle. If no triangles + cover a pixel, then its ID is -1. + + Note: + This function is not differentiable. The gradients should be computed with `edge_grad_estimator` + instead. + """ + _, index_img = th.ops.rasterize_ext.rasterize(v, vi, height, width, wireframe) + return index_img + + +def rasterize_with_depth( + v: th.Tensor, + vi: th.Tensor, + height: int, + width: int, + wireframe: bool = False, +) -> Tuple[th.Tensor, th.Tensor]: + """ + Same as `rasterize` function, additionally returns depth image. Internally it uses the same implementation + as the rasterize function which still computes depth but does not return depth. + + Notes: + The function is not differentiable, including the depth output. + + The split is done intentionally to hide the depth image from the user as it is not differentiable which + may cause errors if assumed otherwise. Instead, the`barycentrics` function should be used instead to + compute the differentiable version of depth. + + However, we still provide `rasterize_with_depth` which returns non-differentiable depth which could allow + to avoid call to `barycentrics` function when differentiability is not required. + + Returns: + The rasterized image of triangle indices of shape [N, H, W] and a depth image of shape [N, H, W]. + Values in of pixels in the depth image not covered by any pixel are 0. + + """ + depth_img, index_img = th.ops.rasterize_ext.rasterize( + v, vi, height, width, wireframe + ) + return depth_img, index_img diff --git a/drtk/rasterize_ext.pyi b/drtk/rasterize_ext.pyi new file mode 100644 index 0000000..8d57dc6 --- /dev/null +++ b/drtk/rasterize_ext.pyi @@ -0,0 +1,16 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import List + +from torch import Tensor + +def rasterize( + v: Tensor, + vi: Tensor, + height: int, + width: int, +) -> List[Tensor]: ... diff --git a/drtk/render.py b/drtk/render.py new file mode 100644 index 0000000..4b14f51 --- /dev/null +++ b/drtk/render.py @@ -0,0 +1,107 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from functools import lru_cache +from typing import Tuple + +import torch as th + +from drtk import render_ext + +th.ops.load_library(render_ext.__file__) + + +def render( + v: th.Tensor, + vi: th.Tensor, + index_img: th.Tensor, +) -> Tuple[th.Tensor, th.Tensor]: + depth_img, bary_img = th.ops.render_ext.render(v, vi, index_img) + return depth_img, bary_img + + +def index(x, I, dim): + target_shape = [*x.shape] + del target_shape[dim] + target_shape[dim:dim] = [*I.shape] + return x.index_select(dim, I.view(-1)).reshape(target_shape) + + +@lru_cache +def _get_grid(width: int, height: int, device: th.device): + return th.stack( + th.meshgrid(th.arange(height, device=device), th.arange(width, device=device))[ + ::-1 + ], + dim=2, + ) + + +def render_ref( + v: th.Tensor, vi: th.Tensor, index_img: th.Tensor +) -> Tuple[th.Tensor, th.Tensor]: + # Run reference implementation in double precision to get as good reference as possible + + orig_dtype = v.dtype + v = v.double() + b = v.shape[0] + mask = th.ne(index_img, -1) + float_mask = mask.float()[:, None] + index_img_clamped = index_img.clamp(min=0).long() + + grid = _get_grid(index_img.shape[-1], index_img.shape[-2], device=v.device) + + # compute barycentric coordinates + vi_img = index(vi, index_img_clamped, 0).long() + v_img0 = th.cat( + [index(v[i], vi_img[i, ..., 0].data, 0)[None, ...] for i in range(b)], dim=0 + ) + v_img1 = th.cat( + [index(v[i], vi_img[i, ..., 1].data, 0)[None, ...] for i in range(b)], dim=0 + ) + v_img2 = th.cat( + [index(v[i], vi_img[i, ..., 2].data, 0)[None, ...] for i in range(b)], dim=0 + ) + + vec01 = v_img1 - v_img0 + vec02 = v_img2 - v_img0 + vec12 = v_img2 - v_img1 + + def epsclamp(x: th.Tensor) -> th.Tensor: + return th.where(x < 0, x.clamp(max=-1e-16), x.clamp(min=1e-16)) + + det = vec01[..., 0] * vec02[..., 1] - vec01[..., 1] * vec02[..., 0] + denominator = epsclamp(det) + + vp0 = grid[None, ...] - v_img0[..., :2] + vp1 = grid[None, ...] - v_img1[..., :2] + + lambda_0 = (vp1[..., 1] * vec12[..., 0] - vp1[..., 0] * vec12[..., 1]) / denominator + lambda_1 = (vp0[..., 0] * vec02[..., 1] - vp0[..., 1] * vec02[..., 0]) / denominator + lambda_2 = (vp0[..., 1] * vec01[..., 0] - vp0[..., 0] * vec01[..., 1]) / denominator + + assert th.allclose(lambda_0 + lambda_1 + lambda_2, th.ones_like(lambda_0)) + + lambda_0_mul_w0 = lambda_0 / epsclamp(v_img0[:, :, :, 2]) + lambda_1_mul_w1 = lambda_1 / epsclamp(v_img1[:, :, :, 2]) + lambda_2_mul_w2 = lambda_2 / epsclamp(v_img2[:, :, :, 2]) + zi = 1.0 / epsclamp(lambda_0_mul_w0 + lambda_1_mul_w1 + lambda_2_mul_w2) + + bary_0 = lambda_0_mul_w0 * zi + bary_1 = lambda_1_mul_w1 * zi + bary_2 = lambda_2_mul_w2 * zi + + bary_img = ( + th.cat( + (bary_0[:, None, :, :], bary_1[:, None, :, :], bary_2[:, None, :, :]), + dim=1, + ) + * float_mask + ) + + depth_img = zi * float_mask[:, 0] + + return depth_img.to(orig_dtype), bary_img.to(orig_dtype) diff --git a/drtk/render_ext.pyi b/drtk/render_ext.pyi new file mode 100644 index 0000000..cdc1bae --- /dev/null +++ b/drtk/render_ext.pyi @@ -0,0 +1,15 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import List + +from torch import Tensor + +def render( + v: Tensor, + vi: Tensor, + index_img: Tensor, +) -> List[Tensor]: ... diff --git a/drtk/screen_space_uv_derivative.py b/drtk/screen_space_uv_derivative.py new file mode 100644 index 0000000..19dc853 --- /dev/null +++ b/drtk/screen_space_uv_derivative.py @@ -0,0 +1,81 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import Optional, Sequence + +import torch as th + +from drtk.interpolate import interpolate + +from drtk.utils import face_dpdt, project_points_grad + + +def screen_space_uv_derivative( + v: th.Tensor, + vt: th.Tensor, + vi: th.Tensor, + vti: th.Tensor, + index_img: th.Tensor, + bary_img: th.Tensor, + mask: th.Tensor, + campos: th.Tensor, + camrot: th.Tensor, + focal: th.Tensor, + dist_mode: Optional[Sequence[str]] = None, + dist_coeff: Optional[th.Tensor] = None, +) -> th.Tensor: + """ + Computes per-pixel uv derivative - vt_dxdy_img with respect to the pixel-space position. + vt_dxdy_img is an image of 2x2 Jacobian matrices of the form: [[du/dx, dv/dx], + [du/dy, dv/dy]] + Shape: N x H x W x 2 x 2 + """ + + dpdt_t, vf = face_dpdt(v, vt, vi.long(), vti.long()) + + # make three of this for each vertex of the face + dpdt_t = dpdt_t[:, :, None] + dpdt_t = dpdt_t.expand(-1, -1, 3, -1, -1) + + # UV grads are not quite well interpolated with barycentrics + # We computes uv grads at per-pixel basis. For faster compute it is better to use CUDA kernel + + # make new index list, because gradients have discontinuities and we do not want to interpolate them + vi_dis = th.arange(0, 3 * vi.shape[0], dtype=th.int32, device=v.device).view(-1, 3) + dpdt_t_img = interpolate( + dpdt_t.reshape(dpdt_t.shape[0], dpdt_t.shape[1] * dpdt_t.shape[2], -1), + vi_dis, + index_img, + bary_img, + ).permute(0, 2, 3, 1) + dpdt_t_img = dpdt_t_img.view(*dpdt_t_img.shape[:3], 2, 3) + vf_img = interpolate( + vf.reshape(vf.shape[0], vf.shape[1] * vf.shape[2], -1), + vi_dis, + index_img, + bary_img, + ).permute(0, 2, 3, 1) + # duplicate vertex position vector for u and v + vf_img = vf_img[:, :, :, None].expand(-1, -1, -1, 2, -1) + # Compute 2D pixel-space gradients (d p_pix / dt)^T. + dp_pix_dt_t_img = project_points_grad( + dpdt_t_img.reshape(v.shape[0], -1, 3), + vf_img.reshape(v.shape[0], -1, 3), + campos, + camrot, + focal, + dist_mode, + dist_coeff, + ) + # Uncollapse dimension. The result is (d p_pix / dt)^T + # Where: dp_pix_dt[..., i, j] = d p_pix[j] / dt[i] + dp_pix_dt_t_img = dp_pix_dt_t_img.view(*dpdt_t_img.shape[:3], 2, 2) + # Inverse Jacobian: (dt / d p_pix)^T = ((d p_pix / dt)^T)^-1 + # Where: dt_dp_pix_t[..., i, j] = dt[j] / dp_pix[i] + vt_dxdy_img, _ = th.linalg.inv_ex(dp_pix_dt_t_img) + # pyre-fixme[16] Undefined attribute: `th.Tensor` has no attribute `__invert__`. + vt_dxdy_img[~mask, :, :] = 0 + return vt_dxdy_img diff --git a/drtk/transform.py b/drtk/transform.py new file mode 100644 index 0000000..ce4b47e --- /dev/null +++ b/drtk/transform.py @@ -0,0 +1,108 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import Optional, Sequence, Tuple + +import torch as th +from drtk.utils import project_points + + +def transform( + v: th.Tensor, + campos: Optional[th.Tensor] = None, + camrot: Optional[th.Tensor] = None, + focal: Optional[th.Tensor] = None, + princpt: Optional[th.Tensor] = None, + K: Optional[th.Tensor] = None, + Rt: Optional[th.Tensor] = None, + distortion_mode: Optional[Sequence[str]] = None, + distortion_coeff: Optional[th.Tensor] = None, + fov: Optional[th.Tensor] = None, +) -> th.Tensor: + """ + v: Tensor, N x V x 3 + Batch of vertex positions for vertices in the mesh. + + campos: Tensor, N x 3 + Camera position. + + camrot: Tensor, N x 3 x 3 + Camera rotation matrix. + + focal: Tensor, N x 2 x 2 + Focal length [[fx, 0], + [0, fy]] + + princpt: Tensor, N x 2 + Principal point [cx, cy] + + K: Tensor, N x 3 x 3 + Camera intrinsic calibration matrix. Either this or both (focal, + princpt) must be provided. + + Rt: Tensor, N x 3 x 4 or N x 4 x 4 + Camera extrinsic matrix. Either this or both (camrot, campos) must be + provided. Camrot is the upper 3x3 of Rt, campos = -R.T @ t. + + distortion_mode: List[str] + Names of the distortion modes. + + distortion_coeff: Tensor, N x 4 + Distortion coefficients. + + fov: Tensor, N x 1 + Valid field of view of the distortion model. + + """ + + v_pix, _ = transform_with_v_cam( + v, campos, camrot, focal, princpt, K, Rt, distortion_mode, distortion_coeff, fov + ) + + return v_pix + + +def transform_with_v_cam( + v: th.Tensor, + campos: Optional[th.Tensor] = None, + camrot: Optional[th.Tensor] = None, + focal: Optional[th.Tensor] = None, + princpt: Optional[th.Tensor] = None, + K: Optional[th.Tensor] = None, + Rt: Optional[th.Tensor] = None, + distortion_mode: Optional[Sequence[str]] = None, + distortion_coeff: Optional[th.Tensor] = None, + fov: Optional[th.Tensor] = None, +) -> Tuple[th.Tensor, th.Tensor]: + """ + Same as transform, but also returns the camera-space coordinates. + In most cases it is not needed, but renderlayer depends on it + """ + + if not ((camrot is not None and campos is not None) ^ (Rt is not None)): + raise ValueError("You must provide exactly one of Rt or (campos, camrot).") + + if not ((focal is not None and princpt is not None) ^ (K is not None)): + raise ValueError("You must provide exactly one of K or (focal, princpt).") + + if campos is None: + assert Rt is not None + camrot = Rt[:, :3, :3] + campos = -(camrot.transpose(-2, -1) @ Rt[:, :3, 3:4])[..., 0] + + if focal is None: + assert K is not None + focal = K[:, :2, :2] + princpt = K[:, :2, 2] + + assert camrot is not None + assert princpt is not None + # Compute camera-space 3D coordinates and 2D pixel-space projections. + v_pix, v_cam = project_points( + v, campos, camrot, focal, princpt, distortion_mode, distortion_coeff, fov + ) + + return v_pix, v_cam diff --git a/drtk/utils/__init__.py b/drtk/utils/__init__.py new file mode 100644 index 0000000..6c221ad --- /dev/null +++ b/drtk/utils/__init__.py @@ -0,0 +1,18 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from drtk.utils.geometry import ( # noqa + face_dpdt, + face_info, + vert_binormals, + vert_normals, +) +from drtk.utils.indexing import index # noqa +from drtk.utils.projection import ( # noqa + DISTORTION_MODES, # noqa + project_points, # noqa + project_points_grad, # noqa +) diff --git a/drtk/utils/geometry.py b/drtk/utils/geometry.py new file mode 100644 index 0000000..5c215a9 --- /dev/null +++ b/drtk/utils/geometry.py @@ -0,0 +1,192 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import Dict, List, Optional, Tuple, Union + +import torch as th +import torch.nn.functional as thf +from drtk.utils.indexing import index +from torch import Tensor + +eps = 1e-8 + + +def face_dpdt( + v: th.Tensor, vt: th.Tensor, vi: th.Tensor, vti: th.Tensor +) -> Tuple[th.Tensor, th.Tensor]: + """ + This function calculates the transposed Jacobian matrix (∂p/∂t)^T for each triangle. + Where: + - p represents the 3D coordinates of a point on the plane of the triangle, + - t denotes the UV coordinates assiciated with the point. + + Args: + v: vertex position tensor + N x V x 3 + + vt: vertex uv tensor + N x T x 2 + + vi: face vertex position index list tensor + F x 3 + + vti: face vertex uv index list tensor + F x 3 + + Jacobian is computed as: + + ∂p/∂t = ∂p / ∂b * (∂t / ∂b)^-1 + + Where b - barycentric coordinates + + However the implementation computes a transposed Jacobian (purely from + practical perspective - fewer permutations are needed), so the above + becomes: + + (∂p/∂t)^T = ((∂t / ∂b)^T)^-1 * (∂p / ∂b)^T + + Returns: + dpdt - transposed Jacobian (∂p/∂t)^T. Shape: N x F x 2 x 3 + Where ∂p∂t[..., i, j] = ∂p[..., j] / ∂t[..., i] + v012 - vertex positions per triangle. Shape: N x F x 3 + + where: N - batch size; F - number of triangles + """ + + v012 = v[:, vi] + vt012 = vt[:, vti] + + dpdb_t = v012[:, :, 1:3] - v012[:, :, 0:1] + dtdb_t = vt012[:, :, 1:3] - vt012[:, :, 0:1] + + # (db / dt)^T = ((dt / db)^T)^-1 + dbdt_t = th.inverse(dtdb_t) + + # (dp / dt)^T = (db / dt)^T) * (dp / db)^T + dpdt_t = dbdt_t @ dpdb_t + return dpdt_t, v012 + + +def face_attribute_to_vert(v: th.Tensor, vi: th.Tensor, attr: th.Tensor) -> Tensor: + """ + For each vertex, computes a summation of the face attributes to which the + vertex belongs. + """ + attr = ( + attr[:, :, None] + .expand(-1, -1, 3, -1) + .reshape(attr.shape[0], -1, attr.shape[-1]) + ) + vi_flat = vi.view(vi.shape[0], -1).expand(v.shape[0], -1) + vattr = th.zeros(v.shape[:-1], dtype=v.dtype, device=v.device) + + vattr = th.stack( + [vattr.scatter_add(1, vi_flat, attr[..., i]) for i in range(attr.shape[-1])], + dim=-1, + ) + return vattr + + +def face_info( + v: th.Tensor, vi: th.Tensor, to_compute: Optional[List[str]] = None +) -> Union[th.Tensor, Dict[str, th.Tensor]]: + """Given a set of vertices ``v`` and indices ``vi`` indexing into ``v`` + defining a set of faces, compute face information (normals, edges, face + areas) for each face. + + Args: + v: Vertex positions, shape [batch_size, n_vertices, 3] + + vi: Vertex indices, shape [n_faces, 3] + + to_compute: list of desired information. Any of: {normals, edges, areas}, defaults to all. + + Returns: + Dict: Face information in the following format:: + + { + "normals": shape [batch_size, n_faces, 3] + "edges": shape [batch_size, n_faces, 3, 3] + "areas": shape [batch_size, n_faces, 1] + } + + or just one of the above values not in a Dict if only one is + requested. + + """ + if to_compute is None: + to_compute = ["normals", "edges", "areas"] + + b = v.shape[0] + vi = vi.expand(b, -1, -1) + + p0 = th.stack([index(v[i], vi[i, :, 0], 0) for i in range(b)]) + p1 = th.stack([index(v[i], vi[i, :, 1], 0) for i in range(b)]) + p2 = th.stack([index(v[i], vi[i, :, 2], 0) for i in range(b)]) + v0 = p1 - p0 + v1 = p0 - p2 + + need_normals = "normals" in to_compute + need_areas = "areas" in to_compute + need_edges = "edges" in to_compute + + output = {} + + if need_normals or need_areas: + normals = th.cross(v1, v0, dim=-1) + norm = th.linalg.vector_norm(normals, dim=-1, keepdim=True) + + if need_areas: + output["areas"] = 0.5 * norm + + if need_normals: + output["normals"] = normals / norm.clamp(min=eps) + + if need_edges: + v2 = p2 - p1 + output["edges"] = th.stack([v0, v1, v2], dim=2) + + if len(to_compute) == 1: + return output[to_compute[0]] + else: + return output + + +def vert_binormals(v: Tensor, vt: Tensor, vi: Tensor, vti: Tensor) -> Tensor: + # Compute (dp/dt)^T + dpdt_t, vf = face_dpdt(v, vt, vi, vti) + + # Take the dp/dt.u part. Produces u vector in 3D world-space which we use for binormal vector + fbnorms = dpdt_t[:, :, 0, :] + + vbnorms = face_attribute_to_vert(v, vi, fbnorms) + return thf.normalize(vbnorms, dim=-1) + + +def vert_normals( + v: th.Tensor, vi: th.Tensor, fnorms: Optional[th.Tensor] = None +) -> th.Tensor: + """Given a set of vertices ``v`` and indices ``vi`` indexing into ``v`` + defining a set of faces, compute normals for each vertex by averaging the + face normals for each face which includes that vertex. + + Args: + v: Vertex positions, shape [batch_size, n_vertices, 3] + + vi: Vertex indices, shape [batch_size, n_faces, 3] + + fnorms: Face normals. Optional, provide them if available, otherwise they will be computed + from `v` and `vi`. Shape [n_faces, 3] + + Returns: + th.Tensor: Vertex normals, shape [batch_size, n_vertices, 3] + + """ + if fnorms is None: + fnorms = face_info(v, vi, ["normals"]) + assert isinstance(fnorms, th.Tensor) + vnorms = face_attribute_to_vert(v, vi, fnorms) + return thf.normalize(vnorms, dim=-1) diff --git a/drtk/utils/indexing.py b/drtk/utils/indexing.py new file mode 100644 index 0000000..070da78 --- /dev/null +++ b/drtk/utils/indexing.py @@ -0,0 +1,26 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +import torch as th + + +def index(x: th.Tensor, idxs: th.Tensor, dim: int) -> th.Tensor: + """Index a tensor along a given dimension using an index tensor, replacing + the shape along the given dimension with the shape of the index tensor. + + Example: + x: [8, 7306, 3] + idxs: [11000, 3] + + y = index(x, idxs, dim=1) -> y: [8, 11000, 3, 3] + with each y[b, i, j, k] = x[b, idxs[i, j], k] + + """ + + target_shape = [*x.shape] + del target_shape[dim] + target_shape[dim:dim] = [*idxs.shape] + return x.index_select(dim, idxs.view(-1)).reshape(target_shape) diff --git a/drtk/utils/projection.py b/drtk/utils/projection.py new file mode 100644 index 0000000..b76fa35 --- /dev/null +++ b/drtk/utils/projection.py @@ -0,0 +1,463 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +from typing import Optional, Sequence, Set, Tuple, Union + +import numpy as np +import torch as th + +DISTORTION_MODES: Set[Optional[str]] = {None, "radial-tangential", "fisheye"} + + +def project_pinhole( + v_cam: th.Tensor, focal: th.Tensor, princpt: th.Tensor +) -> th.Tensor: + """Project camera-space points to pixel-space points with camera + intrinsics. + + v_cam: N x V x 3 + focal: N x 2 x 2 + princpt: N x 2 + """ + + z = v_cam[:, :, 2:3] + z = th.where(z < 0, z.clamp(max=-1e-8), z.clamp(min=1e-8)) + + v_proj = v_cam[:, :, 0:2] / z + v_pix = (focal[:, None] @ v_proj[..., None])[..., 0] + princpt[:, None] + + return v_pix + + +def project_pinhole_distort_rt( + v_cam: th.Tensor, + focal: th.Tensor, + princpt: th.Tensor, + D: th.Tensor, + fov: Optional[th.Tensor] = None, +) -> th.Tensor: + """Project camera-space points to distorted pixel-space using the radial + and tangential model (4 parameters). + + v_cam: N x V x 3 + focal: N x 2 x 2 + princpt: N x 2 + D: N x 4 + fov: N x 1 + """ + + # See https://docs.opencv.org/2.4/doc/tutorials/calib3d/camera_calibration/camera_calibration.html + + if fov is None: + with th.no_grad(): + fov = estimate_rt_fov(D) + + z = v_cam[:, :, 2:3] + z = th.where(z < 0, z.clamp(max=-1e-8), z.clamp(min=1e-8)) + + v_proj = v_cam[:, :, :2] / z + r2 = v_proj.pow(2).sum(-1) + + # Clamp x, y and r to avoid wrapping behavior of the distortion model. + r2 = r2.clamp(max=fov.pow(2)) + v_clamped = v_proj.clamp(min=-fov[..., None], max=fov[..., None]) + + assert D.shape[1] in [4, 5, 8] + + # 4 param: R = (1 + k1 r^2 + k2 r^4) + R = 1 + D[:, 0:1] * r2 + D[:, 1:2] * r2.pow(2) + + # 5 param: R = (1 + k1 r^2 + k2 r^4 + k3 r^6) + if D.shape[1] == 5: + R = R + D[:, 4:5] * r2.pow(3) + + # 8 param: R = (1 + k1 r^2 + k2 r^4 + k3 r^6) / (1 + k4 r^2 + k5 r^4 + k6 r^6) + if D.shape[1] == 8: + R = R + D[:, 4:5] * r2.pow(3) + R = R / (1 + D[:, 5:6] * r2 + D[:, 6:7] * r2.pow(2) + D[:, 7:8] * r2.pow(3)) + + # [x' y'] * R + v_proj_dist = v_proj * R[..., None] + + # [2 p1 x' y', 2 p2 x' y'] + v_proj_dist += ( + 2 + * v_clamped[..., 0:1] + * v_clamped[..., 1:2] + * th.stack((D[:, 2:3], D[:, 3:4]), dim=-1) + ) + # [p2 r^2, p1 r^2] + v_proj_dist += r2[..., None] * th.stack((D[:, 3:4], D[:, 2:3]), dim=-1) + + # [2 p2 x'^2, 2 p1 y'^2] + v_proj_dist += th.stack( + ( + 2 * D[:, 3:4] * v_clamped[..., 0].pow(2), + 2 * D[:, 2:3] * v_clamped[..., 1].pow(2), + ), + dim=-1, + ) + + v_pix_dist = (focal[:, None] @ v_proj_dist[..., None])[..., 0] + princpt[:, None] + + return v_pix_dist + + +def project_fisheye_distort( + v_cam: th.Tensor, + focal: th.Tensor, + princpt: th.Tensor, + D: th.Tensor, + fov: Optional[th.Tensor] = None, +) -> th.Tensor: + """Project camera-space points to distort pixel-space points using the + fisheye distortion model. + + v_cam: N x V x 3 + focal: N x 2 x 2 + princpt: N x 2 + D: N x 4 + fov: N x 1 + """ + + # See https://github.com/opencv/opencv/blob/master/modules/calib3d/src/fisheye.cpp + + if fov is None: + with th.no_grad(): + fov = estimate_fisheye_fov(D) + + z = v_cam[:, :, 2:3] + z = th.where(z < 0, z.clamp(max=-1e-8), z.clamp(min=1e-8)) + + v_proj = v_cam[:, :, :2] / z + r = v_proj.pow(2).sum(-1).sqrt() + r = r.clamp(max=fov, min=1e-8 * th.ones_like(fov)) + theta = th.atan(r) + theta_d = theta * ( + 1 + + D[:, 0:1] * theta.pow(2) + + D[:, 1:2] * theta.pow(4) + + D[:, 2:3] * theta.pow(6) + + D[:, 3:4] * theta.pow(8) + ) + r = th.where(r < 0, r.clamp(max=-1e-8), r.clamp(min=1e-8)) + v_proj_dist = v_proj * (theta_d / r)[..., None] + v_pix_dist = (focal[:, None] @ v_proj_dist[..., None])[..., 0] + princpt[:, None] + + return v_pix_dist + + +def project_fisheye_distort_62( + v_cam: th.Tensor, + focal: th.Tensor, + princpt: th.Tensor, + D: th.Tensor, + fov: Optional[th.Tensor] = None, +) -> th.Tensor: + """Project camera-space points to distort pixel-space points using the + OculusVisionFishEye62 distortion model. + + v_cam: N x V x 3 + focal: N x 2 x 2 + princpt: N x 2 + D: N x 4 + fov: N x 1 + """ + + # See https://www.internalfb.com/code/fbsource/[188bdaeaad64]/arvr/projects/nimble/prod/pynimble/visualization/shaders.py?lines=103-123 + # a more readible version: https://euratom-software.github.io/calcam/html/intro_theory.html + if fov is None: + with th.no_grad(): + fov = estimate_fisheye_fov(D) + + z = v_cam[:, :, 2:3] + z = th.where(z < 0, z.clamp(max=-1e-8), z.clamp(min=1e-8)) + + v_proj = v_cam[:, :, :2] / z + r = v_proj.pow(2).sum(-1).sqrt() # rp + r = r.clamp(max=fov, min=1e-8 * th.ones_like(fov)) + theta = th.atan(r) + theta_d = theta * ( + 1 + + D[:, 0:1] * theta.pow(2) + + D[:, 1:2] * theta.pow(4) + + D[:, 2:3] * theta.pow(6) + + D[:, 3:4] * theta.pow(8) + + D[:, 4:5] * theta.pow(10) + + D[:, 5:6] * theta.pow(12) + ) + r = th.where(r < 0, r.clamp(max=-1e-8), r.clamp(min=1e-8)) + v_proj_dist = v_proj * (theta_d / r)[..., None] + + # Tangential Distortion + x = v_proj_dist[:, :, 0] + y = v_proj_dist[:, :, 1] + + xtan = D[:, 6:7] * (r.pow(2) + 2 * x.pow(2)) + 2 * D[:, 7:8] * x * y + ytan = 2 * D[:, 6:7] * x * y + D[:, 7:8] * (r.pow(2) + 2 * y.pow(2)) + + pTangential = th.cat([xtan[..., None], ytan[..., None]], dim=-1) + + v_proj_dist = v_proj_dist + pTangential + + v_pix_dist = (focal[:, None] @ v_proj_dist[..., None])[..., 0] + princpt[:, None] + + return v_pix_dist + + +def estimate_rt_fov(D: Union[np.ndarray, th.Tensor]) -> th.Tensor: + """Estimate the maximum field of view based on the assumption that the 5th order + polynomial for fish-eye effect is non-decreasing. + + D: N x 4 + """ + + if th.is_tensor(D): + coefs = D.cpu().numpy() + else: + coefs = D + + ones = np.ones_like(coefs[:, 0]) + zeros = np.zeros_like(coefs[:, 0]) + coefs = np.stack( + [ + 5 * coefs[:, 1], + zeros, + 3 * coefs[:, 0], + zeros, + ones, + ], + axis=-1, + ) + + fov = [] + for coef in coefs: + roots = np.roots(coef) + real_valued = roots.real[abs(roots.imag) < 1e-5] + positive_roots = real_valued[real_valued > 0] + if len(positive_roots) == 0: + fov.append(np.inf) + else: + fov.append(positive_roots.min()) + fov = np.asarray(fov, dtype=np.float32)[..., None] + + if th.is_tensor(D): + fov = th.from_numpy(fov).to(D) + + return fov + + +def estimate_fisheye_fov(D: Union[np.ndarray, th.Tensor]) -> th.Tensor: + """Estimate the maximum field of view based on the assumption that the 9th order + polynomial is non-decreasing. + + D: N x 4 + """ + + if th.is_tensor(D): + coefs = D.cpu().numpy() + else: + coefs = D + + ones = np.ones_like(coefs[:, 0]) + zeros = np.zeros_like(coefs[:, 0]) + coefs = np.stack( + [ + 9 * coefs[:, -1], + zeros, + 7 * coefs[:, -2], + zeros, + 5 * coefs[:, -3], + zeros, + 3 * coefs[:, -4], + zeros, + ones, + ], + axis=-1, + ) + + fov = [] + for coef in coefs: + roots = np.roots(coef) + real_valued = roots.real[abs(roots.imag) < 1e-5] + positive_roots = real_valued[real_valued > 0] + if len(positive_roots) == 0: + fov.append(np.pi / 2) + else: + fov.append(min(positive_roots.min(), np.pi / 2)) + fov = np.asarray(np.tan(fov), dtype=np.float32)[..., None] + + if th.is_tensor(D): + fov = th.from_numpy(fov).to(D) + + return fov + + +def project_points( + v: th.Tensor, + campos: th.Tensor, + camrot: th.Tensor, + focal: th.Tensor, + princpt: th.Tensor, + distortion_mode: Optional[Sequence[str]] = None, + distortion_coeff: Optional[th.Tensor] = None, + fov: Optional[th.Tensor] = None, +) -> Tuple[th.Tensor, th.Tensor]: + """Project 3D world-space vertices to pixel-space, optionally applying a + distortion model with provided coefficients. + + Returns v_pix, v_cam, both N x V x 3 since we preserve the camera-space + Z-coordinate for v_pix. + + v: N x V x 3 + camrot: N x 3 x 3 + campos: N x 3 + focal: N x 2 x 2 + princpt: N x 2 + distortion_coeff: N x 4 + fov: N x 1 + """ + + if distortion_mode is not None: + assert distortion_coeff is not None, "Missing distortion coefficients." + + v_cam = (camrot[:, None] @ (v - campos[:, None])[..., None])[..., 0] + + # Fall back to single distortion mode if all the distortion modes are the same. + if isinstance(distortion_mode, (list, tuple)): + modes = list(set(distortion_mode)) + if len(modes) == 0: + distortion_mode = None + elif len(modes) == 1: + distortion_mode = modes[0] + + if distortion_mode is None: + v_pix = project_pinhole(v_cam, focal, princpt) + elif isinstance(distortion_mode, str): + assert distortion_coeff is not None + + # Single distortion model + if distortion_mode == "radial-tangential": + v_pix = project_pinhole_distort_rt( + v_cam, focal, princpt, distortion_coeff, fov + ) + elif distortion_mode == "fisheye": + v_pix = project_fisheye_distort( + v_cam, focal, princpt, distortion_coeff, fov + ) + elif distortion_mode == "fisheye62": + v_pix = project_fisheye_distort_62( + v_cam, focal, princpt, distortion_coeff, fov + ) + else: + raise ValueError( + f"Invalid distortion mode: {distortion_mode}. Valid options: {DISTORTION_MODES}." + ) + elif isinstance(distortion_mode, (list, tuple)): + assert distortion_coeff is not None + + # A mix of multiple distortion modes + modes = set(distortion_mode) + if not modes <= DISTORTION_MODES: + raise ValueError( + f"Invalid distortion mode: {distortion_mode}. Valid options: {DISTORTION_MODES}." + ) + v_pix = th.empty_like(v_cam[..., :2]) + if None in modes: + idx = th.tensor( + [mode is None for mode in distortion_mode], device=v_pix.device + ) + v_pix[idx] = project_pinhole(v_cam[idx], focal[idx], princpt[idx]) + if "radial-tangential" in modes: + idx = th.tensor( + [mode == "radial-tangential" for mode in distortion_mode], + device=v_pix.device, + ) + v_pix[idx] = project_pinhole_distort_rt( + v_cam[idx], + focal[idx], + princpt[idx], + distortion_coeff[idx], + fov[idx] if fov is not None else None, + ) + if "fisheye" in modes: + idx = th.tensor( + [mode == "fisheye" for mode in distortion_mode], device=v_pix.device + ) + v_pix[idx] = project_fisheye_distort( + v_cam[idx], + focal[idx], + princpt[idx], + distortion_coeff[idx], + fov[idx] if fov is not None else None, + ) + else: + raise ValueError( + f"Invalid distortion mode: {distortion_mode}. Valid options: {DISTORTION_MODES}." + ) + + v_pix = th.cat((v_pix[:, :, 0:2], v_cam[:, :, 2:3]), dim=-1) + + return v_pix, v_cam + + +def project_points_grad( + v_grad: th.Tensor, + v: th.Tensor, + campos: th.Tensor, + camrot: th.Tensor, + focal: th.Tensor, + distortion_mode: Optional[Sequence[str]] = None, + distortion_coeff: Optional[th.Tensor] = None, +) -> th.Tensor: + """Computes the gradient of projected (pixel-space) vertex positions with + respect to the 3D world-space vertex positions given the gradient of the 3D + world-space vertex positions. + + project_points_grad(dv, v) = d project_points(v) / dv * dv + + Args: + v_grad: Gradient of 3D world-space vertices. Shape: N x V x 3 + v: 3D world-space vertices. Shape: N x V x 3 + camrot: Camera rotation. Shape: N x 3 x 3 + camrot: Camera position. Shape: N x 3 + focal: Focal length. Shape: N x 2 x 2 + distortion_mode: Distortion currently not implemented and must be None. + distortion_coeff: Distortion currently not implemented and must be None. + + Returns: + Gradient of 2D pixel-space vertices: N x V x 2 + """ + + if distortion_mode is not None: + assert distortion_coeff is not None, "Missing distortion coefficients." + + # d v_cam = d (Rv + T) = Rdv + v_cam_grad = (camrot[:, None] @ v_grad[..., None])[..., 0] + v_cam = (camrot[:, None] @ (v - campos[:, None])[..., None])[..., 0] + + if distortion_mode is None: + z = v_cam[:, :, 2:3] + z_grad = v_cam_grad[:, :, 2:3] + z = th.where(z < 0, z.clamp(max=-1e-8), z.clamp(min=1e-8)) + + # Using quotient rule: + # d (v_cam / z) = (d v_cam * z - v_cam * dz) / z^2 + v_proj_grad = (v_cam_grad[:, :, 0:2] * z - v_cam[:, :, 0:2] * z_grad) / z**2.0 + + # d v_pix = d (Kv + cp) = Kdv + v_pix_grad = (focal[:, None] @ v_proj_grad[..., None])[..., 0] + + elif distortion_mode == "radial-tangential": + raise NotImplementedError + elif distortion_mode == "fisheye": + raise NotImplementedError + else: + raise ValueError( + f"Invalid distortion mode: {distortion_mode}. Valid options: {DISTORTION_MODES}." + ) + + return v_pix_grad diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..3365ac7 --- /dev/null +++ b/setup.py @@ -0,0 +1,161 @@ +# Copyright (c) Meta Platforms, Inc. and affiliates. +# All rights reserved. +# +# This source code is licensed under the license found in the +# LICENSE file in the root directory of this source tree. + +import os +import platform +import re +import sys + +from pkg_resources import DistributionNotFound, get_distribution + +from setuptools import setup +from torch.utils.cpp_extension import BuildExtension, CUDAExtension + + +def main(debug: bool) -> None: + extra_link_args = { + "linux": ["-static-libgcc"] + ([] if debug else ["-flto"]), + "win32": ["/DEBUG"] if debug else [], + } + cxx_args = { + "linux": ["-std=c++17", "-Wall"] + + (["-O0", "-g3", "-DDEBUG"] if debug else ["-O3", "--fast-math"]), + "win32": ( + ["/std:c++17", "/MT", "/GR-", "/EHsc", '/D "NOMINMAX"'] + + ( + ["/Od", '/D "_DEBUG"'] + if debug + else ["/O2", "/fp:fast", "/GL", '/D "NDEBUG"'] + ) + ), + } + + nvcc_args = [ + "-gencode=arch=compute_72,code=sm_72", + "-gencode=arch=compute_75,code=sm_75", + "-gencode=arch=compute_80,code=sm_80", + "-gencode=arch=compute_86,code=sm_86", + "-gencode=arch=compute_90,code=sm_90", + ] + (["-O0", "-g", "-DDEBUG"] if debug else ["-O3", "--use_fast_math"]) + + # There is som issue effecting latest NVCC and pytorch 2.3.0 https://github.com/pytorch/pytorch/issues/122169 + # The workaround is adding -std=c++20 to NVCC args + nvcc_args.append("-std=c++20") + + def get_dist(name): + try: + return get_distribution(name) + except DistributionNotFound: + return None + + root_path = os.path.dirname(os.path.abspath(__file__)) + package_name = "drtk" + with open(os.path.join(root_path, "drtk", "__init__.py")) as f: + init_file = f.read() + + pattern = re.compile(r"__version__\s*=\s*\"(\d*\.\d*.\d*)\"") + groups = pattern.findall(init_file) + assert len(groups) == 1 + version = groups[0] + + if get_dist("torch") is None: + raise RuntimeError("Setup requires torch package to be installed") + + import torch as th + + assert th.cuda.is_available() + + target_os = "none" + + if sys.platform == "darwin": + target_os = "macos" + elif os.name == "posix": + target_os = "linux" + elif platform.system() == "Windows": + target_os = "win32" + + if target_os == "none": + raise RuntimeError("Could not detect platform") + if target_os == "macos": + raise RuntimeError("Platform is not supported") + + include_dir = [os.path.join(root_path, "src", "include")] + + with open("README.md") as f: + readme = f.read() + + pillow = "pillow" if get_dist("pillow-simd") is None else "pillow-simd" + + setup( + name=package_name, + version=version, + author="Reality Labs, Meta", + description="Differentiable Rendering Toolkit", + long_description=readme, + long_description_content_type="text/markdown", + license="MIT", + install_requires=["numpy", "torch", "torchvision", pillow], + ext_modules=[ + CUDAExtension( + name="drtk.rasterize_ext", + sources=[ + "src/rasterize/rasterize_module.cpp", + "src/rasterize/rasterize_kernel.cu", + ], + extra_compile_args={"cxx": cxx_args[target_os], "nvcc": nvcc_args}, + extra_link_args=extra_link_args[target_os], + include_dirs=include_dir, + ), + CUDAExtension( + "drtk.render_ext", + sources=["src/render/render_kernel.cu", "src/render/render_module.cpp"], + extra_compile_args={"cxx": cxx_args[target_os], "nvcc": nvcc_args}, + include_dirs=include_dir, + ), + CUDAExtension( + "drtk.edge_grad_ext", + sources=[ + "src/edge_grad/edge_grad_module.cpp", + "src/edge_grad/edge_grad_kernel.cu", + ], + extra_compile_args={"cxx": cxx_args[target_os], "nvcc": nvcc_args}, + include_dirs=include_dir, + ), + CUDAExtension( + "drtk.mipmap_grid_sampler_ext", + sources=[ + "src/mipmap_grid_sampler/mipmap_grid_sampler_module.cpp", + "src/mipmap_grid_sampler/mipmap_grid_sampler_kernel.cu", + ], + extra_compile_args={"cxx": cxx_args[target_os], "nvcc": nvcc_args}, + include_dirs=include_dir, + ), + CUDAExtension( + "drtk.msi_ext", + sources=[ + "src/msi/msi_module.cpp", + "src/msi/msi_kernel.cu", + ], + extra_compile_args={"cxx": cxx_args[target_os], "nvcc": nvcc_args}, + include_dirs=include_dir, + ), + CUDAExtension( + "drtk.interpolate_ext", + sources=[ + "src/interpolate/interpolate_module.cpp", + "src/interpolate/interpolate_kernel.cu", + ], + extra_compile_args={"cxx": cxx_args[target_os], "nvcc": nvcc_args}, + include_dirs=include_dir, + ), + ], + cmdclass={"build_ext": BuildExtension}, + packages=["drtk", "drtk.utils"], + ) + + +if __name__ == "__main__": + main(any(x in sys.argv for x in ["debug", "-debug", "--debug"])) diff --git a/src/edge_grad/edge_grad_kernel.cu b/src/edge_grad/edge_grad_kernel.cu new file mode 100644 index 0000000..8eb65fc --- /dev/null +++ b/src/edge_grad/edge_grad_kernel.cu @@ -0,0 +1,423 @@ +// Copyright (c) Meta Platforms, Inc. and affiliates. +// All rights reserved. +// +// This source code is licensed under the license found in the +// LICENSE file in the root directory of this source tree. + +#include +#include +#include +#include + +#include +#include "edge_grad_kernel.h" + +using namespace math; + +using at::native::fastAtomicAdd; + +template +struct TriInfo { + typedef typename math::TVec2 scalar2_t; + + const scalar2_t p_0; + const scalar2_t p_1; + const scalar2_t v_01; + const scalar2_t v_02; + const scalar2_t v_12; + const scalar_t denominator; +}; + +template +__device__ bool pix_in_tri(const TriInfo& tri, const int x, const int y) { + typedef typename math::TVec2 scalar2_t; + typedef typename math::TVec3 scalar3_t; + + if (tri.denominator != 0.f) { + const scalar2_t p = {(scalar_t)x, (scalar_t)y}; + + const scalar2_t vp0p = p - tri.p_0; + const scalar2_t vp1p = p - tri.p_1; + + scalar3_t bary = scalar3_t({ + vp1p.y * tri.v_12.x - vp1p.x * tri.v_12.y, + vp0p.x * tri.v_02.y - vp0p.y * tri.v_02.x, + vp0p.y * tri.v_01.x - vp0p.x * tri.v_01.y, + }); + bary *= sign(tri.denominator); + + const bool on_edge_or_inside = (bary.x >= 0.f) && (bary.y >= 0.f) && (bary.z >= 0.f); + + bool on_edge_0 = bary.x == 0.f; + bool on_edge_1 = bary.y == 0.f; + bool on_edge_2 = bary.z == 0.f; + + const bool is_top_left_0 = (tri.denominator > 0) + ? (tri.v_12.y < 0.f || tri.v_12.y == 0.0f && tri.v_12.x > 0.f) + : (tri.v_12.y > 0.f || tri.v_12.y == 0.0f && tri.v_12.x < 0.f); + const bool is_top_left_1 = (tri.denominator > 0) + ? (tri.v_02.y > 0.f || tri.v_02.y == 0.0f && tri.v_02.x < 0.f) + : (tri.v_02.y < 0.f || tri.v_02.y == 0.0f && tri.v_02.x > 0.f); + const bool is_top_left_2 = (tri.denominator > 0) + ? (tri.v_01.y < 0.f || tri.v_01.y == 0.0f && tri.v_01.x > 0.f) + : (tri.v_01.y > 0.f || tri.v_01.y == 0.0f && tri.v_01.x < 0.f); + + const bool is_top_left_or_inside = on_edge_or_inside && + !(on_edge_0 && !is_top_left_0 || on_edge_1 && !is_top_left_1 || + on_edge_2 && !is_top_left_2); + return is_top_left_or_inside; + } + return false; +} + +template +__device__ TriInfo +get_tri_info(const scalar_t* v_ptr, index_t v_sV, index_t v_sC, int3 vi) { + typedef typename math::TVec2 scalar2_t; + const scalar2_t p_0 = {v_ptr[v_sV * vi.x + v_sC * 0], v_ptr[v_sV * vi.x + v_sC * 1]}; + const scalar2_t p_1 = {v_ptr[v_sV * vi.y + v_sC * 0], v_ptr[v_sV * vi.y + v_sC * 1]}; + const scalar2_t p_2 = {v_ptr[v_sV * vi.z + v_sC * 0], v_ptr[v_sV * vi.z + v_sC * 1]}; + + const scalar2_t v_01 = p_1 - p_0; + const scalar2_t v_02 = p_2 - p_0; + const scalar2_t v_12 = p_2 - p_1; + + const scalar_t denominator = v_01.x * v_02.y - v_01.y * v_02.x; + + return {p_0, p_1, v_01, v_02, v_12, denominator}; +} + +template +__device__ math::TVec3 +get_tri_normal(const scalar_t* v_ptr, index_t v_sV, index_t v_sC, int3 vi) { + typedef typename math::TVec3 scalar3_t; + const scalar3_t p_0 = { + v_ptr[v_sV * vi.x + v_sC * 0], v_ptr[v_sV * vi.x + v_sC * 1], v_ptr[v_sV * vi.x + v_sC * 2]}; + const scalar3_t p_1 = { + v_ptr[v_sV * vi.y + v_sC * 0], v_ptr[v_sV * vi.y + v_sC * 1], v_ptr[v_sV * vi.y + v_sC * 2]}; + const scalar3_t p_2 = { + v_ptr[v_sV * vi.z + v_sC * 0], v_ptr[v_sV * vi.z + v_sC * 1], v_ptr[v_sV * vi.z + v_sC * 2]}; + return normalize(cross(p_0 - p_2, p_1 - p_0)); +} + +template +__device__ math::TVec2 get_db_dp( + const math::TVec2& n_varying_, + const math::TVec2& n_fixed_) { + /* + Computes derivative of the point position with respect to edge displacement + Args: + - n_varying_: Projection of the normal of the movable triangle onto the plane of + consideration (XZ or YZ) N x 3 x H x W. + - n_fixed_: Projection of the normal of the fixed triangle onto the plane of consideration + (XZ or YZ) N x 3 x H x W. + Please refer to the paper "Rasterized Edge Gradients: Handling Discontinuities Differentiably" + for details. + */ + typedef typename math::TVec2 scalar2_t; + + const auto n_varying = normalize(n_varying_); + const auto n_fixed = normalize(n_fixed_); + const scalar2_t b = {-n_fixed.y, n_fixed.x}; + const auto b_dot_varyingg = dot(b, n_varying); + return b.x / epsclamp(b_dot_varyingg) * n_varying; +} + +template +__device__ math::TVec3 load_vec3_if_valid( + const scalar_t* __restrict ptr, + index_t stride, + bool valid, + const math::TVec3& def) { + if (valid) { + return {ptr[0 * stride], ptr[1 * stride], ptr[2 * stride]}; + } + return def; +} + +template +C10_LAUNCH_BOUNDS_1(256) +__global__ void edge_grad_backward_kernel( + const index_t nthreads, + TensorInfo v_pix, + TensorInfo img, + TensorInfo index_img, + TensorInfo vi, + TensorInfo grad_output, + TensorInfo grad_v_pix_img, + const index_t memory_span) { + typedef typename math::TVec2 scalar2_t; + typedef typename math::TVec3 scalar3_t; + + const index_t v_pix_sN = v_pix.strides[0]; + const index_t v_pix_sV = v_pix.strides[1]; + const index_t v_pix_sC = v_pix.strides[2]; + + const index_t C = img.sizes[1]; + const index_t H = img.sizes[2]; + const index_t W = img.sizes[3]; + const index_t V = v_pix.sizes[1]; + + const index_t index_img_sN = index_img.strides[0]; + const index_t index_img_sH = index_img.strides[1]; + const index_t index_img_sW = index_img.strides[2]; + + const index_t img_sN = img.strides[0]; + const index_t img_sC = img.strides[1]; + const index_t img_sH = img.strides[2]; + const index_t img_sW = img.strides[3]; + + const index_t grad_output_sN = grad_output.strides[0]; + const index_t grad_output_sC = grad_output.strides[1]; + const index_t grad_output_sH = grad_output.strides[2]; + const index_t grad_output_sW = grad_output.strides[3]; + + const index_t grad_v_pix_img_sN = grad_v_pix_img.strides[0]; + const index_t grad_v_pix_img_sC = grad_v_pix_img.strides[1]; + const index_t grad_v_pix_img_sH = grad_v_pix_img.strides[2]; + const index_t grad_v_pix_img_sW = grad_v_pix_img.strides[3]; + + const index_t vi_sV = vi.strides[0]; + const index_t vi_sF = vi.strides[1]; + + CUDA_KERNEL_LOOP_TYPE(index, nthreads, index_t) { + const index_t x = index % W; + const index_t y = (index / W) % H; + const index_t n = index / (H * W); + + if (x < (W - 1) && y < (H - 1)) { + // center-right-down (CRD) + // + // *--------*--------* + // | center | right | + // | (0, 0) | (1, 0) | + // *--------*--------* + // | down | + // | (0, 1) | + // *--------* + + // Computing indicator variables + // BEGIN + // triangle indices of CRD pixels + const int32_t* __restrict index_img_ptr = index_img.data + n * index_img_sN; + const int32_t center_index = index_img_ptr[(y + 0) * index_img_sH + (x + 0) * index_img_sW]; + const int32_t right_index = index_img_ptr[(y + 0) * index_img_sH + (x + 1) * index_img_sW]; + const int32_t down_index = index_img_ptr[(y + 1) * index_img_sH + (x + 0) * index_img_sW]; + + // valid mask + const bool c_valid = (center_index >= 0); + const bool r_valid = (right_index >= 0); + const bool d_valid = (down_index >= 0); + + // vertex indices of triangles of CRD pixels + // 0,0,0 - if not valid + const int3 vi_pt_center = load_vec3_if_valid( + vi.data + center_index * vi_sV, vi_sF, c_valid, {0, 0, 0}); + const int3 vi_pt_right = load_vec3_if_valid( + vi.data + right_index * vi_sV, vi_sF, r_valid, {0, 0, 0}); + const int3 vi_pt_down = load_vec3_if_valid( + vi.data + down_index * vi_sV, vi_sF, d_valid, {0, 0, 0}); + + // center <-> right differ + const bool lr_diff = (center_index != right_index); + // center <-> down differ + const bool ud_diff = (center_index != down_index); + + // if horizontal pair (vertical edge) composed of two triangles + const bool x_both_valid = c_valid && r_valid; + // if vertical pair (horizontal edge) composed of two triangles + const bool y_both_valid = c_valid && d_valid; + + // Get CRD triangle info + const scalar_t* __restrict v_pix_ptr = v_pix.data + n * v_pix_sN; + const auto tri_center = get_tri_info(v_pix_ptr, v_pix_sV, v_pix_sC, vi_pt_center); + const auto tri_right = get_tri_info(v_pix_ptr, v_pix_sV, v_pix_sC, vi_pt_right); + const auto tri_down = get_tri_info(v_pix_ptr, v_pix_sV, v_pix_sC, vi_pt_down); + + // Compute indicators of edge type + const bool center_pix_in_right_tri = lr_diff && x_both_valid && pix_in_tri(tri_right, x, y); + const bool right_pix_in_center_tri = + lr_diff && x_both_valid && pix_in_tri(tri_center, x + 1, y); + const bool center_pix_in_down_tri = ud_diff && y_both_valid && pix_in_tri(tri_down, x, y); + const bool down_pix_in_center_tri = + ud_diff && y_both_valid && pix_in_tri(tri_center, x, y + 1); + + // Overlap flags + const bool l_over_r = center_pix_in_right_tri && (!right_pix_in_center_tri); + const bool r_over_l = right_pix_in_center_tri && (!center_pix_in_right_tri); + const bool u_over_d = center_pix_in_down_tri && (!down_pix_in_center_tri); + const bool d_over_u = down_pix_in_center_tri && (!center_pix_in_down_tri); + + // Intersection flags + const bool horiz_int = center_pix_in_right_tri && right_pix_in_center_tri; + const bool vert_int = center_pix_in_down_tri && down_pix_in_center_tri; + + // Intersection flags + const bool horiz_adjacent = + lr_diff && x_both_valid && (!center_pix_in_right_tri && !right_pix_in_center_tri); + const bool vert_adjacent = + ud_diff && y_both_valid && (!center_pix_in_down_tri && !down_pix_in_center_tri); + + // END + + // Compute image gradient dot output gradient from backward + // This is computed regardless of the edge type as long as there is an edge (lr_diff or + // ud_diff) BEGIN + const scalar_t* __restrict img_ptr = img.data + img_sN * n; + const scalar_t* __restrict grad_output_ptr = grad_output.data + grad_output_sN * n; + + scalar_t grad_dot_x = 0.f; + scalar_t grad_dot_y = 0.f; + if (lr_diff) { + const scalar_t* __restrict img_ptr_right = img_ptr + y * img_sH + (x + 1) * img_sW; + const scalar_t* __restrict img_ptr_center = img_ptr + y * img_sH + (x + 0) * img_sW; + const scalar_t* __restrict grad_output_ptr_right = + grad_output_ptr + y * grad_output_sH + (x + 1) * grad_output_sW; + const scalar_t* __restrict grad_output_ptr_center = + grad_output_ptr + y * grad_output_sH + (x + 0) * grad_output_sW; + for (size_t c = 0; c < C; ++c) { + grad_dot_x += (img_ptr_right[c * img_sC] - img_ptr_center[c * img_sC]) * + (0.5f * + (grad_output_ptr_right[c * grad_output_sC] + + grad_output_ptr_center[c * grad_output_sC])); + } + } + if (ud_diff) { + const scalar_t* __restrict img_ptr_down = img_ptr + (y + 1) * img_sH + x * img_sW; + const scalar_t* __restrict img_ptr_center = img_ptr + (y + 0) * img_sH + x * img_sW; + const scalar_t* __restrict grad_output_ptr_down = + grad_output_ptr + (y + 1) * grad_output_sH + x * grad_output_sW; + const scalar_t* __restrict grad_output_ptr_center = + grad_output_ptr + (y + 0) * grad_output_sH + x * grad_output_sW; + for (size_t c = 0; c < C; ++c) { + grad_dot_y += (img_ptr_down[c * img_sC] - img_ptr_center[c * img_sC]) * + (0.5f * + (grad_output_ptr_down[c * grad_output_sC] + + grad_output_ptr_center[c * grad_output_sC])); + } + } + // END + + scalar3_t grad_v_pix_center = {0.f, 0.f, 0.f}; + scalar3_t grad_v_pix_right = {0.f, 0.f, 0.f}; + scalar3_t grad_v_pix_down = {0.f, 0.f, 0.f}; + + const scalar3_t center_normal = get_tri_normal(v_pix_ptr, v_pix_sV, v_pix_sC, vi_pt_center); + const scalar3_t right_normal = get_tri_normal(v_pix_ptr, v_pix_sV, v_pix_sC, vi_pt_right); + const scalar3_t down_normal = get_tri_normal(v_pix_ptr, v_pix_sV, v_pix_sC, vi_pt_down); + + if (!horiz_int) { + grad_v_pix_center.x += (!c_valid || r_over_l || horiz_adjacent) ? 0.f : grad_dot_x; + grad_v_pix_right.x += (!r_valid || l_over_r || horiz_adjacent) ? 0.f : grad_dot_x; + } else { + // Center triangle moves, right fixed. + scalar2_t dbx_dp = get_db_dp( + {center_normal.x, center_normal.z}, {right_normal.x, right_normal.z}); + grad_v_pix_center.x += grad_dot_x * dbx_dp.x; + grad_v_pix_center.z += grad_dot_x * dbx_dp.y; + + // Center triangle fixed, right moves. + dbx_dp = get_db_dp( + {right_normal.x, right_normal.z}, {center_normal.x, center_normal.z}); + grad_v_pix_right.x += grad_dot_x * dbx_dp.x; + grad_v_pix_right.z += grad_dot_x * dbx_dp.y; + } + + if (!vert_int) { + grad_v_pix_center.y += (!c_valid || d_over_u || vert_adjacent) ? 0.f : grad_dot_y; + grad_v_pix_down.y += (!d_valid || u_over_d || vert_adjacent) ? 0.f : grad_dot_y; + } else { + // Center triangle moves, lower fixed. + scalar2_t dby_dp = + get_db_dp({center_normal.y, center_normal.z}, {down_normal.y, down_normal.z}); + grad_v_pix_center.y += grad_dot_y * dby_dp.x; + grad_v_pix_center.z += grad_dot_y * dby_dp.x; + + // Center triangle fixed, lower moves. + dby_dp = + get_db_dp({down_normal.y, down_normal.z}, {center_normal.y, center_normal.z}); + grad_v_pix_down.y += grad_dot_y * dby_dp.x; + grad_v_pix_down.z += grad_dot_y * dby_dp.x; + } + + // Writing grads out + // BEGIN + scalar_t* __restrict grad_v_pix_img_ptr = grad_v_pix_img.data + grad_v_pix_img_sN * n; + + // center + auto* ptr_c = grad_v_pix_img_ptr + (y + 0) * grad_v_pix_img_sH + (x + 0) * grad_v_pix_img_sW; + atomicAdd(ptr_c + 0 * grad_v_pix_img_sC, -grad_v_pix_center.x); + atomicAdd(ptr_c + 1 * grad_v_pix_img_sC, -grad_v_pix_center.y); + atomicAdd(ptr_c + 2 * grad_v_pix_img_sC, -grad_v_pix_center.z); + + // right + auto* ptr_r = grad_v_pix_img_ptr + (y + 0) * grad_v_pix_img_sH + (x + 1) * grad_v_pix_img_sW; + atomicAdd(ptr_r + 0 * grad_v_pix_img_sC, -grad_v_pix_right.x); + atomicAdd(ptr_r + 1 * grad_v_pix_img_sC, -grad_v_pix_right.y); + atomicAdd(ptr_r + 2 * grad_v_pix_img_sC, -grad_v_pix_right.z); + + // down + auto* ptr_d = grad_v_pix_img_ptr + (y + 1) * grad_v_pix_img_sH + (x + 0) * grad_v_pix_img_sW; + atomicAdd(ptr_d + 0 * grad_v_pix_img_sC, -grad_v_pix_down.x); + atomicAdd(ptr_d + 1 * grad_v_pix_img_sC, -grad_v_pix_down.y); + atomicAdd(ptr_d + 2 * grad_v_pix_img_sC, -grad_v_pix_down.z); + // END + } + } +} + +template +void edge_grad_estimator_cuda_backward_( + const int64_t count, + const torch::Tensor& v_pix, + const torch::Tensor& img, + const torch::Tensor& index_img, + const torch::Tensor& vi, + const torch::Tensor& grad_outputs, + const torch::Tensor& grad_v_pix_img) { + edge_grad_backward_kernel + <<>>( + static_cast(count), + getTensorInfo(v_pix), + getTensorInfo(img), + getTensorInfo(index_img), + getTensorInfo(vi), + getTensorInfo(grad_outputs), + getTensorInfo(grad_v_pix_img), + grad_v_pix_img.numel()); + C10_CUDA_KERNEL_LAUNCH_CHECK(); +} + +torch::Tensor edge_grad_estimator_cuda_backward( + const torch::Tensor& v_pix, + const torch::Tensor& img, + const torch::Tensor& index_img, + const torch::Tensor& vi, + const torch::Tensor& grad_outputs) { + const at::cuda::OptionalCUDAGuard device_guard(device_of(img)); + + const auto N = img.sizes()[0]; + const auto C = img.sizes()[1]; + const auto H = img.sizes()[2]; + const auto W = img.sizes()[3]; + const auto V = v_pix.sizes()[1]; + const auto count = N * H * W; + + auto grad_v_pix_img = torch::zeros({N, 3, H, W}, v_pix.options()); + + if (count > 0) { + AT_DISPATCH_FLOATING_TYPES(v_pix.scalar_type(), "edge_grad_estimator_kernel", [&] { + if (at::native::canUse32BitIndexMath(v_pix) && at::native::canUse32BitIndexMath(img) && + at::native::canUse32BitIndexMath(index_img) && at::native::canUse32BitIndexMath(vi) && + at::native::canUse32BitIndexMath(grad_outputs) && + at::native::canUse32BitIndexMath(grad_v_pix_img)) { + edge_grad_estimator_cuda_backward_( + count, v_pix, img, index_img, vi, grad_outputs, grad_v_pix_img); + } else { + edge_grad_estimator_cuda_backward_( + count, v_pix, img, index_img, vi, grad_outputs, grad_v_pix_img); + } + }); + } + return grad_v_pix_img; +} diff --git a/src/edge_grad/edge_grad_kernel.h b/src/edge_grad/edge_grad_kernel.h new file mode 100644 index 0000000..7221a39 --- /dev/null +++ b/src/edge_grad/edge_grad_kernel.h @@ -0,0 +1,14 @@ +// Copyright (c) Meta Platforms, Inc. and affiliates. +// All rights reserved. +// +// This source code is licensed under the license found in the +// LICENSE file in the root directory of this source tree. + +#pragma once + +torch::Tensor edge_grad_estimator_cuda_backward( + const torch::Tensor& v_pix, + const torch::Tensor& img, + const torch::Tensor& index_img, + const torch::Tensor& vi, + const torch::Tensor& grad_outputs); diff --git a/src/edge_grad/edge_grad_module.cpp b/src/edge_grad/edge_grad_module.cpp new file mode 100644 index 0000000..9cc27c3 --- /dev/null +++ b/src/edge_grad/edge_grad_module.cpp @@ -0,0 +1,194 @@ +// Copyright (c) Meta Platforms, Inc. and affiliates. +// All rights reserved. +// +// This source code is licensed under the license found in the +// LICENSE file in the root directory of this source tree. + +#include + +#include + +#ifndef NO_PYBIND +#include +#endif + +#include "edge_grad_kernel.h" + +// Dispatch function +torch::Tensor edge_grad_estimator( + const torch::Tensor& v_pix, + const torch::Tensor& v_pix_img, + const torch::Tensor& vi, + const torch::Tensor& img, + const torch::Tensor& index_img) { + static auto op = torch::Dispatcher::singleton() + .findSchemaOrThrow("edge_grad_ext::edge_grad_estimator", "") + .typed(); + return op.call(v_pix, v_pix_img, vi, img, index_img); +} + +torch::Tensor edge_grad_estimator_fwd( + const torch::Tensor& v_pix, + const torch::Tensor& v_pix_img, + const torch::Tensor& vi, + const torch::Tensor& img, + const torch::Tensor& index_img) { + TORCH_CHECK( + v_pix.defined() && v_pix_img.defined() && vi.defined() && img.defined() && + index_img.defined(), + "edge_grad_estimator(): expected all inputs to be defined"); + TORCH_CHECK( + (v_pix.device() == v_pix_img.device()) && (v_pix.device() == vi.device()) && + (v_pix.device() == img.device()) && (v_pix.device() == index_img.device()) && + (v_pix.is_cuda()), + "edge_grad_estimator(): expected all inputs to be on same cuda device"); + TORCH_CHECK( + v_pix.is_floating_point() && v_pix_img.is_floating_point() && img.is_floating_point(), + "edge_grad_estimator(): expected v_pix, v_pix_img, and img to have floating point type, but v_pix has ", + v_pix.dtype(), + " v_pix has ", + v_pix_img.dtype(), + " img has ", + img.dtype()); + TORCH_CHECK( + vi.dtype() == torch::kInt32, + "edge_grad_estimator(): expected vi to have int32 type, but vi has ", + vi.dtype()); + TORCH_CHECK( + index_img.dtype() == torch::kInt32, + "edge_grad_estimator(): expected index_img to have int32 type, but index_img has ", + index_img.dtype()); + TORCH_CHECK( + v_pix.layout() == torch::kStrided && v_pix_img.layout() == torch::kStrided && + vi.layout() == torch::kStrided && img.layout() == torch::kStrided && + index_img.layout() == torch::kStrided, + "edge_grad_estimator(): expected all inputs to have torch.strided layout"); + TORCH_CHECK( + (v_pix.dim() == 3) && (v_pix_img.dim() == 4) && (vi.dim() == 2) && (img.dim() == 4) && + (index_img.dim() == 3), + "edge_grad_estimator(): expected v_pix.ndim == 3, v_pix_img.ndim == 4, vi.ndim == 2, img.ndim == 4, index_img.ndim == 3, " + "but got v_pix with sizes ", + v_pix.sizes(), + " and v_pix_img with sizes ", + v_pix_img.sizes(), + " and vi with sizes ", + vi.sizes(), + " and img with sizes ", + img.sizes(), + " and index_img with sizes ", + index_img.sizes()); + TORCH_CHECK( + v_pix.size(0) == v_pix_img.size(0) && v_pix.size(0) == img.size(0) && + v_pix.size(0) == index_img.size(0), + "edge_grad_estimator(): expected v and index_img to have same batch size, " + "but got v_pix with sizes ", + v_pix.sizes(), + ", v_pix_img with sizes ", + v_pix_img.sizes(), + ", img with sizes ", + img.sizes(), + " and index_img with sizes ", + index_img.sizes()); + TORCH_CHECK( + v_pix.size(2) == 3 && v_pix_img.size(1) == 3 && vi.size(1) == 3, + "edge_grad_estimator(): expected third dim of v_pix to be of size 3, and second dim of vi to be of size 3, but got ", + v_pix.size(2), + " in the third dim of v_pix, and ", + v_pix_img.size(1), + " in the second dim of v_pix_img, and ", + vi.size(1), + " in the second dim of vi"); + TORCH_CHECK( + v_pix_img.size(3) == img.size(3) && v_pix_img.size(3) == index_img.size(2) && + v_pix_img.size(2) == img.size(2) && v_pix_img.size(2) == index_img.size(1), + "edge_grad_estimator(): expected width and height of v_pix_img, img, and index_img to match, but got size of v_pix_img: ", + v_pix_img.sizes(), + ", size of img: ", + img.sizes(), + ", size of index_img: ", + index_img.sizes()); + return img; +} + +// Ideally we would need to turn off autograd handling and re-dispatch, but we just call +// cuda kernels directly +class EdgeGradEstimatorFunction : public torch::autograd::Function { + public: + static torch::autograd::tensor_list forward( + torch::autograd::AutogradContext* ctx, + const torch::Tensor& v_pix, + const torch::Tensor& v_pix_img, + const torch::Tensor& vi, + const torch::Tensor& img, + const torch::Tensor& index_img) { + ctx->set_materialize_grads(false); + ctx->save_for_backward({v_pix, img, index_img, vi}); + ctx->saved_data["v_pix_img_requires_grad"] = v_pix_img.requires_grad(); + return {img}; + } + + static torch::autograd::tensor_list backward( + torch::autograd::AutogradContext* ctx, + torch::autograd::tensor_list grad_outputs) { + // If v_pix_img doesn't require grad, we don't need to do anything. + if (!ctx->saved_data["v_pix_img_requires_grad"].toBool()) { + return {torch::Tensor(), torch::Tensor(), torch::Tensor(), grad_outputs[0], torch::Tensor()}; + } + const auto saved = ctx->get_saved_variables(); + const auto& v_pix = saved[0]; + const auto& img = saved[1]; + const auto& index_img = saved[2]; + const auto& vi = saved[3]; + + auto grad_v_pix_img = + edge_grad_estimator_cuda_backward(v_pix, img, index_img, vi, grad_outputs[0]); + return {torch::Tensor(), grad_v_pix_img, torch::Tensor(), grad_outputs[0], torch::Tensor()}; + } +}; + +torch::Tensor edge_grad_estimator_autograd( + const torch::Tensor& v_pix, + const torch::Tensor& v_pix_img, + const torch::Tensor& vi, + const torch::Tensor& img, + const torch::Tensor& index_img) { + return EdgeGradEstimatorFunction::apply(v_pix, v_pix_img, vi, img, index_img)[0]; +} + +torch::Tensor edge_grad_estimator_autocast( + const torch::Tensor& v_pix, + const torch::Tensor& v_pix_img, + const torch::Tensor& vi, + const torch::Tensor& img, + const torch::Tensor& index_img) { + c10::impl::ExcludeDispatchKeyGuard no_autocast(c10::DispatchKey::Autocast); + return edge_grad_estimator( + at::autocast::cached_cast(torch::kFloat32, v_pix), + at::autocast::cached_cast(torch::kFloat32, v_pix_img), + vi, + at::autocast::cached_cast(torch::kFloat32, img), + index_img)[0]; +} + +#ifndef NO_PYBIND +// Just so that we can import this file as a Python module to get the path and +// import the Torch ops. +PYBIND11_MODULE(edge_grad_ext, m) {} +#endif + +TORCH_LIBRARY(edge_grad_ext, m) { + m.def( + "edge_grad_estimator(Tensor v_pix, Tensor v_pix_img, Tensor vi, Tensor img, Tensor index_img) -> Tensor"); +} + +TORCH_LIBRARY_IMPL(edge_grad_ext, Autograd, m) { + m.impl("edge_grad_estimator", &edge_grad_estimator_autograd); +} + +TORCH_LIBRARY_IMPL(edge_grad_ext, Autocast, m) { + m.impl("edge_grad_estimator", edge_grad_estimator_autocast); +} + +TORCH_LIBRARY_IMPL(edge_grad_ext, CUDA, m) { + m.impl("edge_grad_estimator", &edge_grad_estimator_fwd); +} diff --git a/src/include/cuda_math_helper.h b/src/include/cuda_math_helper.h new file mode 100644 index 0000000..c8cc781 --- /dev/null +++ b/src/include/cuda_math_helper.h @@ -0,0 +1,1133 @@ +// Copyright (c) Meta Platforms, Inc. and affiliates. + +// This source code is licensed under the MIT license found in the +// LICENSE file in the root directory of this source tree. + +#pragma once + +#include +#include +#include +#include + +// The header provides uniform GLSL-like math API for the following three cases: +// - non-NVCC compiler +// - NVCC compiler, host code +// - NVCC compiler, device code +// Designed to be a more flexible replacement of similar header from NVidia +#ifndef __CUDACC__ +#define MH_NON_NVCC +#else +#define MH_NVCC +#ifdef __CUDA_ARCH__ +#define MH_NVCC_DEVICE +#else +#define MH_NVCC_HOST +#endif +#endif + +#if defined(MH_NVCC_HOST) || defined(MH_NON_NVCC) +#define HOST_DEVICE_DISPATCH(HOST_CODE, DEVICE_CODE) (HOST_CODE) +#elif defined(MH_NVCC_DEVICE) +#define HOST_DEVICE_DISPATCH(HOST_CODE, DEVICE_CODE) (DEVICE_CODE) +#else +#error Dispatch failed +#endif + +// if not NVCC, need to include cmath, since certain builtin NVCC functions have +// equivalent ones in cmath +#ifdef MH_NON_NVCC +#include +#endif + +#define CHD_FUNC constexpr inline __host__ __device__ +#define HD_FUNC inline __host__ __device__ + +namespace math { + +template +struct epsilon; + +template <> +struct epsilon { + static constexpr float value = 1e-8f; +}; + +template <> +struct epsilon { + static constexpr double value = 1e-16; +}; + +// Host and device version of saturate +// Note that unfortunately `__saturatef` aka `saturate` is a device only +// function. If you do `using namespace math` you would still have to use math +// namespace for scalars: math::saturate +HD_FUNC float saturate(float a) { + return HOST_DEVICE_DISPATCH( + fminf(fmaxf(a, 0.0f), 1.0f), + __saturatef(a) // __saturatef is a device only function + ); +} +// There is no CUDA intrinsic for saturate for double type +HD_FUNC double saturate(double a) { + return fmin(fmax(a, 0.0), 1.0); +} + +// If NVCC then use builtin abs/max/min/sqrt/rsqrt. +// All of them have overloads for ints, floats, and doubles,defined in +// `cuda/crt/math_functions.hpp` thus no need for explicit usage of e.g. fabsf +#if defined(MH_NVCC) +using ::abs; +using ::max; +using ::min; +using ::rsqrt; +using ::sqrt; +#else +// Otherwise use the ones from cmath +using std::abs; +using std::max; +using std::min; +using std::sqrt; + +inline double rsqrt(double v) { + return 1.0 / std::sqrt(v); +} +inline float rsqrt(float v) { + return 1.0f / std::sqrt(v); +} +#endif + +namespace detail { +// Provide overloads of norm3d/norm4d for floats and doubles +HD_FUNC float norm3d(float a, float b, float c) { + return HOST_DEVICE_DISPATCH( + sqrt(a * a + b * b + c * c), ::norm3df(a, b, c) // norm3df is device only + ); +} +HD_FUNC double norm3d(double a, double b, double c) { + return HOST_DEVICE_DISPATCH( + sqrt(a * a + b * b + c * c), ::norm3d(a, b, c) // norm3d is device only + ); +} +HD_FUNC float rnorm3d(float a, float b, float c) { + return HOST_DEVICE_DISPATCH( + 1.0f / sqrt(a * a + b * b + c * c), ::rnorm3df(a, b, c) // rnorm3df is device only + ); +} +HD_FUNC double rnorm3d(double a, double b, double c) { + return HOST_DEVICE_DISPATCH( + 1.0 / sqrt(a * a + b * b + c * c), ::rnorm3d(a, b, c) // rnorm3d is device only + ); +} +HD_FUNC float norm4d(float a, float b, float c, float d) { + return HOST_DEVICE_DISPATCH( + sqrt(a * a + b * b + c * c + d * d), ::norm4df(a, b, c, d) // norm4df is device only + ); +} +HD_FUNC double norm4d(double a, double b, double c, double d) { + return HOST_DEVICE_DISPATCH( + sqrt(a * a + b * b + c * c + d * d), ::norm4d(a, b, c, d) // norm4d is device only + ); +} +HD_FUNC float rnorm4d(float a, float b, float c, float d) { + return HOST_DEVICE_DISPATCH( + 1.0f / sqrt(a * a + b * b + c * c + d * d), ::rnorm4df(a, b, c, d) // rnorm4df is device only + ); +} +HD_FUNC double rnorm4d(double a, double b, double c, double d) { + return HOST_DEVICE_DISPATCH( + 1.0 / sqrt(a * a + b * b + c * c + d * d), ::rnorm4d(a, b, c, d) // rnorm4d is device only + ); +} +} // namespace detail + +// Unary operators +#define UNARY_OP(T, T2, T3, T4) \ + CHD_FUNC T2 operator+(T2 const& v) { \ + return v; \ + } \ + CHD_FUNC T2 operator-(T2 const& v) { \ + return {-v.x, -v.y}; \ + } \ + CHD_FUNC T3 operator+(T3 const& v) { \ + return v; \ + } \ + CHD_FUNC T3 operator-(T3 const& v) { \ + return {-v.x, -v.y, -v.z}; \ + } \ + CHD_FUNC T4 operator+(T4 const& v) { \ + return v; \ + } \ + CHD_FUNC T4 operator-(T4 const& v) { \ + return {-v.x, -v.y, -v.z, -v.w}; \ + } + +// -- Binary arithmetic operators -- +#define BINARY_ARITHM_OP(T, T2, T3, T4) \ + CHD_FUNC T2 operator+(T2 const& v, T scalar) { \ + return {v.x + scalar, v.y + scalar}; \ + } \ + CHD_FUNC T2 operator+(T scalar, T2 const& v) { \ + return {scalar + v.x, scalar + v.y}; \ + } \ + CHD_FUNC T2 operator+(T2 const& v1, T2 const& v2) { \ + return {v1.x + v2.x, v1.y + v2.y}; \ + } \ + CHD_FUNC T2 operator+=(T2& v, T scalar) { \ + v.x += scalar; \ + v.y += scalar; \ + return v; \ + } \ + CHD_FUNC T2 operator+=(T2& v, T2 const& v2) { \ + v.x += v2.x; \ + v.y += v2.y; \ + return v; \ + } \ + CHD_FUNC T2 operator-(T2 const& v, T scalar) { \ + return {v.x - scalar, v.y - scalar}; \ + } \ + CHD_FUNC T2 operator-(T scalar, T2 const& v) { \ + return {scalar - v.x, scalar - v.y}; \ + } \ + CHD_FUNC T2 operator-(T2 const& v1, T2 const& v2) { \ + return {v1.x - v2.x, v1.y - v2.y}; \ + } \ + CHD_FUNC T2 operator-=(T2& v, T scalar) { \ + v.x -= scalar; \ + v.y -= scalar; \ + return v; \ + } \ + CHD_FUNC T2 operator-=(T2& v, T2 const& v2) { \ + v.x -= v2.x; \ + v.y -= v2.y; \ + return v; \ + } \ + CHD_FUNC T2 operator*(T2 const& v, T scalar) { \ + return {v.x * scalar, v.y * scalar}; \ + } \ + CHD_FUNC T2 operator*(T scalar, T2 const& v) { \ + return {scalar * v.x, scalar * v.y}; \ + } \ + CHD_FUNC T2 operator*(T2 const& v1, T2 const& v2) { \ + return {v1.x * v2.x, v1.y * v2.y}; \ + } \ + CHD_FUNC T2 operator*=(T2& v, T scalar) { \ + v.x *= scalar; \ + v.y *= scalar; \ + return v; \ + } \ + CHD_FUNC T2 operator*=(T2& v, T2 const& v2) { \ + v.x *= v2.x; \ + v.y *= v2.y; \ + return v; \ + } \ + CHD_FUNC T2 operator/(T2 const& v, T scalar) { \ + return {v.x / scalar, v.y / scalar}; \ + } \ + CHD_FUNC T2 operator/(T scalar, T2 const& v) { \ + return {scalar / v.x, scalar / v.y}; \ + } \ + CHD_FUNC T2 operator/(T2 const& v1, T2 const& v2) { \ + return {v1.x / v2.x, v1.y / v2.y}; \ + } \ + CHD_FUNC T2 operator/=(T2& v, T scalar) { \ + v.x /= scalar; \ + v.y /= scalar; \ + return v; \ + } \ + CHD_FUNC T2 operator/=(T2& v, T2 const& v2) { \ + v.x /= v2.x; \ + v.y /= v2.y; \ + return v; \ + } \ + CHD_FUNC T3 operator+(T3 const& v, T scalar) { \ + return {v.x + scalar, v.y + scalar, v.z + scalar}; \ + } \ + CHD_FUNC T3 operator+(T scalar, T3 const& v) { \ + return {scalar + v.x, scalar + v.y, scalar + v.z}; \ + } \ + CHD_FUNC T3 operator+(T3 const& v1, T3 const& v2) { \ + return {v1.x + v2.x, v1.y + v2.y, v1.z + v2.z}; \ + } \ + CHD_FUNC T3 operator+=(T3& v, T scalar) { \ + v.x += scalar; \ + v.y += scalar; \ + v.z += scalar; \ + return v; \ + } \ + CHD_FUNC T3 operator+=(T3& v, T3 const& v2) { \ + v.x += v2.x; \ + v.y += v2.y; \ + v.z += v2.z; \ + return v; \ + } \ + CHD_FUNC T3 operator-(T3 const& v, T scalar) { \ + return {v.x - scalar, v.y - scalar, v.z - scalar}; \ + } \ + CHD_FUNC T3 operator-(T scalar, T3 const& v) { \ + return {scalar - v.x, scalar - v.y, scalar - v.z}; \ + } \ + CHD_FUNC T3 operator-(T3 const& v1, T3 const& v2) { \ + return {v1.x - v2.x, v1.y - v2.y, v1.z - v2.z}; \ + } \ + CHD_FUNC T3 operator-=(T3& v, T scalar) { \ + v.x -= scalar; \ + v.y -= scalar; \ + v.z -= scalar; \ + return v; \ + } \ + CHD_FUNC T3 operator-=(T3& v, T3 const& v2) { \ + v.x -= v2.x; \ + v.y -= v2.y; \ + v.z -= v2.z; \ + return v; \ + } \ + CHD_FUNC T3 operator*(T3 const& v, T scalar) { \ + return {v.x * scalar, v.y * scalar, v.z * scalar}; \ + } \ + CHD_FUNC T3 operator*(T scalar, T3 const& v) { \ + return {scalar * v.x, scalar * v.y, scalar * v.z}; \ + } \ + CHD_FUNC T3 operator*(T3 const& v1, T3 const& v2) { \ + return {v1.x * v2.x, v1.y * v2.y, v1.z * v2.z}; \ + } \ + CHD_FUNC T3 operator*=(T3& v, T scalar) { \ + v.x *= scalar; \ + v.y *= scalar; \ + v.z *= scalar; \ + return v; \ + } \ + CHD_FUNC T3 operator*=(T3& v, T3 const& v2) { \ + v.x *= v2.x; \ + v.y *= v2.y; \ + v.z *= v2.z; \ + return v; \ + } \ + CHD_FUNC T3 operator/(T3 const& v, T scalar) { \ + return {v.x / scalar, v.y / scalar, v.z / scalar}; \ + } \ + CHD_FUNC T3 operator/(T scalar, T3 const& v) { \ + return {scalar / v.x, scalar / v.y, scalar / v.z}; \ + } \ + CHD_FUNC T3 operator/(T3 const& v1, T3 const& v2) { \ + return {v1.x / v2.x, v1.y / v2.y, v1.z / v2.z}; \ + } \ + CHD_FUNC T3 operator/=(T3& v, T scalar) { \ + v.x /= scalar; \ + v.y /= scalar; \ + v.z /= scalar; \ + return v; \ + } \ + CHD_FUNC T3 operator/=(T3& v, T3 const& v2) { \ + v.x /= v2.x; \ + v.y /= v2.y; \ + v.z /= v2.z; \ + return v; \ + } \ + CHD_FUNC T4 operator+(T4 const& v, T scalar) { \ + return {v.x + scalar, v.y + scalar, v.z + scalar, v.w + scalar}; \ + } \ + CHD_FUNC T4 operator+(T scalar, T4 const& v) { \ + return {scalar + v.x, scalar + v.y, scalar + v.z, scalar + v.w}; \ + } \ + CHD_FUNC T4 operator+(T4 const& v1, T4 const& v2) { \ + return {v1.x + v2.x, v1.y + v2.y, v1.z + v2.z, v1.w + v2.w}; \ + } \ + CHD_FUNC T4 operator+=(T4& v, T scalar) { \ + v.x += scalar; \ + v.y += scalar; \ + v.z += scalar; \ + v.w += scalar; \ + return v; \ + } \ + CHD_FUNC T4 operator+=(T4& v, T4 const& v2) { \ + v.x += v2.x; \ + v.y += v2.y; \ + v.z += v2.z; \ + v.w += v2.w; \ + return v; \ + } \ + CHD_FUNC T4 operator-(T4 const& v, T scalar) { \ + return {v.x - scalar, v.y - scalar, v.z - scalar, v.w - scalar}; \ + } \ + CHD_FUNC T4 operator-(T scalar, T4 const& v) { \ + return {scalar - v.x, scalar - v.y, scalar - v.z, scalar - v.w}; \ + } \ + CHD_FUNC T4 operator-(T4 const& v1, T4 const& v2) { \ + return {v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.w - v2.w}; \ + } \ + CHD_FUNC T4 operator-=(T4& v, T scalar) { \ + v.x -= scalar; \ + v.y -= scalar; \ + v.z -= scalar; \ + v.w -= scalar; \ + return v; \ + } \ + CHD_FUNC T4 operator-=(T4& v, T4 const& v2) { \ + v.x -= v2.x; \ + v.y -= v2.y; \ + v.z -= v2.z; \ + v.w -= v2.w; \ + return v; \ + } \ + CHD_FUNC T4 operator*(T4 const& v, T scalar) { \ + return {v.x * scalar, v.y * scalar, v.z * scalar, v.w * scalar}; \ + } \ + CHD_FUNC T4 operator*(T scalar, T4 const& v) { \ + return {scalar * v.x, scalar * v.y, scalar * v.z, scalar * v.w}; \ + } \ + CHD_FUNC T4 operator*(T4 const& v1, T4 const& v2) { \ + return {v1.x * v2.x, v1.y * v2.y, v1.z * v2.z, v1.w * v2.w}; \ + } \ + CHD_FUNC T4 operator*=(T4& v, T scalar) { \ + v.x *= scalar; \ + v.y *= scalar; \ + v.z *= scalar; \ + v.w *= scalar; \ + return v; \ + } \ + CHD_FUNC T4 operator*=(T4& v, T4 const& v2) { \ + v.x *= v2.x; \ + v.y *= v2.y; \ + v.z *= v2.z; \ + v.w *= v2.w; \ + return v; \ + } \ + CHD_FUNC T4 operator/(T4 const& v, T scalar) { \ + return {v.x / scalar, v.y / scalar, v.z / scalar, v.w / scalar}; \ + } \ + CHD_FUNC T4 operator/(T scalar, T4 const& v) { \ + return {scalar / v.x, scalar / v.y, scalar / v.z, scalar / v.w}; \ + } \ + CHD_FUNC T4 operator/(T4 const& v1, T4 const& v2) { \ + return {v1.x / v2.x, v1.y / v2.y, v1.z / v2.z, v1.w / v2.w}; \ + } \ + CHD_FUNC T4 operator/=(T4& v, T scalar) { \ + v.x /= scalar; \ + v.y /= scalar; \ + v.z /= scalar; \ + v.w /= scalar; \ + return v; \ + } \ + CHD_FUNC T4 operator/=(T4& v, T4 const& v2) { \ + v.x /= v2.x; \ + v.y /= v2.y; \ + v.z /= v2.z; \ + v.w /= v2.w; \ + return v; \ + } + +#define BINARY_INT_OP(T, T2, T3, T4) \ + CHD_FUNC T2 operator%(T2 const& v, T scalar) { \ + return {(T)(v.x % scalar), (T)(v.y % scalar)}; \ + } \ + CHD_FUNC T2 operator%(T scalar, T2 const& v) { \ + return {(T)(scalar % v.x), (T)(scalar % v.y)}; \ + } \ + CHD_FUNC T2 operator%(T2 const& v1, T2 const& v2) { \ + return {(T)(v1.x % v2.x), (T)(v1.y % v2.y)}; \ + } \ + CHD_FUNC T3 operator%(T3 const& v, T scalar) { \ + return {(T)(v.x % scalar), (T)(v.y % scalar), (T)(v.z % scalar)}; \ + } \ + CHD_FUNC T3 operator%(T scalar, T3 const& v) { \ + return {(T)(scalar % v.x), (T)(scalar % v.y), (T)(scalar % v.z)}; \ + } \ + CHD_FUNC T3 operator%(T3 const& v1, T3 const& v2) { \ + return {(T)(v1.x % v2.x), (T)(v1.y % v2.y), (T)(v1.z % v2.z)}; \ + } \ + CHD_FUNC T4 operator%(T4 const& v, T scalar) { \ + return {(T)(v.x % scalar), (T)(v.y % scalar), (T)(v.z % scalar), (T)(v.w % scalar)}; \ + } \ + CHD_FUNC T4 operator%(T scalar, T4 const& v) { \ + return {(T)(scalar % v.x), (T)(scalar % v.y), (T)(scalar % v.z), (T)(scalar % v.w)}; \ + } \ + CHD_FUNC T4 operator%(T4 const& v1, T4 const& v2) { \ + return {(T)(v1.x % v2.x), (T)(v1.y % v2.y), (T)(v1.z % v2.z), (T)(v1.w % v2.w)}; \ + } + +// -- Binary bit operators -- +#define BINARY_BIT_OP(T, T2, T3, T4) \ + CHD_FUNC T2 operator&(T2 const& v, T scalar) { \ + return {v.x & scalar, v.y & scalar}; \ + } \ + CHD_FUNC T2 operator&(T scalar, T2 const& v) { \ + return {scalar & v.x, scalar & v.y}; \ + } \ + CHD_FUNC T2 operator&(T2 const& v1, T2 const& v2) { \ + return {v1.x & v2.x, v1.y & v2.y}; \ + } \ + CHD_FUNC T2 operator|(T2 const& v, T scalar) { \ + return {v.x | scalar, v.y | scalar}; \ + } \ + CHD_FUNC T2 operator|(T scalar, T2 const& v) { \ + return {scalar | v.x, scalar | v.y}; \ + } \ + CHD_FUNC T2 operator|(T2 const& v1, T2 const& v2) { \ + return {v1.x | v2.x, v1.y | v2.y}; \ + } \ + CHD_FUNC T2 operator^(T2 const& v, T scalar) { \ + return {v.x ^ scalar, v.y ^ scalar}; \ + } \ + CHD_FUNC T2 operator^(T scalar, T2 const& v) { \ + return {scalar ^ v.x, scalar ^ v.y}; \ + } \ + CHD_FUNC T2 operator^(T2 const& v1, T2 const& v2) { \ + return {v1.x ^ v2.x, v1.y ^ v2.y}; \ + } \ + CHD_FUNC T2 operator<<(T2 const& v, T scalar) { \ + return {v.x << scalar, v.y << scalar}; \ + } \ + CHD_FUNC T2 operator<<(T scalar, T2 const& v) { \ + return {scalar << v.x, scalar << v.y}; \ + } \ + CHD_FUNC T2 operator<<(T2 const& v1, T2 const& v2) { \ + return {v1.x << v2.x, v1.y << v2.y}; \ + } \ + CHD_FUNC T2 operator>>(T2 const& v, T scalar) { \ + return {v.x >> scalar, v.y >> scalar}; \ + } \ + CHD_FUNC T2 operator>>(T scalar, T2 const& v) { \ + return {scalar >> v.x, scalar >> v.y}; \ + } \ + CHD_FUNC T2 operator>>(T2 const& v1, T2 const& v2) { \ + return {v1.x >> v2.x, v1.y >> v2.y}; \ + } \ + CHD_FUNC T2 operator~(T2 const& v) { \ + return {~v.x, ~v.y}; \ + } \ + CHD_FUNC T3 operator&(T3 const& v, T scalar) { \ + return {v.x & scalar, v.y & scalar, v.z & scalar}; \ + } \ + CHD_FUNC T3 operator&(T scalar, T3 const& v) { \ + return {scalar & v.x, scalar & v.y, scalar & v.z}; \ + } \ + CHD_FUNC T3 operator&(T3 const& v1, T3 const& v2) { \ + return {v1.x & v2.x, v1.y & v2.y, v1.z & v2.z}; \ + } \ + CHD_FUNC T3 operator|(T3 const& v, T scalar) { \ + return {v.x | scalar, v.y | scalar, v.z | scalar}; \ + } \ + CHD_FUNC T3 operator|(T scalar, T3 const& v) { \ + return {scalar | v.x, scalar | v.y, scalar | v.z}; \ + } \ + CHD_FUNC T3 operator|(T3 const& v1, T3 const& v2) { \ + return {v1.x | v2.x, v1.y | v2.y, v1.z | v2.z}; \ + } \ + CHD_FUNC T3 operator^(T3 const& v, T scalar) { \ + return {v.x ^ scalar, v.y ^ scalar, v.z ^ scalar}; \ + } \ + CHD_FUNC T3 operator^(T scalar, T3 const& v) { \ + return {scalar ^ v.x, scalar ^ v.y, scalar ^ v.z}; \ + } \ + CHD_FUNC T3 operator^(T3 const& v1, T3 const& v2) { \ + return {v1.x ^ v2.x, v1.y ^ v2.y, v1.z ^ v2.z}; \ + } \ + CHD_FUNC T3 operator<<(T3 const& v, T scalar) { \ + return {v.x << scalar, v.y << scalar, v.z << scalar}; \ + } \ + CHD_FUNC T3 operator<<(T scalar, T3 const& v) { \ + return {scalar << v.x, scalar << v.y, scalar << v.z}; \ + } \ + CHD_FUNC T3 operator<<(T3 const& v1, T3 const& v2) { \ + return {v1.x << v2.x, v1.y << v2.y, v1.z << v2.z}; \ + } \ + CHD_FUNC T3 operator>>(T3 const& v, T scalar) { \ + return {v.x >> scalar, v.y >> scalar, v.z >> scalar}; \ + } \ + CHD_FUNC T3 operator>>(T scalar, T3 const& v) { \ + return {scalar >> v.x, scalar >> v.y, scalar >> v.z}; \ + } \ + CHD_FUNC T3 operator>>(T3 const& v1, T3 const& v2) { \ + return {v1.x >> v2.x, v1.y >> v2.y, v1.z >> v2.z}; \ + } \ + CHD_FUNC T3 operator~(T3 const& v) { \ + return {~v.x, ~v.y, ~v.z}; \ + } \ + CHD_FUNC T4 operator&(T4 const& v, T scalar) { \ + return {v.x & scalar, v.y & scalar, v.z & scalar, v.w & scalar}; \ + } \ + CHD_FUNC T4 operator&(T scalar, T4 const& v) { \ + return {scalar & v.x, scalar & v.y, scalar & v.z, scalar & v.w}; \ + } \ + CHD_FUNC T4 operator&(T4 const& v1, T4 const& v2) { \ + return {v1.x & v2.x, v1.y & v2.y, v1.z & v2.z, v1.w & v2.w}; \ + } \ + CHD_FUNC T4 operator|(T4 const& v, T scalar) { \ + return {v.x | scalar, v.y | scalar, v.z | scalar, v.w | scalar}; \ + } \ + CHD_FUNC T4 operator|(T scalar, T4 const& v) { \ + return {scalar | v.x, scalar | v.y, scalar | v.z, scalar | v.w}; \ + } \ + CHD_FUNC T4 operator|(T4 const& v1, T4 const& v2) { \ + return {v1.x | v2.x, v1.y | v2.y, v1.z | v2.z, v1.w | v2.w}; \ + } \ + CHD_FUNC T4 operator^(T4 const& v, T scalar) { \ + return {v.x ^ scalar, v.y ^ scalar, v.z ^ scalar, v.w ^ scalar}; \ + } \ + CHD_FUNC T4 operator^(T scalar, T4 const& v) { \ + return {scalar ^ v.x, scalar ^ v.y, scalar ^ v.z, scalar ^ v.w}; \ + } \ + CHD_FUNC T4 operator^(T4 const& v1, T4 const& v2) { \ + return {v1.x ^ v2.x, v1.y ^ v2.y, v1.z ^ v2.z, v1.w ^ v2.w}; \ + } \ + CHD_FUNC T4 operator<<(T4 const& v, T scalar) { \ + return {v.x << scalar, v.y << scalar, v.z << scalar, v.w << scalar}; \ + } \ + CHD_FUNC T4 operator<<(T scalar, T4 const& v) { \ + return {scalar << v.x, scalar << v.y, scalar << v.z, scalar << v.w}; \ + } \ + CHD_FUNC T4 operator<<(T4 const& v1, T4 const& v2) { \ + return {v1.x << v2.x, v1.y << v2.y, v1.z << v2.z, v1.w << v2.w}; \ + } \ + CHD_FUNC T4 operator>>(T4 const& v, T scalar) { \ + return {v.x >> scalar, v.y >> scalar, v.z >> scalar, v.w >> scalar}; \ + } \ + CHD_FUNC T4 operator>>(T scalar, T4 const& v) { \ + return {scalar >> v.x, scalar >> v.y, scalar >> v.z, scalar >> v.w}; \ + } \ + CHD_FUNC T4 operator>>(T4 const& v1, T4 const& v2) { \ + return {v1.x >> v2.x, v1.y >> v2.y, v1.z >> v2.z, v1.w >> v2.w}; \ + } \ + CHD_FUNC T4 operator~(T4 const& v) { \ + return {~v.x, ~v.y, ~v.z, ~v.w}; \ + } + +#define BINARY_EQ_OP(T, T2, T3, T4) \ + CHD_FUNC bool operator==(T2 const& v1, T2 const& v2) { \ + return v1.x == v2.x && v1.y == v2.y; \ + } \ + CHD_FUNC bool operator!=(T2 const& v1, T2 const& v2) { \ + return !(v1 == v2); \ + } \ + CHD_FUNC bool operator==(T3 const& v1, T3 const& v2) { \ + return v1.x == v2.x && v1.y == v2.y && v1.z == v2.z; \ + } \ + CHD_FUNC bool operator!=(T3 const& v1, T3 const& v2) { \ + return !(v1 == v2); \ + } \ + CHD_FUNC bool operator==(T4 const& v1, T4 const& v2) { \ + return v1.x == v2.x && v1.y == v2.y && v1.z == v2.z && v1.w == v2.w; \ + } \ + CHD_FUNC bool operator!=(T4 const& v1, T4 const& v2) { \ + return !(v1 == v2); \ + } + +// These apply for all types +#define OTHER_FUNC_ALL(T, T2, T3, T4) \ + CHD_FUNC bool all_less(T2 const& v1, T2 const& v2) { \ + return (v1.x < v2.x) && (v1.y < v2.y); \ + } \ + CHD_FUNC bool all_less_or_eq(T2 const& v1, T2 const& v2) { \ + return (v1.x <= v2.x) && (v1.y <= v2.y); \ + } \ + CHD_FUNC bool all_greater(T2 const& v1, T2 const& v2) { \ + return (v1.x > v2.x) && (v1.y > v2.y); \ + } \ + CHD_FUNC bool all_greater_or_eq(T2 const& v1, T2 const& v2) { \ + return (v1.x >= v2.x) && (v1.y >= v2.y); \ + } \ + CHD_FUNC bool all_less(T3 const& v1, T3 const& v2) { \ + return (v1.x < v2.x) && (v1.y < v2.y) && (v1.z < v2.z); \ + } \ + CHD_FUNC bool all_less_or_eq(T3 const& v1, T3 const& v2) { \ + return (v1.x <= v2.x) && (v1.y <= v2.y) && (v1.z <= v2.z); \ + } \ + CHD_FUNC bool all_greater(T3 const& v1, T3 const& v2) { \ + return (v1.x > v2.x) && (v1.y > v2.y) && (v1.z > v2.z); \ + } \ + CHD_FUNC bool all_greater_or_eq(T3 const& v1, T3 const& v2) { \ + return (v1.x >= v2.x) && (v1.y >= v2.y) && (v1.z >= v2.z); \ + } \ + CHD_FUNC bool all_less(T4 const& v1, T4 const& v2) { \ + return (v1.x < v2.x) && (v1.y < v2.y) && (v1.z < v2.z) && (v1.w < v2.w); \ + } \ + CHD_FUNC bool all_less_or_eq(T4 const& v1, T4 const& v2) { \ + return (v1.x <= v2.x) && (v1.y <= v2.y) && (v1.z <= v2.z) && (v1.w <= v2.w); \ + } \ + CHD_FUNC bool all_greater(T4 const& v1, T4 const& v2) { \ + return (v1.x > v2.x) && (v1.y > v2.y) && (v1.z > v2.z) && (v1.w > v2.w); \ + } \ + CHD_FUNC bool all_greater_or_eq(T4 const& v1, T4 const& v2) { \ + return (v1.x >= v2.x) && (v1.y >= v2.y) && (v1.z >= v2.z) && (v1.w >= v2.w); \ + } \ + CHD_FUNC bool any_less(T2 const& v1, T2 const& v2) { \ + return (v1.x < v2.x) || (v1.y < v2.y); \ + } \ + CHD_FUNC bool any_less_or_eq(T2 const& v1, T2 const& v2) { \ + return (v1.x <= v2.x) || (v1.y <= v2.y); \ + } \ + CHD_FUNC bool any_greater(T2 const& v1, T2 const& v2) { \ + return (v1.x > v2.x) || (v1.y > v2.y); \ + } \ + CHD_FUNC bool any_greater_or_eq(T2 const& v1, T2 const& v2) { \ + return (v1.x >= v2.x) || (v1.y >= v2.y); \ + } \ + CHD_FUNC bool any_less(T3 const& v1, T3 const& v2) { \ + return (v1.x < v2.x) || (v1.y < v2.y) || (v1.z < v2.z); \ + } \ + CHD_FUNC bool any_less_or_eq(T3 const& v1, T3 const& v2) { \ + return (v1.x <= v2.x) || (v1.y <= v2.y) || (v1.z <= v2.z); \ + } \ + CHD_FUNC bool any_greater(T3 const& v1, T3 const& v2) { \ + return (v1.x > v2.x) || (v1.y > v2.y) || (v1.z > v2.z); \ + } \ + CHD_FUNC bool any_greater_or_eq(T3 const& v1, T3 const& v2) { \ + return (v1.x >= v2.x) || (v1.y >= v2.y) || (v1.z >= v2.z); \ + } \ + CHD_FUNC bool any_less(T4 const& v1, T4 const& v2) { \ + return (v1.x < v2.x) || (v1.y < v2.y) || (v1.z < v2.z) || (v1.w < v2.w); \ + } \ + CHD_FUNC bool any_less_or_eq(T4 const& v1, T4 const& v2) { \ + return (v1.x <= v2.x) || (v1.y <= v2.y) || (v1.z <= v2.z) || (v1.w <= v2.w); \ + } \ + CHD_FUNC bool any_greater(T4 const& v1, T4 const& v2) { \ + return (v1.x > v2.x) || (v1.y > v2.y) || (v1.z > v2.z) || (v1.w > v2.w); \ + } \ + CHD_FUNC bool any_greater_or_eq(T4 const& v1, T4 const& v2) { \ + return (v1.x >= v2.x) || (v1.y >= v2.y) || (v1.z >= v2.z) || (v1.w >= v2.w); \ + } \ + HD_FUNC T2 max(T2 const& v1, T const& v2) { \ + return {max(v1.x, v2), max(v1.y, v2)}; \ + } \ + HD_FUNC T2 max(T2 const& v1, T2 const& v2) { \ + return {max(v1.x, v2.x), max(v1.y, v2.y)}; \ + } \ + HD_FUNC T2 min(T2 const& v1, T const& v2) { \ + return {min(v1.x, v2), min(v1.y, v2)}; \ + } \ + HD_FUNC T2 min(T2 const& v1, T2 const& v2) { \ + return {min(v1.x, v2.x), min(v1.y, v2.y)}; \ + } \ + HD_FUNC T3 max(T3 const& v1, T const& v2) { \ + return {max(v1.x, v2), max(v1.y, v2), max(v1.z, v2)}; \ + } \ + HD_FUNC T3 max(T3 const& v1, T3 const& v2) { \ + return {max(v1.x, v2.x), max(v1.y, v2.y), max(v1.z, v2.z)}; \ + } \ + HD_FUNC T3 min(T3 const& v1, T const& v2) { \ + return {min(v1.x, v2), min(v1.y, v2), min(v1.z, v2)}; \ + } \ + HD_FUNC T3 min(T3 const& v1, T3 const& v2) { \ + return {min(v1.x, v2.x), min(v1.y, v2.y), min(v1.z, v2.z)}; \ + } \ + HD_FUNC T4 max(T4 const& v1, T const& v2) { \ + return {max(v1.x, v2), max(v1.y, v2), max(v1.z, v2), max(v1.w, v2)}; \ + } \ + HD_FUNC T4 max(T4 const& v1, T4 const& v2) { \ + return {max(v1.x, v2.x), max(v1.y, v2.y), max(v1.z, v2.z), max(v1.w, v2.w)}; \ + } \ + HD_FUNC T4 min(T4 const& v1, T const& v2) { \ + return {min(v1.x, v2), min(v1.y, v2), min(v1.z, v2), min(v1.w, v2)}; \ + } \ + HD_FUNC T4 min(T4 const& v1, T4 const& v2) { \ + return {min(v1.x, v2.x), min(v1.y, v2.y), min(v1.z, v2.z), min(v1.w, v2.w)}; \ + } \ + HD_FUNC T clamp(T v, T _min, T _max) { \ + return min(max(v, _min), _max); \ + } \ + HD_FUNC T2 clamp(T2 v, T2 _min, T2 _max) { \ + return min(max(v, _min), _max); \ + } \ + HD_FUNC T3 clamp(T3 v, T3 _min, T3 _max) { \ + return min(max(v, _min), _max); \ + } \ + HD_FUNC T4 clamp(T4 v, T4 _min, T4 _max) { \ + return min(max(v, _min), _max); \ + } \ + HD_FUNC T2 clamp(T2 v, T _min, T _max) { \ + return min(max(v, _min), _max); \ + } \ + HD_FUNC T3 clamp(T3 v, T _min, T _max) { \ + return min(max(v, _min), _max); \ + } \ + HD_FUNC T4 clamp(T4 v, T _min, T _max) { \ + return min(max(v, _min), _max); \ + } \ + CHD_FUNC T mix(T v1, T v2, bool a) { \ + return a ? v2 : v1; \ + } \ + CHD_FUNC T2 mix(T2 v1, T2 v2, bool a) { \ + return a ? v2 : v1; \ + } \ + CHD_FUNC T3 mix(T3 v1, T3 v2, bool a) { \ + return a ? v2 : v1; \ + } \ + CHD_FUNC T4 mix(T4 v1, T4 v2, bool a) { \ + return a ? v2 : v1; \ + } + +// These apply for all types, but unsigned ones +#define ABS_FUNC(T, T2, T3, T4) \ + HD_FUNC T2 abs(T2 const& v) { \ + return {abs(v.x), abs(v.y)}; \ + } \ + HD_FUNC T3 abs(T3 const& v) { \ + return {abs(v.x), abs(v.y), abs(v.z)}; \ + } \ + HD_FUNC T4 abs(T4 const& v) { \ + return {abs(v.x), abs(v.y), abs(v.z), abs(v.w)}; \ + } + +// Make functions +#define MAKE_FUNC(T, T2, T3, T4) \ + HD_FUNC T2 make_##T2(T scalar) { \ + return {scalar, scalar}; \ + } \ + HD_FUNC T3 make_##T3(T scalar) { \ + return {scalar, scalar, scalar}; \ + } \ + HD_FUNC T4 make_##T4(T scalar) { \ + return {scalar, scalar, scalar, scalar}; \ + } \ + HD_FUNC T3 make_##T3(T2 const& v, T scalar) { \ + return {v.x, v.y, scalar}; \ + } \ + HD_FUNC T3 make_##T3(T scalar, T2 const& v) { \ + return {scalar, v.x, v.y}; \ + } \ + HD_FUNC T4 make_##T4(T2 const& v1, T2 const& v2) { \ + return {v1.x, v1.y, v2.x, v2.y}; \ + } \ + HD_FUNC T4 make_##T4(T2 const& v, T scalar1, T scalar2) { \ + return {v.x, v.y, scalar1, scalar2}; \ + } \ + HD_FUNC T4 make_##T4(T scalar1, T scalar2, T2 const& v) { \ + return {scalar1, scalar2, v.x, v.y}; \ + } \ + HD_FUNC T4 make_##T4(T scalar1, T2 const& v, T scalar2) { \ + return {scalar1, v.x, v.y, scalar2}; \ + } \ + HD_FUNC T4 make_##T4(T3 const& v, T scalar) { \ + return {v.x, v.y, v.z, scalar}; \ + } \ + HD_FUNC T4 make_##T4(T scalar, T3 const& v) { \ + return {scalar, v.x, v.y, v.z}; \ + } \ + HD_FUNC T2 make_##T2(T3 const& v) { \ + return {v.x, v.y}; \ + } \ + HD_FUNC T2 make_##T2(T4 const& v) { \ + return {v.x, v.y}; \ + } \ + HD_FUNC T3 make_##T3(T4 const& v) { \ + return {v.x, v.y, v.z}; \ + } + +#define OTHER_FUNC_INT(T, T2, T3, T4) \ + CHD_FUNC T floor_div(T a, T b) { \ + T t = 1 - a / b; \ + return (a + t * b) / b - t; \ + } \ + CHD_FUNC T2 floor_div(T2 const& v1, T2 const& v2) { \ + return {floor_div(v1.x, v2.x), floor_div(v1.y, v2.y)}; \ + } \ + CHD_FUNC T2 floor_div(T2 const& v1, T v2) { \ + return {floor_div(v1.x, v2), floor_div(v1.y, v2)}; \ + } \ + CHD_FUNC T3 floor_div(T3 const& v1, T3 const& v2) { \ + return {floor_div(v1.x, v2.x), floor_div(v1.y, v2.y), floor_div(v1.z, v2.z)}; \ + } \ + CHD_FUNC T3 floor_div(T3 const& v1, T v2) { \ + return {floor_div(v1.x, v2), floor_div(v1.y, v2), floor_div(v1.z, v2)}; \ + } \ + CHD_FUNC T4 floor_div(T4 const& v1, T4 const& v2) { \ + return { \ + floor_div(v1.x, v2.x), \ + floor_div(v1.y, v2.y), \ + floor_div(v1.z, v2.z), \ + floor_div(v1.w, v2.w)}; \ + } \ + CHD_FUNC T4 floor_div(T4 const& v1, T v2) { \ + return {floor_div(v1.x, v2), floor_div(v1.y, v2), floor_div(v1.z, v2), floor_div(v1.w, v2)}; \ + } + +#define OTHER_FUNC_FP(T, T2, T3, T4) \ + CHD_FUNC T dot(T2 a, T2 b) { \ + return a.x * b.x + a.y * b.y; \ + } \ + CHD_FUNC T dot(T3 a, T3 b) { \ + return a.x * b.x + a.y * b.y + a.z * b.z; \ + } \ + CHD_FUNC T dot(T4 a, T4 b) { \ + return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; \ + } \ + CHD_FUNC T cross(T2 a, T2 b) { \ + return a.x * b.y - a.y * b.x; \ + } \ + CHD_FUNC T3 cross(T3 a, T3 b) { \ + return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; \ + } \ + HD_FUNC T norm(T2 a) { \ + return sqrt(dot(a, a)); \ + } \ + HD_FUNC T norm(T3 a) { \ + return detail::norm3d(a.x, a.y, a.z); \ + } \ + HD_FUNC T norm(T4 a) { \ + return detail::norm4d(a.x, a.y, a.z, a.w); \ + } \ + HD_FUNC T rnorm(T2 a) { \ + return rsqrt(dot(a, a)); \ + } \ + HD_FUNC T rnorm(T3 a) { \ + return detail::rnorm3d(a.x, a.y, a.z); \ + } \ + HD_FUNC T rnorm(T4 a) { \ + return detail::rnorm4d(a.x, a.y, a.z, a.w); \ + } \ + HD_FUNC T2 normalize(T2 v) { \ + T invLen = rnorm(v); \ + return v * invLen; \ + } \ + HD_FUNC T3 normalize(T3 v) { \ + T invLen = rnorm(v); \ + return v * invLen; \ + } \ + HD_FUNC T4 normalize(T4 v) { \ + T invLen = rnorm(v); \ + return v * invLen; \ + } \ + HD_FUNC T2 saturate(T2 v) { \ + return {saturate(v.x), saturate(v.y)}; \ + } \ + HD_FUNC T3 saturate(T3 v) { \ + return {saturate(v.x), saturate(v.y), saturate(v.z)}; \ + } \ + HD_FUNC T4 saturate(T4 v) { \ + return {saturate(v.x), saturate(v.y), saturate(v.z), saturate(v.w)}; \ + } \ + CHD_FUNC T sign(T v) { \ + return v > 0 ? 1 : (v < 0 ? -1 : 0); \ + } \ + CHD_FUNC T2 sign(T2 v) { \ + return {sign(v.x), sign(v.y)}; \ + } \ + CHD_FUNC T3 sign(T3 v) { \ + return {sign(v.x), sign(v.y), sign(v.z)}; \ + } \ + CHD_FUNC T4 sign(T4 v) { \ + return {sign(v.x), sign(v.y), sign(v.z), sign(v.w)}; \ + } \ + CHD_FUNC T mix(T v1, T v2, T a) { \ + return v1 * (T(1.0) - a) + v2 * a; \ + } \ + CHD_FUNC T2 mix(T2 v1, T2 v2, T a) { \ + return v1 * (T(1.0) - a) + v2 * a; \ + } \ + CHD_FUNC T3 mix(T3 v1, T3 v2, T a) { \ + return v1 * (T(1.0) - a) + v2 * a; \ + } \ + CHD_FUNC T4 mix(T4 v1, T4 v2, T a) { \ + return v1 * (T(1.0) - a) + v2 * a; \ + } \ + CHD_FUNC T2 mix(T2 v1, T2 v2, T2 a) { \ + return v1 * (T(1.0) - a) + v2 * a; \ + } \ + CHD_FUNC T3 mix(T3 v1, T3 v2, T3 a) { \ + return v1 * (T(1.0) - a) + v2 * a; \ + } \ + CHD_FUNC T4 mix(T4 v1, T4 v2, T4 a) { \ + return v1 * (T(1.0) - a) + v2 * a; \ + } \ + CHD_FUNC T sum(T2 const& v) { \ + return v.x + v.y; \ + } \ + CHD_FUNC T sum(T3 const& v) { \ + return v.x + v.y + v.z; \ + } \ + CHD_FUNC T sum(T4 const& v) { \ + return v.x + v.y + v.z + v.w; \ + } \ + HD_FUNC T epsclamp(T v, T eps) { \ + return (v < 0) ? min(v, -eps) : max(v, eps); \ + } \ + HD_FUNC T epsclamp(T v) { \ + return epsclamp(v, epsilon::value); \ + } \ + HD_FUNC T2 epsclamp(T2 v) { \ + return {epsclamp(v.x), epsclamp(v.y)}; \ + } \ + HD_FUNC T3 epsclamp(T3 v) { \ + return {epsclamp(v.x), epsclamp(v.y), epsclamp(v.z)}; \ + } \ + HD_FUNC T4 epsclamp(T4 v) { \ + return {epsclamp(v.x), epsclamp(v.y), epsclamp(v.z), epsclamp(v.w)}; \ + } \ + HD_FUNC T2 epsclamp(T2 v, T eps) { \ + return {epsclamp(v.x, eps), epsclamp(v.y, eps)}; \ + } \ + HD_FUNC T3 epsclamp(T3 v, T eps) { \ + return {epsclamp(v.x, eps), epsclamp(v.y, eps), epsclamp(v.z, eps)}; \ + } \ + HD_FUNC T4 epsclamp(T4 v, T eps) { \ + return {epsclamp(v.x, eps), epsclamp(v.y, eps), epsclamp(v.z, eps), epsclamp(v.w, eps)}; \ + } \ + CHD_FUNC void inverse(const T2(&m)[2], T2(&out)[2]) { \ + T det_m = T(1.0) / (m[0].x * m[1].y - m[0].y * m[1].x); \ + out[0] = det_m * T2({m[1].y, -m[0].y}); \ + out[1] = det_m * T2({-m[1].x, m[0].x}); \ + } \ + CHD_FUNC void inverse(const T3(&m)[3], T3(&out)[3]) { \ + T det_m = T(1.0) / \ + (+m[0].x * (m[1].y * m[2].z - m[1].z * m[2].y) - \ + m[0].y * (m[1].x * m[2].z - m[1].z * m[2].x) + \ + m[0].z * (m[1].x * m[2].y - m[1].y * m[2].x)); \ + out[0] = det_m * \ + T3({ \ + +(m[1].y * m[2].z - m[2].y * m[1].z), \ + -(m[0].y * m[2].z - m[2].y * m[0].z), \ + +(m[0].y * m[1].z - m[1].y * m[0].z), \ + }); \ + out[1] = det_m * \ + T3({ \ + -(m[1].x * m[2].z - m[2].x * m[1].z), \ + +(m[0].x * m[2].z - m[2].x * m[0].z), \ + -(m[0].x * m[1].z - m[1].x * m[0].z), \ + }); \ + out[2] = det_m * \ + T3({ \ + +(m[1].x * m[2].y - m[2].x * m[1].y), \ + -(m[0].x * m[2].y - m[2].x * m[0].y), \ + +(m[0].x * m[1].y - m[1].x * m[0].y), \ + }); \ + } \ + CHD_FUNC T2 mul(const T2(&r)[2], T2 v) { \ + return T2({dot(r[0], v), dot(r[1], v)}); \ + } \ + CHD_FUNC T3 mul(const T3(&r)[3], T3 v) { \ + return T3({dot(r[0], v), dot(r[1], v), dot(r[2], v)}); \ + } \ + CHD_FUNC T4 mul(const T4(&r)[4], T4 v) { \ + return T4({dot(r[0], v), dot(r[1], v), dot(r[2], v), dot(r[3], v)}); \ + } \ + CHD_FUNC void mul(const T2(&a)[2], const T2(&b)[2], T2(&out)[2]) { \ + out[0] = T2({dot(a[0], T2({b[0].x, b[1].x})), dot(a[0], T2({b[0].y, b[1].y}))}); \ + out[1] = T2({dot(a[1], T2({b[0].x, b[1].x})), dot(a[1], T2({b[0].y, b[1].y}))}); \ + } \ + CHD_FUNC void mul(const T2(&a)[2], const T3(&b)[2], T3(&out)[2]) { \ + out[0] = \ + T3({dot(a[0], T2({b[0].x, b[1].x})), \ + dot(a[0], T2({b[0].y, b[1].y})), \ + dot(a[0], T2({b[0].z, b[1].z}))}); \ + out[1] = \ + T3({dot(a[1], T2({b[0].x, b[1].x})), \ + dot(a[1], T2({b[0].y, b[1].y})), \ + dot(a[1], T2({b[0].z, b[1].z}))}); \ + } \ + CHD_FUNC void mul(const T3(&a)[3], const T3(&b)[3], T3(&out)[3]) { \ + out[0] = \ + T3({dot(a[0], T3({b[0].x, b[1].x, b[2].x})), \ + dot(a[0], T3({b[0].y, b[1].y, b[2].y})), \ + dot(a[0], T3({b[0].z, b[1].z, b[2].z}))}); \ + out[1] = \ + T3({dot(a[1], T3({b[0].x, b[1].x, b[2].x})), \ + dot(a[1], T3({b[0].y, b[1].y, b[2].y})), \ + dot(a[1], T3({b[0].z, b[1].z, b[2].z}))}); \ + out[2] = \ + T3({dot(a[2], T3({b[0].x, b[1].x, b[2].x})), \ + dot(a[2], T3({b[0].y, b[1].y, b[2].y})), \ + dot(a[2], T3({b[0].z, b[1].z, b[2].z}))}); \ + } + +#define DEFINE_FUNC_FOR_UNSIGNED_INT(T, T2, T3, T4) \ + UNARY_OP(T, T2, T3, T4) \ + BINARY_ARITHM_OP(T, T2, T3, T4) \ + BINARY_BIT_OP(T, T2, T3, T4) \ + BINARY_EQ_OP(T, T2, T3, T4) \ + BINARY_INT_OP(T, T2, T3, T4) \ + OTHER_FUNC_ALL(T, T2, T3, T4) \ + OTHER_FUNC_INT(T, T2, T3, T4) \ + MAKE_FUNC(T, T2, T3, T4) + +#define DEFINE_FUNC_FOR_SIGNED_INT(T, T2, T3, T4) \ + DEFINE_FUNC_FOR_UNSIGNED_INT(T, T2, T3, T4) \ + ABS_FUNC(T, T2, T3, T4) + +#define DEFINE_FUNC_FOR_FLOAT(T, T2, T3, T4) \ + UNARY_OP(T, T2, T3, T4) \ + BINARY_ARITHM_OP(T, T2, T3, T4) \ + BINARY_EQ_OP(T, T2, T3, T4) \ + OTHER_FUNC_ALL(T, T2, T3, T4) \ + OTHER_FUNC_FP(T, T2, T3, T4) \ + ABS_FUNC(T, T2, T3, T4) \ + MAKE_FUNC(T, T2, T3, T4) + +DEFINE_FUNC_FOR_UNSIGNED_INT(unsigned int, uint2, uint3, uint4); +DEFINE_FUNC_FOR_SIGNED_INT(int, int2, int3, int4); +DEFINE_FUNC_FOR_FLOAT(float, float2, float3, float4); +DEFINE_FUNC_FOR_FLOAT(double, double2, double3, double4); + +namespace detail { +template +struct VecType; +} + +// Type inference utils for writing templates +// +// Derive vector type given the scalar type: +// math::TVec2 a; // `a` is of type `float2` +// math::TVec3 b; // `b` is of type `int3`; +// +// Derive vector type given the vector size and scalar type: +// math::TVec c; // `c` is of type `double4`; +template +using TVec1 = typename detail::VecType::scalar1_t; + +template +using TVec2 = typename detail::VecType::scalar2_t; + +template +using TVec3 = typename detail::VecType::scalar3_t; + +template +using TVec4 = typename detail::VecType::scalar4_t; + +template +using TVec = typename detail::VecType::template dim::type; + +namespace detail { +template class Vec, typename scalar_t> +struct VecD; + +template