crocodile (empty) → 0.1
raw patch · 40 files changed
+3256/−0 lines, 40 filesdep +HUnitdep +basedep +bmpsetup-changed
Dependencies added: HUnit, base, bmp, bytestring, deepseq, ghc-prim, heap, mersenne-random-pure64, mtl, parallel
Files
- LICENSE +339/−0
- README +72/−0
- Setup.hs +2/−0
- app/src/BoundingBox.hs +65/−0
- app/src/Camera.hs +22/−0
- app/src/Colour.hs +121/−0
- app/src/CornellBox.hs +162/−0
- app/src/Distribution.hs +55/−0
- app/src/IrradianceCache.hs +75/−0
- app/src/KDTree.hs +58/−0
- app/src/Light.hs +110/−0
- app/src/Light.hs-boot +29/−0
- app/src/Main.hs +148/−0
- app/src/Material.hs +21/−0
- app/src/Matrix.hs +55/−0
- app/src/Misc.hs +31/−0
- app/src/Octree.hs +133/−0
- app/src/PhotonMap.hs +261/−0
- app/src/PhotonMap.hs-boot +28/−0
- app/src/Primitive.hs +298/−0
- app/src/Primitive.hs-boot +69/−0
- app/src/Ray.hs +32/−0
- app/src/RayTrace.hs +254/−0
- app/src/RayTrace.hs-boot +22/−0
- app/src/RenderContext.hs +23/−0
- app/src/RenderContext.hs-boot +20/−0
- app/src/SceneGraph.hs +48/−0
- app/src/SceneGraphTest.hs +30/−0
- app/src/Shader.hs +50/−0
- app/src/Tests/BoundingBoxTest.hs +22/−0
- app/src/Tests/ColourTest.hs +84/−0
- app/src/Tests/UnitTests.hs +10/−0
- app/src/Tests/VectorTest.hs +144/−0
- app/src/ToneMap.hs +61/−0
- app/src/Vector.hs +233/−0
- app/src/scripts/build +4/−0
- app/src/scripts/clean +2/−0
- app/src/scripts/run-hlint +3/−0
- app/src/scripts/unit_test +2/−0
- crocodile.cabal +58/−0
+ LICENSE view
@@ -0,0 +1,339 @@+ GNU GENERAL PUBLIC LICENSE+ Version 2, June 1991++ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++ Preamble++ The licenses for most software are designed to take away your+freedom to share and change it. By contrast, the GNU General Public+License is intended to guarantee your freedom to share and change free+software--to make sure the software is free for all its users. This+General Public License applies to most of the Free Software+Foundation's software and to any other program whose authors commit to+using it. (Some other Free Software Foundation software is covered by+the GNU Lesser General Public License instead.) You can apply it to+your programs, too.++ When we speak of free software, we are referring to freedom, not+price. Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+this service if you wish), that you receive source code or can get it+if you want it, that you can change the software or use pieces of it+in new free programs; and that you know you can do these things.++ To protect your rights, we need to make restrictions that forbid+anyone to deny you these rights or to ask you to surrender the rights.+These restrictions translate to certain responsibilities for you if you+distribute copies of the software, or if you modify it.++ For example, if you distribute copies of such a program, whether+gratis or for a fee, you must give the recipients all the rights that+you have. You must make sure that they, too, receive or can get the+source code. And you must show them these terms so they know their+rights.++ We protect your rights with two steps: (1) copyright the software, and+(2) offer you this license which gives you legal permission to copy,+distribute and/or modify the software.++ Also, for each author's protection and ours, we want to make certain+that everyone understands that there is no warranty for this free+software. If the software is modified by someone else and passed on, we+want its recipients to know that what they have is not the original, so+that any problems introduced by others will not reflect on the original+authors' reputations.++ Finally, any free program is threatened constantly by software+patents. We wish to avoid the danger that redistributors of a free+program will individually obtain patent licenses, in effect making the+program proprietary. To prevent this, we have made it clear that any+patent must be licensed for everyone's free use or not licensed at all.++ The precise terms and conditions for copying, distribution and+modification follow.++ GNU GENERAL PUBLIC LICENSE+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION++ 0. This License applies to any program or other work which contains+a notice placed by the copyright holder saying it may be distributed+under the terms of this General Public License. The "Program", below,+refers to any such program or work, and a "work based on the Program"+means either the Program or any derivative work under copyright law:+that is to say, a work containing the Program or a portion of it,+either verbatim or with modifications and/or translated into another+language. (Hereinafter, translation is included without limitation in+the term "modification".) Each licensee is addressed as "you".++Activities other than copying, distribution and modification are not+covered by this License; they are outside its scope. The act of+running the Program is not restricted, and the output from the Program+is covered only if its contents constitute a work based on the+Program (independent of having been made by running the Program).+Whether that is true depends on what the Program does.++ 1. You may copy and distribute verbatim copies of the Program's+source code as you receive it, in any medium, provided that you+conspicuously and appropriately publish on each copy an appropriate+copyright notice and disclaimer of warranty; keep intact all the+notices that refer to this License and to the absence of any warranty;+and give any other recipients of the Program a copy of this License+along with the Program.++You may charge a fee for the physical act of transferring a copy, and+you may at your option offer warranty protection in exchange for a fee.++ 2. You may modify your copy or copies of the Program or any portion+of it, thus forming a work based on the Program, and copy and+distribute such modifications or work under the terms of Section 1+above, provided that you also meet all of these conditions:++ a) You must cause the modified files to carry prominent notices+ stating that you changed the files and the date of any change.++ b) You must cause any work that you distribute or publish, that in+ whole or in part contains or is derived from the Program or any+ part thereof, to be licensed as a whole at no charge to all third+ parties under the terms of this License.++ c) If the modified program normally reads commands interactively+ when run, you must cause it, when started running for such+ interactive use in the most ordinary way, to print or display an+ announcement including an appropriate copyright notice and a+ notice that there is no warranty (or else, saying that you provide+ a warranty) and that users may redistribute the program under+ these conditions, and telling the user how to view a copy of this+ License. (Exception: if the Program itself is interactive but+ does not normally print such an announcement, your work based on+ the Program is not required to print an announcement.)++These requirements apply to the modified work as a whole. If+identifiable sections of that work are not derived from the Program,+and can be reasonably considered independent and separate works in+themselves, then this License, and its terms, do not apply to those+sections when you distribute them as separate works. But when you+distribute the same sections as part of a whole which is a work based+on the Program, the distribution of the whole must be on the terms of+this License, whose permissions for other licensees extend to the+entire whole, and thus to each and every part regardless of who wrote it.++Thus, it is not the intent of this section to claim rights or contest+your rights to work written entirely by you; rather, the intent is to+exercise the right to control the distribution of derivative or+collective works based on the Program.++In addition, mere aggregation of another work not based on the Program+with the Program (or with a work based on the Program) on a volume of+a storage or distribution medium does not bring the other work under+the scope of this License.++ 3. You may copy and distribute the Program (or a work based on it,+under Section 2) in object code or executable form under the terms of+Sections 1 and 2 above provided that you also do one of the following:++ a) Accompany it with the complete corresponding machine-readable+ source code, which must be distributed under the terms of Sections+ 1 and 2 above on a medium customarily used for software interchange; or,++ b) Accompany it with a written offer, valid for at least three+ years, to give any third party, for a charge no more than your+ cost of physically performing source distribution, a complete+ machine-readable copy of the corresponding source code, to be+ distributed under the terms of Sections 1 and 2 above on a medium+ customarily used for software interchange; or,++ c) Accompany it with the information you received as to the offer+ to distribute corresponding source code. (This alternative is+ allowed only for noncommercial distribution and only if you+ received the program in object code or executable form with such+ an offer, in accord with Subsection b above.)++The source code for a work means the preferred form of the work for+making modifications to it. For an executable work, complete source+code means all the source code for all modules it contains, plus any+associated interface definition files, plus the scripts used to+control compilation and installation of the executable. However, as a+special exception, the source code distributed need not include+anything that is normally distributed (in either source or binary+form) with the major components (compiler, kernel, and so on) of the+operating system on which the executable runs, unless that component+itself accompanies the executable.++If distribution of executable or object code is made by offering+access to copy from a designated place, then offering equivalent+access to copy the source code from the same place counts as+distribution of the source code, even though third parties are not+compelled to copy the source along with the object code.++ 4. You may not copy, modify, sublicense, or distribute the Program+except as expressly provided under this License. Any attempt+otherwise to copy, modify, sublicense or distribute the Program is+void, and will automatically terminate your rights under this License.+However, parties who have received copies, or rights, from you under+this License will not have their licenses terminated so long as such+parties remain in full compliance.++ 5. You are not required to accept this License, since you have not+signed it. However, nothing else grants you permission to modify or+distribute the Program or its derivative works. These actions are+prohibited by law if you do not accept this License. Therefore, by+modifying or distributing the Program (or any work based on the+Program), you indicate your acceptance of this License to do so, and+all its terms and conditions for copying, distributing or modifying+the Program or works based on it.++ 6. Each time you redistribute the Program (or any work based on the+Program), the recipient automatically receives a license from the+original licensor to copy, distribute or modify the Program subject to+these terms and conditions. You may not impose any further+restrictions on the recipients' exercise of the rights granted herein.+You are not responsible for enforcing compliance by third parties to+this License.++ 7. If, as a consequence of a court judgment or allegation of patent+infringement or for any other reason (not limited to patent issues),+conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License. If you cannot+distribute so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you+may not distribute the Program at all. For example, if a patent+license would not permit royalty-free redistribution of the Program by+all those who receive copies directly or indirectly through you, then+the only way you could satisfy both it and this License would be to+refrain entirely from distribution of the Program.++If any portion of this section is held invalid or unenforceable under+any particular circumstance, the balance of the section is intended to+apply and the section as a whole is intended to apply in other+circumstances.++It is not the purpose of this section to induce you to infringe any+patents or other property right claims or to contest validity of any+such claims; this section has the sole purpose of protecting the+integrity of the free software distribution system, which is+implemented by public license practices. Many people have made+generous contributions to the wide range of software distributed+through that system in reliance on consistent application of that+system; it is up to the author/donor to decide if he or she is willing+to distribute software through any other system and a licensee cannot+impose that choice.++This section is intended to make thoroughly clear what is believed to+be a consequence of the rest of this License.++ 8. If the distribution and/or use of the Program is restricted in+certain countries either by patents or by copyrighted interfaces, the+original copyright holder who places the Program under this License+may add an explicit geographical distribution limitation excluding+those countries, so that distribution is permitted only in or among+countries not thus excluded. In such case, this License incorporates+the limitation as if written in the body of this License.++ 9. The Free Software Foundation may publish revised and/or new versions+of the General Public License from time to time. Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++Each version is given a distinguishing version number. If the Program+specifies a version number of this License which applies to it and "any+later version", you have the option of following the terms and conditions+either of that version or of any later version published by the Free+Software Foundation. If the Program does not specify a version number of+this License, you may choose any version ever published by the Free Software+Foundation.++ 10. If you wish to incorporate parts of the Program into other free+programs whose distribution conditions are different, write to the author+to ask for permission. For software which is copyrighted by the Free+Software Foundation, write to the Free Software Foundation; we sometimes+make exceptions for this. Our decision will be guided by the two goals+of preserving the free status of all derivatives of our free software and+of promoting the sharing and reuse of software generally.++ NO WARRANTY++ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,+REPAIR OR CORRECTION.++ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE+POSSIBILITY OF SUCH DAMAGES.++ END OF TERMS AND CONDITIONS++ How to Apply These Terms to Your New Programs++ If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++ To do so, attach the following notices to the program. It is safest+to attach them to the start of each source file to most effectively+convey the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ This program is free software; you can redistribute it and/or modify+ it under the terms of the GNU General Public License as published by+ the Free Software Foundation; either version 2 of the License, or+ (at your option) any later version.++ This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License for more details.++ You should have received a copy of the GNU General Public License along+ with this program; if not, write to the Free Software Foundation, Inc.,+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.++Also add information on how to contact you by electronic and paper mail.++If the program is interactive, make it output a short notice like this+when it starts in an interactive mode:++ Gnomovision version 69, Copyright (C) year name of author+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+ This is free software, and you are welcome to redistribute it+ under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License. Of course, the commands you use may+be called something other than `show w' and `show c'; they could even be+mouse-clicks or menu items--whatever suits your program.++You should also get your employer (if you work as a programmer) or your+school, if any, to sign a "copyright disclaimer" for the program, if+necessary. Here is a sample; alter the names:++ Yoyodyne, Inc., hereby disclaims all copyright interest in the program+ `Gnomovision' (which makes passes at compilers) written by James Hacker.++ <signature of Ty Coon>, 1 April 1989+ Ty Coon, President of Vice++This General Public License does not permit incorporating your program into+proprietary programs. If your program is a subroutine library, you may+consider it more useful to permit linking proprietary applications with the+library. If this is what you want to do, use the GNU Lesser General+Public License instead of this License.
+ README view
@@ -0,0 +1,72 @@+Haskell Ray Tracer and Photon Mapper v0.0+-----------------------------------------++This is a parallel ray tracer and partially parallel photon mapper written in Haskell. ++Features:++Fully recursive raytracer with reflection and refraction+* Sphere primitive+* Plane primitive+* Triangle mesh primitive+* Parallelised ray tracing+* Distributed ray tracing giving depth of field and anti-aliasing+* Photon mapping+* Irradiance caching+* Gamma correction+* Tone mapping++Disclaimer+----------++I am a novice Haskell programmer. This code is not intended to be representative of best-practice Haskell programming.++If you see code that could be improved in terms of style, correctness, clarity, flexibility or efficiency, please, let me know! I'm eager to learn and I'd love to hear the opinion of those far more capable than me.++Why call it "crocodile"?+------------------------++My two-year old son currently has quite a penchant for crocodiles.++Usage+-----++Currently there is a default hard-coded scene of the Cornell Box. I build the program with:++scripts/build++and execute it with:++time ./crocodile -p -i -c +RTS -N -RTS++Cabal also works, of course++This will output a file called test.bmp containing the resulting image. Depending on the number of photons emitted, the photon gathering radius and the maximum number of photons gathered, this could take quite some time. You can tune the parameters in Main.hs.++Options+-------++-i outputs intermediate renderings for fast debug feedback (recommended)+-p enables the photon mapping pass+-v directly visualises the photon map+-c enable the irradiance cache++Bugs+----++The Photon Mapping is largely functional though it currently has some minor artefacts. Work is ongoing to fix them.+There are numerous TODO issues noted in the code.++Future planned work+-------------------++* Gradients for the irradiance cache+* Parsing of scene data files+* Optimisation! Particularly of photon mapping+* Parallelisation of photon gathering+* Extended shader model+* Alternative GI code paths - e.g. path tracing+* Caustic photon mapping++Tom Hammersley 12/5/2011+tomhammersley@gmail.com
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/src/BoundingBox.hs view
@@ -0,0 +1,65 @@+-- Bounding box code+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE BangPatterns #-}++module BoundingBox where++import Vector+import GHC.Types+import GHC.Prim++type AABB = (Vector, Vector)++boundingBoxUnion :: AABB -> AABB -> AABB+boundingBoxUnion (min1, max1) (min2, max2) = (Vector.min min1 min2, Vector.max max1 max2)++boundingBoxValid :: AABB -> Bool+boundingBoxValid (boxMin, boxMax) = vecX boxMin <= vecX boxMax &&+ vecY boxMin <= vecY boxMax &&+ vecZ boxMin <= vecZ boxMax++boundingBoxOverlaps :: AABB -> AABB -> Bool+boundingBoxOverlaps box1 box2 = overlaps box1 box2 || overlaps box2 box1+ where+ overlaps (min1, max1) (min2, max2) = vecX min1 <= vecX max2 && (vecX max1 >= vecX min2) &&+ vecY min1 <= vecY max2 && (vecY max1 >= vecY min2) &&+ vecZ min1 <= vecZ max2 && (vecZ max1 >= vecZ min2)++-- Enlarge a bounding box to include a point+enlargeBoundingBox :: Position -> AABB -> AABB+enlargeBoundingBox pos(boxMin, boxMax) = (Vector.min boxMin pos, Vector.max boxMax pos)++-- Linearly scale a box+scaleBoundingBox :: AABB -> Double -> AABB+scaleBoundingBox (boxMin, boxMax) k = (setWTo1 $ boxMin <*> k, setWTo1 $ boxMax <*> k)++-- Give a bounding box a buffer of a certain distance all the way around+growBoundingBox :: AABB -> Double -> AABB+growBoundingBox (Vector x1 y1 z1 _, Vector x2 y2 z2 _) k = (Vector (x1 - k) (y1 - k) (z1 - k) 1, Vector (x2 + k) (y2 + k) (z2 + k) 1)++-- This is a dummy box that is used initially. If anything is intersected with it, it becomes valid. Else it is an invalid box that can be tested for+initialInvalidBox :: AABB+initialInvalidBox = (Vector bigNumber bigNumber bigNumber 1, Vector smallNumber smallNumber smallNumber 1)+ where+ !bigNumber = 10000000+ !smallNumber = -10000000++-- These functions are useful for finding the greatest or smallest part of a box relative to a normal+selectMinBoxComponent :: (Vector -> Double) -> Vector -> AABB -> Double+selectMinBoxComponent f norm (boxMin, boxMax) = if f norm > 0 then f boxMin else f boxMax++selectMaxBoxComponent :: (Vector -> Double) -> Vector -> AABB -> Double+selectMaxBoxComponent f norm (boxMin, boxMax) = if f norm > 0 then f boxMax else f boxMin++-- Does a box contain a point?+contains :: AABB -> Position -> Bool+{-# SPECIALIZE INLINE contains :: AABB -> Position -> Bool #-}+contains (Vector !(D# minX) !(D# minY) !(D# minZ) _, Vector !(D# maxX) !(D# maxY) !(D# maxZ) _) (Vector !(D# x) !(D# y) !(D# z) _) = + x >=## minX && x <=## maxX &&+ y >=## minY && y <=## maxY &&+ z >=## minZ && z <=## maxZ++overlapsSphere :: AABB -> Position -> Double -> Bool+overlapsSphere (boxMin, boxMax) p r = all insideInterval [vecX, vecY, vecZ]+ where + insideInterval f = f p >= (f boxMin - r) && f p <= (f boxMax + r)
+ app/src/Camera.hs view
@@ -0,0 +1,22 @@+-- Camera module++module Camera where++import Vector+import Matrix++data Camera = Camera { worldToCamera :: Matrix, fieldOfView :: !Double, position :: Vector } deriving (Show)++lookAt :: Position -> Position -> Direction -> Double -> Camera+lookAt pos target up fov = + Camera matrix fov pos+ where+ forward = normalise $ target - pos+ right = normalise $ up `cross` forward+ up' = right `cross` forward+ matrix = buildMatrix right up' forward pos++withVectors :: Position -> Direction -> Direction -> Direction -> Double -> Camera+withVectors pos basisX basisY basisZ fov = Camera matrix fov pos+ where+ matrix = buildMatrix basisX basisY basisZ (Vector.negate pos)
+ app/src/Colour.hs view
@@ -0,0 +1,121 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module Colour where+import Vector hiding (min, max)+import Misc+import Data.Word+import Control.DeepSeq++-- Normalised RGBA colour+data Colour = Colour { red :: {-# UNPACK #-} !Double, + green :: {-# UNPACK #-} !Double, + blue :: {-# UNPACK #-} !Double, + alpha :: {-# UNPACK #-} !Double } deriving (Show, Read, Ord, Eq)++instance Num Colour where+ {-# SPECIALIZE INLINE (+) :: Colour -> Colour -> Colour #-}+ (Colour !r1 !g1 !b1 !a1) + (Colour !r2 !g2 !b2 !a2) = Colour (r1 + r2) (g1 + g2) (b1 + b2) (a1 + a2)+ {-# SPECIALIZE INLINE (-) :: Colour -> Colour -> Colour #-}+ (Colour !r1 !g1 !b1 !a1) - (Colour !r2 !g2 !b2 !a2) = Colour (r1 - r2) (g1 - g2) (b1 - b2) (a1 - a2)+ {-# SPECIALIZE INLINE (*) :: Colour -> Colour -> Colour #-}+ (Colour !r1 !g1 !b1 !a1) * (Colour !r2 !g2 !b2 !a2) = Colour (r1 * r2) (g1 * g2) (b1 * b2) (a1 * a2)+ abs (Colour r g b a) = Colour (abs r) (abs g) (abs b) (abs a)+ signum (Colour r g b a) = Colour (signum r) (signum g) (signum b) (signum a)+ fromInteger x = Colour (fromInteger x) (fromInteger x) (fromInteger x) (fromInteger x)++instance Fractional Colour where+ {-# SPECIALIZE INLINE (/) :: Colour -> Colour -> Colour #-}+ (Colour !r1 !g1 !b1 !a1) / (Colour !r2 !g2 !b2 !a2) = Colour (r1 / r2) (g1 / g2) (b1 / b2) (a1 / a2)+ fromRational x = Colour (fromRational x) (fromRational x) (fromRational x) (fromRational x)++instance NFData Colour where+ rnf (Colour r g b a) = rnf r `seq` rnf g `seq` rnf b `seq` rnf a++infixl 7 <*>+infixl 7 </>+infixl 6 <+>+infixl 6 <->++(<*>) :: Colour -> Double -> Colour+(Colour !r !g !b !a) <*> k = Colour (r * k) (g * k) (b * k) (a * k)++(</>) :: Colour -> Double -> Colour+(Colour !r !g !b !a) </> k = Colour (r / k) (g / k) (b / k) (a / k)++(<+>) :: Colour -> Double -> Colour+(Colour !r !g !b !a) <+> k = Colour (r + k) (g + k) (b + k) (a + k)++(<->) :: Colour -> Double -> Colour+(Colour !r !g !b !a) <-> k = Colour (r - k) (g - k) (b - k) (a - k)++clamp :: Colour -> Colour+clamp (Colour !r !g !b !a) = Colour (max 0 (min r 1)) (max 0 (min g 1)) (max 0 (min b 1)) (max 0 (min a 1))++fold :: (Double -> Double -> Double) -> Colour -> Double -> Colour+fold f (Colour !r !g !b !a) k = Colour (f r k) (f g k) (f b k) (f a k)++-- Basic colours+colRed :: Colour+colRed = Colour 1 0 0 1++colGreen :: Colour+colGreen = Colour 0 1 0 1++colBlue :: Colour+colBlue = Colour 0 0 1 1++colWhite :: Colour+colWhite = Colour 1 1 1 1++colBlack :: Colour+colBlack = Colour 0 0 0 1++colGrey :: Colour+colGrey = Colour 0.5 0.5 0.5 1++colYellow :: Colour+colYellow = Colour 1 1 0 1++gamma :: Double+gamma = 2.2++invGamma :: Double+invGamma = 1.0 / gamma++-- Gamma correct a colour+gammaCorrect :: Colour -> Colour+gammaCorrect (Colour !r !g !b !a) = Colour (r ** gamma) (g ** gamma) (b ** gamma) (a ** gamma)++invGammaCorrect ::Colour -> Colour+invGammaCorrect (Colour !r !g !b !a) = Colour (r ** invGamma) (g ** invGamma) (b ** invGamma) (a ** invGamma)++-- Colour encode a normal+encodeNormal :: Vector -> Colour+encodeNormal (Vector !x !y !z _) = gammaCorrect $ Colour (saturate $ x * 0.5 + 0.5) (saturate $ y * 0.5 + 0.5) (saturate $ z * 0.5 + 0.5) 1++-- Convert a list of colours to a list of Word8s+convertColoursToPixels :: [Colour] -> [Word8]+convertColoursToPixels (col:cols) = r : g : b : 255 : convertColoursToPixels cols+ where+ r = truncate (red col * 255.0)+ g = truncate (green col * 255.0)+ b = truncate (blue col * 255.0)+convertColoursToPixels [] = []++-- Measure overall magnitude of a colour+magnitude :: Colour -> Double+magnitude (Colour r g b _) = r * 0.3 + g * 0.6 + b * 0.1++-- Convert to a list+toListRGBA :: Colour -> [Double]+toListRGBA (Colour r g b a) = [r, g, b, a]++toListRGB :: Colour -> [Double]+toListRGB (Colour r g b _) = [r, g, b]++luminance :: Colour -> Double+luminance (Colour !r !g !b _) = r * 0.3 + g * 0.6 + b * 0.1++logLuminance :: Colour -> Double+logLuminance = log . max 1e-5 . luminance
+ app/src/CornellBox.hs view
@@ -0,0 +1,162 @@+-- Cornell box reference data+{-# LANGUAGE MagicHash #-}++module CornellBox(cornellBox, cornellBoxCamera, cornellBoxLights) where++import Vector+import Primitive+import Camera+import Material+import Colour+import Shader+import Matrix+import Light+--import GHC.Types++cornellBoxLights :: [Light]++cornellBoxCamera :: Camera++cameraPosition :: Vector++floorObject :: Object+leftWallObject :: Object+rightWallObject :: Object+frontWallObject :: Object+ceilingObject :: Object+backWallObject :: Object+tallBlockObject :: Object+shortBlockObject :: Object++leftWallVertices :: [Vector]+rightWallVertices :: [Vector]+backWallVertices :: [Vector]+frontWallVertices :: [Vector]+ceilingVertices :: [Vector]+--lightVertices :: [Vector]+floorVertices :: [Vector]+tallBlockVertices :: [Vector]+shortBlockVertices :: [Vector]++whiteMaterial :: Material+redMaterial :: Material+greenMaterial :: Material++cornellBoxLights = [ + QuadLight (CommonLightData (Colour 500 500 500 0) True) (Vector 213.0 548.0 227.0 1.0) 600 (Vector 130.0 0.0 0.0 0.0) (Vector 0.0 0.0 105.0 0.0)+ ]++cameraPosition = Vector 278.0 273.0 (-800.0) 1.0+cornellBoxCamera = withVectors cameraPosition xaxis yaxis zaxis 45.0++whiteMaterial = Material (Colour 0.5 0.5 0.5 1) (Colour 0.5 0.5 0.5 1) colBlack 0 0 0 iorAir NullShader+redMaterial = Material (Colour 0.5 0.0 0.0 1) (Colour 0.5 0.0 0.0 1) colBlack 0 0 0 iorAir NullShader+greenMaterial = Material (Colour 0.0 0.5 0.0 1) (Colour 0.0 0.5 0.0 1) colBlack 0 0 0 iorAir NullShader++floorVertices = [+ Vector 556.0 0.0 0.0 1.0,+ Vector 0.0 0.0 0.0 1.0,+ Vector 0.0 0.0 559.2 1.0,+ Vector 556.0 0.0 559.2 1.0+ ]++ceilingVertices = [+ Vector 556.0 548.8 0.0 1.0,+ Vector 556.0 548.8 559.2 1.0,+ Vector 0.0 548.8 559.2 1.0,+ Vector 0.0 548.8 0.0 1.0+ ]++backWallVertices = [+ Vector 556.0 0.0 559.2 1.0,+ Vector 0.0 0.0 559.2 1.0,+ Vector 0.0 548.8 559.2 1.0,+ Vector 556.0 548.8 559.2 1.0+ ]++frontWallVertices = [+ Vector 556.0 548.8 0.0 1.0,+ Vector 0.0 548.8 0.0 1.0,+ Vector 0.0 0.0 0.0 1.0,+ Vector 556.0 0.0 0.0 1.0+ ]++rightWallVertices = [+ Vector 0.0 0.0 559.2 1.0,+ Vector 0.0 0.0 0.0 1.0,+ Vector 0.0 548.8 0.0 1.0,+ Vector 0.0 548.8 559.2 1.0+ ]++leftWallVertices = [+ Vector 556.0 0.0 0.0 1.0,+ Vector 556.0 0.0 559.2 1.0,+ Vector 556.0 548.8 559.2 1.0,+ Vector 556.0 548.8 0.0 1.0+ ]++shortBlockVertices = [+ Vector 130.0 165.0 65.0 1.0,+ Vector 82.0 165.0 225.0 1.0,+ Vector 240.0 165.0 272.0 1.0,+ Vector 290.0 165.0 114.0 1.0,++ Vector 290.0 0.0 114.0 1.0,+ Vector 290.0 165.0 114.0 1.0,+ Vector 240.0 165.0 272.0 1.0,+ Vector 240.0 0.0 272.0 1.0,++ Vector 130.0 0.0 65.0 1.0,+ Vector 130.0 165.0 65.0 1.0,+ Vector 290.0 165.0 114.0 1.0,+ Vector 290.0 0.0 114.0 1.0,++ Vector 82.0 0.0 225.0 1.0,+ Vector 82.0 165.0 225.0 1.0,+ Vector 130.0 165.0 65.0 1.0,+ Vector 130.0 0.0 65.0 1.0,++ Vector 240.0 0.0 272.0 1.0,+ Vector 240.0 165.0 272.0 1.0,+ Vector 82.0 165.0 225.0 1.0,+ Vector 82.0 0.0 225.0 1.0+ ]++tallBlockVertices = [+ Vector 423.0 330.0 247.0 1.0,+ Vector 265.0 330.0 296.0 1.0,+ Vector 314.0 330.0 456.0 1.0,+ Vector 472.0 330.0 406.0 1.0,++ Vector 423.0 0.0 247.0 1.0,+ Vector 423.0 330.0 247.0 1.0,+ Vector 472.0 330.0 406.0 1.0,+ Vector 472.0 0.0 406.0 1.0,++ Vector 472.0 0.0 406.0 1.0,+ Vector 472.0 330.0 406.0 1.0,+ Vector 314.0 330.0 456.0 1.0,+ Vector 314.0 0.0 456.0 1.0,++ Vector 314.0 0.0 456.0 1.0,+ Vector 314.0 330.0 456.0 1.0,+ Vector 265.0 330.0 296.0 1.0,+ Vector 265.0 0.0 296.0 1.0,++ Vector 265.0 0.0 296.0 1.0,+ Vector 265.0 330.0 296.0 1.0,+ Vector 423.0 330.0 247.0 1.0,+ Vector 423.0 0.0 247.0 1.0+ ]++floorObject = Object (TriangleMesh (quadsToTriangles floorVertices)) whiteMaterial identity+frontWallObject = Object (TriangleMesh (quadsToTriangles frontWallVertices)) whiteMaterial identity+leftWallObject = Object (TriangleMesh (quadsToTriangles leftWallVertices)) redMaterial identity+rightWallObject = Object (TriangleMesh (quadsToTriangles rightWallVertices)) greenMaterial identity+ceilingObject = Object (TriangleMesh (quadsToTriangles ceilingVertices)) whiteMaterial identity+backWallObject = Object (TriangleMesh (quadsToTriangles backWallVertices)) whiteMaterial identity+shortBlockObject = Object (TriangleMesh (quadsToTriangles shortBlockVertices)) whiteMaterial identity+tallBlockObject = Object (TriangleMesh (quadsToTriangles tallBlockVertices)) whiteMaterial identity++cornellBox :: [Object]+cornellBox = [ceilingObject, floorObject, leftWallObject, rightWallObject, backWallObject, frontWallObject, tallBlockObject, shortBlockObject]
+ app/src/Distribution.hs view
@@ -0,0 +1,55 @@+-- Module for generating sample patterns for distributed ray tracing+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module Distribution (generatePointsOnSphere, generatePointsOnQuad, generatePointsOnHemisphere) where++import Vector+import System.Random.Mersenne.Pure64+import Control.Monad.State++type GeneratorState = State PureMT++-- Generate a pair of random normalised floats+randomUV :: GeneratorState (Double, Double)+randomUV = do generator <- get+ let !(u, generator') = randomDouble generator+ let !(v, generator'') = randomDouble generator'+ put generator''+ return (u, v)++-- Generate a list of N random UVs+generateRandomUVs :: Int -> GeneratorState [(Double, Double)]+generateRandomUVs n = replicateM n randomUV++-- Generate a list of random points on a sphere+generatePointsOnSphere :: Int -> Double -> Int -> [Position]+generatePointsOnSphere numPoints r seed = map uvToPosition randomUVs+ where+ randomUVs = evalState (generateRandomUVs numPoints) (pureMT (fromIntegral seed))+ uvToPosition (!u, !v) = Vector (r * x) (r * y) (r * z) 1+ where+ !z = 2 * u - 1+ !t = 2 * pi * v+ !w = sqrt (1 - z * z)+ !x = w * cos t+ !y = w * sin t++-- Generate a list of random points on a hemisphere (z > 0)+generatePointsOnHemisphere :: Int -> Double -> Int -> [Position]+generatePointsOnHemisphere numPoints r seed = map uvToPosition randomUVs+ where+ randomUVs = evalState (generateRandomUVs numPoints) (pureMT (fromIntegral seed))+ uvToPosition (!u, !v) = Vector (r * x) (r * y) (r * z) 1+ where+ !y = 2 * u - 1+ !t = pi * v+ !w = sqrt (1 - y * y)+ !x = w * cos t+ !z = w * sin t++generatePointsOnQuad :: Position -> Direction -> Direction -> Int -> Int -> [Position]+generatePointsOnQuad pos deltaU deltaV numPoints seed = map uvToPosition randomUVs+ where+ randomUVs = evalState (generateRandomUVs numPoints) (pureMT (fromIntegral seed))+ uvToPosition (u, v) = pos + (deltaU <*> u) + (deltaV <*> v)
+ app/src/IrradianceCache.hs view
@@ -0,0 +1,75 @@+-- The irradiance cache+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module IrradianceCache (IrradianceCache, query, initialiseCache) where++import Vector+import Colour+import BoundingBox+import Octree+import SceneGraph++-- Irradiance gradient using a central-differencing approach+data IrradianceGradient = CentralDifferenceGradient {-# UNPACK #-} !(Colour, Colour, Colour)++-- Direction of normal, colour, radius+data CacheSample = CacheSample {-# UNPACK #-} !(Normal, Colour, Double)++type IrradianceCache = OctTree CacheSample++-- Pretty printer for cache samples+instance Show CacheSample where+ show (CacheSample (dir, col, r)) = "\tDirection: " ++ show dir ++ "\n\tColour: " ++ show col ++ "\n\tRadius: " ++ show r ++ "\n"++-- This gives an initial empty cache that will later be populated+initialiseCache :: SceneGraph -> IrradianceCache+initialiseCache sceneGraph = OctTreeNode slightlyEnlargedBox $ map OctTreeDummy (generateOctreeBoxList slightlyEnlargedBox)+ where+ -- Create the initial irradiance cache tree. This is a box a little larger than the world so that we fit any points offset along the normal etc+ slightlyEnlargedBox = growBoundingBox (finiteBox sceneGraph) 10++-- Quantify the error if we use a given sample to shade a point+-- The bigger the number, the better the estimate+errorWeight :: (Position, Direction) -> (Position, CacheSample) -> Double+{-# SPECIALIZE INLINE errorWeight :: (Position, Direction) -> (Position, CacheSample) -> Double #-}+errorWeight (pos', dir') (pos, CacheSample (dir, _, r)) = 1 / ((pos `distance` pos') / r + sqrt (1 + (dir `sdot3` dir')))++-- This slightly convoluted version is written to be tail recursive. I effectively have to maintain a software stack of the+-- nodes remaining to be traversed+findSamples :: (Position, Direction) -> [IrradianceCache] -> [(Vector, CacheSample, Double)] -> [(Vector, CacheSample, Double)]+findSamples posDir@(!pos, _) (OctTreeNode !box nodeChildren : xs) !acc+ | box `contains` pos = findSamples posDir (nodeChildren ++ xs) acc+ | otherwise = findSamples posDir xs acc+findSamples posDir@(pos, _) (OctTreeLeaf _ (samplePos, sample) : xs) !acc+ | (pos `distanceSq` samplePos) <= sampleR * sampleR && weight > minimumWeight = findSamples posDir xs ((samplePos, sample, weight) : acc)+ | otherwise = findSamples posDir xs acc+ where+ !weight = errorWeight posDir (samplePos, sample)+ (CacheSample (_, _, !sampleR)) = sample+ minimumWeight = 0.4 -- This is approximately the lower bound of the weight at the radius of the sample+findSamples posDir (OctTreeDummy _ : xs) !acc = findSamples posDir xs acc+findSamples _ [] !acc = acc++-- Sum together a list of samples and error weights+sumSamples :: [(Vector, CacheSample, Double)] -> Colour+sumSamples !samples = colourSum Colour.</> weightSum+ where+ sumSamples' !(!colAcc, !weightAcc) ((_, CacheSample (_, !col, _), !weight):xs) = sumSamples' (colAcc + col Colour.<*> weight, weightAcc + weight) xs+ sumSamples' !(!colAcc, !weightAcc) [] = (colAcc, weightAcc)+ !(!colourSum, !weightSum) = sumSamples' (colBlack, 0) samples++-- Query the irradiance given a point+-- Supplied function supplies the irradiance colour at a surface location along with the radius it is valid for+query :: IrradianceCache -> SurfaceLocation -> (SurfaceLocation -> (Colour, Double)) -> (Colour, IrradianceCache)+query irrCache !posTanSpace f = case findSamples (position, normal) [irrCache] [] of+ -- Insert a new cache sample+ [] -> let (!colour, !r) = f posTanSpace + !sample = CacheSample (normal, colour, r)+ in (colour, Octree.insert (fst posTanSpace) sample irrCache)+ -- Re-use existing cache samples+ list -> (sumSamples list, irrCache)+ where+ !position = fst posTanSpace+ !tanSpace = snd posTanSpace+ !normal = tsNormal tanSpace
+ app/src/KDTree.hs view
@@ -0,0 +1,58 @@+-- This is a module for constructing bounding volume hierarchies using a kdtree+{-# LANGUAGE MagicHash #-}++module KDTree(generateSceneGraphUsingKDTree, makeSplittingPlane, degenerateSplitList, findSplittingPlane) where+import Vector+import Primitive+import Data.List+import BoundingBox++-- This stuff is object specific ++-- What side of a plane is an object on?+onPositiveSide :: (Vector, Double) -> Object -> Bool+onPositiveSide (planeNormal, planeDist) obj = planeDist + (planeNormal `dot3` objBoxCentre) > 0.01+ where+ Just (boxMin, boxMax) = primitiveBoundingBox (primitive obj) obj+ objBoxCentre = (boxMin + boxMax) <*> 0.5++-- This stuff is generic++-- Generate a plane to split the objects along+makeSplittingPlane :: AABB -> Int -> (Vector, Double)+makeSplittingPlane (boxMin, boxMax) buildCycle = case nthLargestAxis (boxMax - boxMin) buildCycle of+ 0 -> (xaxis, -(vecX midPoint))+ 1 -> (yaxis, -(vecY midPoint))+ 2 -> (zaxis, -(vecZ midPoint))+ _ -> error "Undefined value"+ where+ midPoint = (boxMin + boxMax) <*> 0.5++-- Find a working splitting plane+findSplittingPlane :: AABB -> Int -> [t] -> ((Vector, Double) -> t -> Bool) -> Maybe (Vector, Double) +findSplittingPlane box buildCycle objs partitionFunc+ | buildCycle > 2 = Nothing+ | otherwise = if length leftObjects > 0 && length rightObjects > 0 + then Just candidateSplittingPlane+ else findSplittingPlane box (buildCycle + 1) objs partitionFunc+ where+ candidateSplittingPlane = makeSplittingPlane box buildCycle+ (leftObjects, rightObjects) = partition (partitionFunc candidateSplittingPlane) objs++-- We use this dysfunctional strategy where all our smarter ideas run out+degenerateSplitList :: (Eq t) => [t] -> ([t], [t])+degenerateSplitList objs = ([x | x <- objs, case x `elemIndex` objs of+ Just index -> odd index+ Nothing -> False], + [x | x <- objs, case x `elemIndex` objs of+ Just index -> even index+ Nothing -> False])++-- Make children using a kd tree+generateSceneGraphUsingKDTree :: [Object] -> [[Object]]+generateSceneGraphUsingKDTree objs = [leftObjects, rightObjects]+ where+ objBox = objectListBoundingBox objs+ (leftObjects, rightObjects) = case findSplittingPlane objBox 0 objs onPositiveSide of+ Nothing -> degenerateSplitList objs+ Just splittingPlane -> partition (onPositiveSide splittingPlane) objs
+ app/src/Light.hs view
@@ -0,0 +1,110 @@+-- Module for lights+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module Light (applyLight, + surfaceEpsilon, + Light(PointLight, AmbientLight, QuadLight),+ CommonLightData(CommonLightData), + LightingResult,+ position, + colour, + range, + deltaU, + deltaV, + addToPhotonMap, + common) where++import Vector+import Colour+import Ray+import Material+import Shader+import SceneGraph+import Misc+import RayTrace (findAnyIntersection)++data CommonLightData = CommonLightData { colour :: !Colour,+ addToPhotonMap :: !Bool } deriving (Show)++data Light = PointLight { common :: CommonLightData, position :: !Position, range :: !Double }+ | AmbientLight { common :: CommonLightData }+ | QuadLight { common :: CommonLightData, position :: !Position, range :: !Double, deltaU :: !Direction, deltaV :: !Direction } deriving (Show)++type LightingResult = (Colour, Colour, Colour) -- Ambient, diffuse, specular++-- Value for the surface epsilon+surfaceEpsilon :: Double+surfaceEpsilon = 0.1++-- Find the attenuation for a light source+lightAttenuation :: Vector -> Vector -> Double -> Double+lightAttenuation !lightPos !shadePos !lightRange =+ let dist = lightPos `Vector.distance` shadePos+ in if dist < lightRange then 1 - dist / lightRange else 0++-- Apply phong lighting to an object+phongLighting :: SurfaceLocation -> Light -> Material -> SceneGraph -> Direction -> Colour+phongLighting (!shadePos, !tanSpace) (PointLight (CommonLightData !lightColour !inPhotonMap') !lightPos !lightRange) objMaterial sceneGraph !viewDirection + | (lightPos `distanceSq` shadePos) < (lightRange * lightRange) && dotProd > 0 = case findAnyIntersection sceneGraph (rayWithPoints intersectionPlusEpsilon lightPos) of+ Just _ -> colBlack -- An object is closer to our point of consideration than the light, so occluded+ Nothing -> (lightColour * lightingSum) Colour.<*> attenuation+ where+ !lightingSum = diffuseLighting + specularLighting+ !attenuation = lightAttenuation lightPos shadePos lightRange+ !specularCorrection = (specularPower objMaterial + 2) / (2 * pi)+ !specularLighting = specular objMaterial Colour.<*> (specularCorrection * saturate (reflection `dot3` Vector.negate viewDirection) ** specularPower objMaterial)+ !reflection = reflect incoming normal+ !diffuseLighting = if inPhotonMap'+ then colBlack+ else shaderDiffuse * diffuse objMaterial Colour.<*> saturate dotProd+ !shaderDiffuse = evaluateDiffuse (shader objMaterial) shadePos tanSpace+ | otherwise = colBlack+ where + !intersectionPlusEpsilon = shadePos + (normal Vector.<*> surfaceEpsilon)+ !incoming = normalise (lightPos - shadePos)+ !dotProd = normal `dot3` incoming+ !normal = thr tanSpace+phongLighting _ (AmbientLight (CommonLightData _ _)) _ _ _ = error "phongLighting: Do not know how to handle AmbientLight"+phongLighting (!shadePos, !tanSpace) (QuadLight (CommonLightData !lightColour !inPhotonMap') !lightPos !lightRange !du !dv) objMaterial sceneGraph !viewDirection + | (lightCentre `distanceSq` shadePos) < (lightRange * lightRange) && dotProd > 0 = case findAnyIntersection sceneGraph (rayWithPoints intersectionPlusEpsilon lightCentre) of+ Just _ -> colBlack -- An object is closer to our point of consideration than the light, so occluded+ Nothing -> (lightColour * lightingSum) Colour.<*> attenuation+ where+ !lightingSum = diffuseLighting + specularLighting+ !attenuation = lightAttenuation lightCentre shadePos lightRange+ !specularCorrection = (specularPower objMaterial + 2) / (2 * pi)+ !specularLighting = specular objMaterial Colour.<*> (specularCorrection * saturate (reflection `dot3` Vector.negate viewDirection) ** specularPower objMaterial)+ !reflection = reflect incoming normal+ !diffuseLighting = if inPhotonMap'+ then colBlack+ else shaderDiffuse * diffuse objMaterial Colour.<*> saturate dotProd+ !shaderDiffuse = evaluateDiffuse (shader objMaterial) shadePos tanSpace+ | otherwise = colBlack+ where + !lightCentre = lightPos + du Vector.<*> 0.5 + dv Vector.<*> 0.5+ !intersectionPlusEpsilon = shadePos + (normal Vector.<*> surfaceEpsilon)+ !incoming = normalise (lightCentre - shadePos)+ !dotProd = normal `dot3` incoming+ !normal = thr tanSpace++-- For a given surface point, work out the lighting, including occlusion+applyLight :: SceneGraph -> SurfaceLocation -> Material -> Direction -> Light -> Colour+applyLight sceneGraph !intersectionPointNormal !objMaterial !viewDirection light@(PointLight (CommonLightData _ _) _ _)+ = phongLighting + intersectionPointNormal + light+ objMaterial + sceneGraph+ viewDirection+applyLight _ (!intersectionPoint, !intersectionTanSpace) !objMaterial _ (AmbientLight (CommonLightData !ambientColour _)) = + let shaderAmbient = evaluateAmbient (shader objMaterial) intersectionPoint intersectionTanSpace+ materialAmbient = ambient objMaterial+ in ambientColour * shaderAmbient * materialAmbient+applyLight sceneGraph !intersectionPointNormal !objMaterial !viewDirection light@(QuadLight (CommonLightData _ _) _ _ _ _)+ = phongLighting + intersectionPointNormal + light+ objMaterial + sceneGraph+ viewDirection
+ app/src/Light.hs-boot view
@@ -0,0 +1,29 @@+module Light (applyLight, + surfaceEpsilon, + Light(PointLight, AmbientLight, QuadLight), + CommonLightData(CommonLightData), + LightingResult, + position, + colour, + range, + deltaU, + deltaV, + addToPhotonMap, + common) where++import Vector+import Colour+import Material+import SceneGraph++data CommonLightData = CommonLightData { colour :: !Colour,+ addToPhotonMap :: !Bool } ++data Light = PointLight { common :: CommonLightData, position :: !Position, range :: !Double }+ | AmbientLight { common :: CommonLightData }+ | QuadLight { common :: CommonLightData, position :: !Position, range :: !Double, deltaU :: !Direction, deltaV :: !Direction } ++type LightingResult = (Colour, Colour, Colour) -- Ambient, diffuse, specular++applyLight :: SceneGraph -> (Position, TangentSpace) -> Material -> Direction -> Light -> Colour+surfaceEpsilon :: Double
+ app/src/Main.hs view
@@ -0,0 +1,148 @@+-- Main module of raytracer++import Data.Bits+import Data.ByteString hiding (map)+import System.Console.GetOpt+import System.Environment+import RayTrace+import Colour+import SceneGraph+import KDTree+import CornellBox+import GHC.Conc (numCapabilities)+import Codec.BMP+import PhotonMap+import RenderContext+import Light+import ToneMap+import Control.Arrow++-- Command line option support+data Option+ = ShowIntermediate -- -i+ | PhotonMap -- -p+ | DirectPhotonMapVisualisation -- -v+ | DistributedRayTracing -- d+ | IrradianceCaching -- c+ deriving (Eq, Ord, Enum, Show, Bounded)++options :: [OptDescr Option]+options = [+ Option ['i'] [] (NoArg ShowIntermediate) "Show intermediates",+ Option ['p'] [] (NoArg PhotonMap) "Photon map",+ Option ['v'] [] (NoArg DirectPhotonMapVisualisation) "Direct photon map visualisation",+ Option ['d'] [] (NoArg DistributedRayTracing) "Distributed ray tracing",+ Option ['c'] [] (NoArg IrradianceCaching) "Irradiance caching"+ ]++parsedOptions :: [String] -> [Option]+parsedOptions argv = case getOpt Permute options argv of+ (args,_,[]) -> args+ (_,_,_) -> []++-- Some hardcoded values, at present+renderWidth :: Int -> Int+renderWidth mipLevel = 1280 `shiftR` mipLevel++renderHeight :: Int -> Int+renderHeight mipLevel = 720 `shiftR` mipLevel++-- This returns a list of colours of pixels+renderImage :: Int -> RenderContext -> Maybe PhotonMap -> [Colour]+renderImage mipLevel renderSettings photonMap = finalImage+ where+ rawImageOutput = rayTraceImage renderSettings cornellBoxCamera (renderWidth mipLevel) (renderHeight mipLevel) photonMap+ exposedImage = exposeImage imageAverageLogLuminance rawImageOutput 4+ toneMappedImage = toneMapImage toneMapHejlBurgessDawson exposedImage+ finalImage = map (clamp . invGammaCorrect) toneMappedImage++-- In the interest of rapid developer feedback, this functions writes a progressively-increasing image+-- So, we get quick feedback on the intermediate results, but will still ultimately get the final image+-- Note this does no re-use, so it'll be slower overall+writeRaytracedImage :: [Int] -> Maybe PhotonMap -> RenderContext -> IO ()+writeRaytracedImage [] photonMap renderSettings = do+ let imageData = renderImage 0 renderSettings photonMap+ let rgba = Data.ByteString.pack (convertColoursToPixels imageData)+ let bmp = packRGBA32ToBMP (renderWidth 0) (renderHeight 0) rgba+ Prelude.putStrLn "Performing final render"+ writeBMP "test.bmp" bmp+writeRaytracedImage (mipLevel:mipLevels) photonMap renderSettings = do+ let imageData = renderImage mipLevel renderSettings photonMap+ let rgba = Data.ByteString.pack (convertColoursToPixels imageData)+ let bmp = packRGBA32ToBMP (renderWidth mipLevel) (renderHeight mipLevel) rgba+ let filename = "test-intermediate-" ++ show mipLevel ++ ".bmp"+ Prelude.putStrLn filename+ writeBMP filename bmp+ writeRaytracedImage mipLevels photonMap renderSettings++-- Strip off the photon map flag from a light+notInPhotonMap :: Light -> Light+notInPhotonMap (PointLight (CommonLightData colour' _) position' range') = PointLight (CommonLightData colour' False) position' range'+notInPhotonMap (AmbientLight (CommonLightData colour' _)) = AmbientLight (CommonLightData colour' False)+notInPhotonMap (QuadLight (CommonLightData colour' _) position' range' deltaU' deltaV') = QuadLight (CommonLightData colour' False) position' range' deltaU' deltaV'++-- Main function+main :: IO ()+main = do + args <- getArgs+ let opts = parsedOptions args++ let renderSettings = RenderContext + numDistributedSamples + (buildSceneGraph cornellBox generateSceneGraphUsingKDTree) + cornellBoxLights + maxRayDepth + reflectionDistance + refractionDistance + (PhotonMapContext photonGatherDistance maxGatherPhotons coneFilterConstant directPhotonMapVisualisation) + rayOriginDistribution'+ depthOfFieldFocalDistance'+ renderMode'+ enableIrradianceCache+ where+ -- Ray trace constants+ numDistributedSamples = if DistributedRayTracing `Prelude.elem` opts + then 64 + else 1+ maxRayDepth = 5+ reflectionDistance = 1000+ refractionDistance = 1000+ -- Photon constants+ photonGatherDistance = 100+ maxGatherPhotons = 200+ coneFilterConstant = 2+ -- Depth of field constants+ rayOriginDistribution' = 0.5+ depthOfFieldFocalDistance' = 400+ renderMode'+ | PhotonMap `Prelude.elem` opts = PhotonMapper+ | otherwise = RayTrace+ directPhotonMapVisualisation = DirectPhotonMapVisualisation `Prelude.elem` opts+ enableIrradianceCache = IrradianceCaching `Prelude.elem` opts++ -- Display hardware capabilities+ Prelude.putStrLn $ "Running on " ++ show numCapabilities ++ " cores"++ -- Create a photon map, if necessary+ let doPhotonMapping = PhotonMap `Prelude.elem` opts+ let photonMapMessage = if doPhotonMapping + then if DirectPhotonMapVisualisation `Prelude.elem` opts+ then "Directly visualising photon map"+ else "Creating photon map..." + else "Photon mapping disabled"+ Prelude.putStrLn photonMapMessage+ let thousand = 1000+ let numPhotons = 200 * thousand+ let (photonMap, lights')+ | doPhotonMapping = Control.Arrow.first Just $ + buildPhotonMap (sceneGraph renderSettings) cornellBoxLights numPhotons+ | otherwise = (Nothing, map notInPhotonMap (lights renderSettings))++ -- Render the image+ let renderSettings' = renderSettings { lights = lights' }+ let maxMipLevel = 8+ let intermediateMipLevels = if ShowIntermediate `Prelude.elem` opts+ then Prelude.reverse [1..maxMipLevel]+ else []+ Prelude.putStrLn "Rendering image..."+ writeRaytracedImage intermediateMipLevels photonMap renderSettings'
+ app/src/Material.hs view
@@ -0,0 +1,21 @@+-- Materials of an object++module Material where++import Colour+import Shader++data Material = Material { ambient :: {-# UNPACK #-} !Colour, + diffuse :: {-# UNPACK #-} !Colour, + specular :: {-# UNPACK #-} !Colour, + specularPower :: {-# UNPACK #-} !Double,+ reflectivity :: {-# UNPACK #-} !Double,+ transmit :: {-# UNPACK #-} !Double,+ indexOfRefraction :: {-# UNPACK #-} !Double,+ shader :: Shader } deriving (Show, Eq)++iorAir :: Double+iorAir = 1.000293++iorWater :: Double+iorWater = 1.3330
+ app/src/Matrix.hs view
@@ -0,0 +1,55 @@+-- 4D Matrix Library+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module Matrix where++import Vector++data Matrix = Matrix ![Double] deriving (Show, Read, Eq)++-- Just pass back the identity matrix+identity :: Matrix+identity = Matrix [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1]++-- This multiplies one element of the matrix against+mulElement :: Matrix -> Matrix -> Int -> Int -> Double+mulElement (Matrix mat1) (Matrix mat2) i j = sum [(mat1 !! (i * 4 + k)) * (mat2 !! (k * 4 + j)) | k <- [0..3]]++-- Multiply together two matrices+-- TODO - Turn mulElement into a lambda function?+mul :: Matrix -> Matrix -> Matrix+mul a b = Matrix [mulElement a b i j | j <- [0..3], i <- [0..3]]++-- Need matrix inversion code++-- Vector * Matrix+transformVector :: Matrix -> Vector -> Vector+transformVector (Matrix mat) (Vector !x !y !z !w) = Vector x' y' z' w'+ where+ !vec = [x, y, z, w]+ !xvector = take 4 mat+ !yvector = take 4 (drop 4 mat)+ !zvector = take 4 (drop 8 mat)+ !wvector = take 4 (drop 12 mat)+ !x' = sum $ zipWith (*) vec xvector+ !y' = sum $ zipWith (*) vec yvector+ !z' = sum $ zipWith (*) vec zvector+ !w' = sum $ zipWith (*) vec wvector++-- Build a matrix from 4 vectors+buildMatrix :: Vector -> Vector -> Vector -> Vector -> Matrix+buildMatrix (Vector xx xy xz _) (Vector yx yy yz _) (Vector zx zy zz _) (Vector px py pz _) = Matrix ([xx, xy, xz] ++ [px] +++ [yx, yy, yz] ++ [py] +++ [zx, zy, zz] ++ [pz] +++ [0, 0, 0, 1])+getTranslation :: Matrix -> Vector+{-# SPECIALIZE INLINE getTranslation :: Matrix -> Vector #-}+getTranslation (Matrix a) = Vector x y z 1+ where+ [!x, !y, !z, _] = drop 12 a++translationMatrix :: Double -> Double -> Double -> Matrix+translationMatrix !x !y !z = Matrix [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, x, y, z, 1]+translationMatrix' :: Vector -> Matrix+translationMatrix' (Vector !x !y !z _) = Matrix [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, x, y, z, 1]
+ app/src/Misc.hs view
@@ -0,0 +1,31 @@+-- Various assorted bits and pieces+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE BangPatterns #-}++module Misc where+import GHC.Prim+import GHC.Types+import Data.List++degreesToRadians :: Double -> Double+degreesToRadians x = x * pi / 180++xor :: Bool -> Bool -> Bool+xor True a = not a+xor False a = a++thr :: (x, y, z) -> z+thr (_, _, c) = c++-- Little helper for saturation+saturate :: (Num t, Ord t) => t -> t+saturate x = Prelude.max 0 (Prelude.min x 1)++saturate## :: Double# -> Double#+saturate## !x = value+ where+ !(D# !value) = Prelude.max 0 (Prelude.min (D# x) 1)++harmonicMean :: (Num t, Fractional t) => [t] -> t+harmonicMean (x:xs) = fromIntegral (length (x:xs)) / foldl' (\a b -> b + 1 / a) 0 (x:xs)+harmonicMean [] = 0
+ app/src/Octree.hs view
@@ -0,0 +1,133 @@+-- This is a module for constructing bounding volume hierarchies using an octree approach+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module Octree(generateSceneGraphUsingOctree, generateOctreeBoxList, OctTree(OctTreeNode, OctTreeLeaf, OctTreeDummy), create, Octree.insert, gather) where++import Vector+import Primitive+import BoundingBox+import Control.Parallel.Strategies++data OctTree a = OctTreeDummy !AABB+ | OctTreeNode !AABB [OctTree a]+ | OctTreeLeaf !AABB !(Vector, a) deriving (Eq)++instance Show a => Show (OctTree a) where+ show = display 0++tabs :: String+tabs = '\t' : tabs++display :: (Show a) => Int -> OctTree a -> String+display level (OctTreeDummy box) = take level tabs ++ "[Dummy] box=" ++ show box ++ "\n"+display level (OctTreeNode box children) = take level tabs ++ "[Node] box=" ++ show box ++ "\n" ++ concatMap (display (level + 1)) children ++ "\n"+display level (OctTreeLeaf box (pos, value)) = take level tabs ++ "[Leaf] box=" ++ show box ++ " pos=" ++ show pos ++ " value=" ++ show value ++ "\n"++create :: AABB -> OctTree a+create box = OctTreeNode box $ map OctTreeDummy (generateOctreeBoxList box)++-- This performs a map, and passes through the state of the completed operation to the next recursion+-- Couldn't work out the equivalent using the state monad etc+mapS :: (a -> s -> (b, s)) -> [a] -> s -> ([b], s)+mapS f xs state = mapS' f xs state []++mapS' :: (a -> s -> (b, s)) -> [a] -> s -> [b] -> ([b], s)+mapS' !f !(x:xs) !state !acc = seq (result, state') $ mapS' f xs state' (result : acc)+ where (!result, !state') = f x state `using` rseq+mapS' _ [] !state !acc = (acc, state)++-- Insert into an octree+insert :: Vector -> a -> OctTree a -> OctTree a+insert pos a oct = fst $ insert' pos oct (Just a)++insert' :: Vector -> OctTree a -> Maybe a -> (OctTree a, Maybe a)+insert' pos oct@(OctTreeDummy box) state = case state of+ -- If we have been passed some state then attempt to consume it+ Just !value -> if box `contains` pos+ then (OctTreeLeaf box (pos, value), Nothing)+ else (oct, state)+ _ -> (oct, state)++insert' pos oct@(OctTreeNode box nodeChildren) !state = if box `contains` pos+ then let (nodeChildren', state') = mapS (insert' pos) nodeChildren state + in (OctTreeNode box nodeChildren', state')+ else (oct, state)++insert' pos oct@(OctTreeLeaf box (pos', a')) state = if box `contains` pos + then + -- First up, we turn this leaf into a node with 8 children+ -- Discard result of mapS - we assume it returns Nothing+ -- Then, re-insert the original value into our nascent octree+ let (!newChildren, _) = mapS (insert' pos) (map OctTreeDummy (generateOctreeBoxList box)) state+ (!octTree', !state') = insert' pos' (OctTreeNode box newChildren) (Just a')+ in (octTree', state')+ else (oct, state)++-- Gather data within a sphere from an octree+gather :: Position -> Double -> OctTree a -> [(a, Double)]+gather pos r (OctTreeNode box nodeChildren) = if overlapsSphere box pos r+ then concatMap (gather pos r) nodeChildren+ else []+gather pos r (OctTreeLeaf _ (pos', a))+ | dSq <= r * r = [(a, dSq)]+ | otherwise = []+ where !dSq = pos `distanceSq` pos'+gather _ _ (OctTreeDummy _) = []++-- Generate a scene graph using an octree. Refactor this to just be an octree later+generateOctreeBoxList :: AABB -> [AABB]+generateOctreeBoxList (Vector xmin ymin zmin _, Vector xmax ymax zmax _) =+ [+ (Vector xmin ymin zmin 1, Vector centreX centreY centreZ 1),+ (Vector centreX ymin zmin 1, Vector xmax centreY centreZ 1),+ (Vector xmin centreY zmin 1, Vector centreX ymax centreZ 1),+ (Vector centreX centreY zmin 1, Vector xmax ymax centreZ 1),++ (Vector xmin ymin centreZ 1, Vector centreX centreY zmax 1),+ (Vector centreX ymin centreZ 1, Vector xmax centreY zmax 1),+ (Vector xmin centreY centreZ 1, Vector centreX ymax zmax 1),+ (Vector centreX centreY centreZ 1, Vector xmax ymax zmax 1)+ ]+ where+ !centreX = (xmin + xmax) * 0.5+ !centreY = (ymin + ymax) * 0.5+ !centreZ = (zmin + zmax) * 0.5++-- Octree code that's spilt out from other modules... this is scene graph specific helper code rather than self-contained octree stuff++-- Take a list of objects and split it into a list of objects that intersect a box, and those that don't+objectsIntersectingBox :: [Object] -> AABB -> ([Object], [Object])+objectsIntersectingBox objects box = objectsIntersectingBox' objects box ([], [])++objectsIntersectingBox' :: [Object] -> AABB -> ([Object], [Object]) -> ([Object], [Object])+objectsIntersectingBox' (obj:objs) box (currentHit, currentMiss) = if intersectsBox (primitive obj) (transform obj) box+ then objectsIntersectingBox' objs box (obj : currentHit, currentMiss) + else objectsIntersectingBox' objs box (currentHit, obj : currentMiss)+objectsIntersectingBox' [] _ (currentHit, currentMiss) = (currentHit, currentMiss)++-- Iterator function. Match up objects to this box, and then iterate with the remainder+assignObjectsToOctreeBoxes' :: [Object] -> [AABB] -> [[Object]] -> [[Object]]+assignObjectsToOctreeBoxes' objs (box:boxes) (x:xs) = assignObjectsToOctreeBoxes' remainingObjects boxes (matchedObjects : x : xs)+ where+ (matchedObjects, remainingObjects) = objectsIntersectingBox objs box+assignObjectsToOctreeBoxes' _ [] currentList = currentList+assignObjectsToOctreeBoxes' objs (box:boxes) [] = assignObjectsToOctreeBoxes' remainingObjects boxes [matchedObjects]+ where+ (matchedObjects, remainingObjects) = objectsIntersectingBox objs box++-- Generate the list of objects for each bounding box+assignObjectsToOctreeBoxes :: [Object] -> [AABB] -> [[Object]]+assignObjectsToOctreeBoxes objects boxes = assignObjectsToOctreeBoxes' objects boxes []++-- Make children using an octree algorithm+generateSceneGraphUsingOctree :: [Object] -> [[Object]]+generateSceneGraphUsingOctree (obj:objs) + | not (boundingBoxValid nodeBox) = error "Invalid bounding box"+ | otherwise = onlyPopulatedBoxes+ where+ nodeBox = objectListBoundingBox (obj:objs)+ octreeBoxes = generateOctreeBoxList nodeBox+ objsPerOctreeBox = assignObjectsToOctreeBoxes (obj:objs) octreeBoxes+ onlyPopulatedBoxes = filter (\x -> length x > 0) objsPerOctreeBox+generateSceneGraphUsingOctree [] = []
+ app/src/PhotonMap.hs view
@@ -0,0 +1,261 @@+-- Photon mapping+{-# LANGUAGE BangPatterns #-}++module PhotonMap(buildPhotonMap, PhotonMap(photonList), irradiance, PhotonMapContext(PhotonMapContext)) where++import {-# SOURCE #-} Light hiding (position)+import Vector+import Distribution+import Material+import Colour+import SceneGraph+import RayTrace+import Ray hiding (direction)+import Control.Monad.State+import BoundingBox+import KDTree+import Debug.Trace+import Misc+import Control.Parallel.Strategies+import Control.DeepSeq+import Data.Heap hiding (partition)+import System.Random.Mersenne.Pure64+import Data.List hiding (union, insert)+import Primitive++type GeneratorState = State PureMT++data PhotonMapContext = PhotonMapContext {+ photonGatherDistance :: Double,+ maxGatherPhotons :: Int,+ coneFilterK :: Double,+ directVisualisation :: Bool }++data Photon = Photon { power :: {-# UNPACK #-} !Colour, posDir :: {-# UNPACK #-} !(Position, Direction) } deriving (Show, Eq, Ord)++data PhotonMapTree = PhotonMapNode {-# UNPACK #-} !Int {-# UNPACK #-} !Double PhotonMapTree PhotonMapTree+ | PhotonMapLeaf {-# UNPACK #-} !Photon deriving (Show, Eq)++data PhotonMap = PhotonMap { photonList :: [Photon],+ photonMapTree :: PhotonMapTree } deriving(Show, Eq)++data PhotonChoice = DiffuseReflect | SpecularReflect | Absorb++instance NFData Photon where+ rnf (Photon power' posDir') = rnf power' `seq` rnf posDir'++seedToRefactor :: Int+seedToRefactor = 12345++-- Generate a list of photon position and direction tuples to emit+-- I zip up each pos,dir tuple with a random number generator to give each photon a different sequence of random values+-- Helps parallelisation...+-- TODO Eliminate magic number seeds from here+emitPhotons :: Light -> Int -> [(Position, Direction, PureMT, Colour)]+emitPhotons (PointLight (CommonLightData lightPower True) pos _) numPhotons = zipWith (\dir num -> (pos, dir, pureMT (fromIntegral num), flux)) (generatePointsOnSphere numPhotons 1 seedToRefactor) [1..numPhotons]+ where+ flux = lightPower Colour.<*> (1.0 / fromIntegral numPhotons)+emitPhotons (QuadLight (CommonLightData lightPower True) corner _ du dv) numPhotons = zipWith3 (\pos dir num -> (pos, transformDir dir tanSpace, pureMT (fromIntegral num), flux)) randomPoints randomDirs [1..numPhotons]+ where+ randomPoints = generatePointsOnQuad corner du dv numPhotons seedToRefactor+ randomDirs = generatePointsOnHemisphere numPhotons 1 (seedToRefactor * 10)+ area = Vector.magnitude (du `cross` dv)+ flux = lightPower Colour.<*> (area / fromIntegral numPhotons)+ tanSpace = (normalise du, normalise dv, normalise (du `cross` dv))+emitPhotons _ _ = []++-- Compute russian roulette coefficients+russianRouletteCoefficients :: Material -> (Double, Double)+russianRouletteCoefficients mat = (diffuseP, specularP)+ where+ (Colour diffuseR diffuseG diffuseB _) = Material.diffuse mat+ (Colour specularR specularG specularB _) = Material.specular mat+ diffuseP = (diffuseR + diffuseG + diffuseB) / 3+ specularP = (specularR + specularG + specularB) / 3++-- Decide what to do with a photon+choosePhotonFate :: (Double, Double) -> GeneratorState PhotonChoice+choosePhotonFate (diffuseP, specularP) = do+ generator <- get+ let (p, generator') = randomDouble generator+ let result | p < diffuseP = DiffuseReflect+ | p < (diffuseP + specularP) = SpecularReflect+ | otherwise = Absorb+ put generator'+ return $! result++-- Compute new power for a photon+computeNewPhotonPower :: PhotonChoice -> (Double, Double) -> Colour -> Material -> Colour+computeNewPhotonPower fate (diffuseP, specularP) photonPower mat = case fate of+ DiffuseReflect -> photonPower * diffuse mat Colour.</> diffuseP+ SpecularReflect -> photonPower * specular mat Colour.</> specularP+ Absorb -> colBlack++-- Compute a new diffuse reflection in spherical co-ordinates+generateUV :: GeneratorState (Double, Double)+generateUV = do generator <- get+ let (u, generator') = randomDouble generator+ let (v, generator'') = randomDouble generator'+ put generator''+ return $! (u, v)++-- Find a diffuse reflection direction in the hemisphere of the normal+-- Realistic Image Synthesis Using Photon Mapping - Eq 2.24+diffuseReflectionDirection :: PureMT -> TangentSpace -> (Direction, PureMT)+diffuseReflectionDirection stdGen tanSpace = (transformDir dir tanSpace, stdGen')+ where+ ((u, v), stdGen') = runState generateUV stdGen+ theta = acos (sqrt u)+ phi = 2 * pi * v+ dir = sphericalToDirection theta phi++-- Main working photon tracing function+-- Realistic Image Synthesis Using Photon Mapping p60+tracePhoton :: [Photon] -> Photon -> SceneGraph -> PureMT -> (Int, Int) -> [Photon]+tracePhoton currentPhotons (Photon photonPower photonPosDir) sceneGraph rndState (bounce, maxBounces) = + -- See if the photon intersects a surfaces+ case findNearestIntersection sceneGraph ray of+ Nothing -> currentPhotons+ Just (obj, t, subId) -> case photonFate of+ -- Diffuse reflection. Here, we store the photon that got reflected, and trace a new photon - but only if it's bright enough to be worthwhile+ DiffuseReflect -> if Colour.magnitude newPhotonPower > brightnessEpsilon && (bounce + 1) <= maxBounces+ then tracePhoton (storedPhoton : currentPhotons) reflectedPhoton sceneGraph rndState'' (bounce + 1, maxBounces)+ else storedPhoton : currentPhotons+ where+ reflectedPhoton = Photon newPhotonPower (surfacePos, reflectedDir)+ (reflectedDir, rndState'') = diffuseReflectionDirection rndState' tanSpace++ -- Specular reflection. Here, we reflect the photon in the fashion that the surface would reflect towards the viewer and+ -- aim to absorb it somewhere else in the photon map+ SpecularReflect -> if Colour.magnitude newPhotonPower > brightnessEpsilon && (bounce + 1) <= maxBounces+ then tracePhoton currentPhotons reflectedPhoton sceneGraph rndState' (bounce + 1, maxBounces)+ else currentPhotons+ where+ reflectedPhoton = Photon newPhotonPower (surfacePos, reflectedDir)+ reflectedDir = Vector.negate (snd photonPosDir) `reflect` normal++ -- Absorb. The photon simply gets absorbed into the map+ Absorb -> storedPhoton : currentPhotons+ where+ (photonFate, rndState') = runState (choosePhotonFate coefficients) rndState+ coefficients = russianRouletteCoefficients (material obj)+ newPhotonPower = computeNewPhotonPower photonFate coefficients photonPower (material obj)+ tanSpace = primitiveTangentSpace (primitive obj) subId hitPosition obj+ normal = thr tanSpace+ hitPosition = pointAlongRay ray t+ surfacePos = hitPosition + (normal Vector.<*> surfaceEpsilon)+ brightnessEpsilon = 0.1+ storedPhoton = Photon photonPower (surfacePos, snd photonPosDir)+ where+ ray = rayWithPosDir photonPosDir 10000++-- Build a list of photons for a light source+tracePhotonsForLight :: Int -> SceneGraph -> Light -> [Photon]+tracePhotonsForLight numPhotons sceneGraph light = concat (map (\(pos, dir, rndState, flux) -> tracePhoton [] (Photon flux (pos, dir)) sceneGraph rndState (0, maxBounces)) posDirGens `using` parListChunk photonsPerChunk rdeepseq)+ where+ posDirGens = emitPhotons light numPhotons -- Positions, directions, random number generators+ maxBounces = 500+ photonsPerChunk = 256++-- High-level function to build a photon map+buildPhotonMap :: SceneGraph -> [Light] -> Int -> (PhotonMap, [Light])+buildPhotonMap sceneGraph lights numPhotonsPerLight = photons `seq` kdTree `seq` (PhotonMap photons kdTree, lightsNotForPhotonMap)+ where+ (lightsForPhotonMap, lightsNotForPhotonMap) = partition (addToPhotonMap . common) lights+ photons = concatMap (tracePhotonsForLight numPhotonsPerLight sceneGraph) lightsForPhotonMap+ kdTree = buildKDTree photons++-- Make a bounding box of a list of photons+photonsBoundingBox :: [Photon] -> AABB+photonsBoundingBox = foldl' (\box photon -> enlargeBoundingBox (fst . posDir $ photon) box) initialInvalidBox++-- Construct a balanced kd tree of photons+-- Realistic Image Synthesis Using Photon Mapping p72+buildKDTree :: [Photon] -> PhotonMapTree+buildKDTree (x:[]) = PhotonMapLeaf x+buildKDTree [] = error "buildKDTree [] should never get called"+buildKDTree photons = let (boxMin, boxMax) = photonsBoundingBox photons+ axis = largestAxis (boxMax - boxMin)+ numPhotons = fromIntegral (length photons)+ photonsMedian = foldl' (\box photon -> box + (fst . posDir $ photon)) zeroVector photons Vector.</> numPhotons+ value = component photonsMedian axis+ photonsGT = Prelude.filter (\p -> component ((fst . posDir) p) axis > value) photons+ photonsLE = Prelude.filter (\p -> component ((fst . posDir) p) axis <= value) photons+ in if length photonsGT > 0 && length photonsLE > 0+ then let gtTree = buildKDTree photonsGT+ leTree = buildKDTree photonsLE+ in gtTree `seq` leTree `seq` PhotonMapNode axis value gtTree leTree+ else let (photons0', photons1') = trace "Using degenerate case" $ degenerateSplitList photons in PhotonMapNode axis value (buildKDTree photons0') (buildKDTree photons1')++-- Use a max heap to make it easy to eliminate distant photons+data GatheredPhoton = GatheredPhoton Double Photon deriving (Show)+type PhotonHeap = MaxHeap GatheredPhoton++instance Ord GatheredPhoton where+ compare (GatheredPhoton dist1 _) (GatheredPhoton dist2 _) = dist1 `compare` dist2++instance Eq GatheredPhoton where+ (GatheredPhoton dist1 _) == (GatheredPhoton dist2 _) = dist1 == dist2++instance NFData GatheredPhoton where+ rnf (GatheredPhoton dist photon) = rnf dist `seq` rnf photon++-- Return the minimum squared search radius from that specified, versus the furthest photon in the heap+-- We don't want to locate any photons further away than our current furthest - we're looking for the closest ones, after all+minimalSearchRadius :: Double -> PhotonHeap -> Double+minimalSearchRadius !rSq photonHeap = case viewHead photonHeap of+ Nothing -> rSq+ Just (GatheredPhoton !dSq _) -> Prelude.min rSq dSq++-- Gather photons for irradiance computations+-- Algorithm adapted from Realistic Image Synthesis Using Photon Mapping p73+gatherPhotons :: PhotonMapTree -> Position -> Double -> PhotonHeap -> Int -> PhotonHeap+gatherPhotons (PhotonMapNode axis value gtChild leChild) pos rSq photonHeap maxPhotons+ -- In this case, the split plane bisects the search sphere - search both halves of tree+ | (value - posComponent) ** 2 <= rSq = let heap1 = gatherPhotons gtChild pos rSq' photonHeap maxPhotons+ rSq'' = minimalSearchRadius rSq' heap1+ heap2 = gatherPhotons leChild pos rSq'' photonHeap maxPhotons+ newHeap = union heap1 heap2+ in heap1 `seq` heap2 `seq` newHeap `seq` Data.Heap.drop (size newHeap - maxPhotons) newHeap++ -- One side of the tree...+ | posComponent > value = gatherPhotons gtChild pos rSq' photonHeap maxPhotons++ -- ... or the other+ | posComponent <= value = gatherPhotons leChild pos rSq' photonHeap maxPhotons++ -- Prolapse+ | otherwise = error "gatherPhotons: unexplained/unexpected case here"+ where+ posComponent = component pos axis+ rSq' = minimalSearchRadius rSq photonHeap -- Refine search radius as we go down tree to search no further than closest allowed photon+gatherPhotons (PhotonMapLeaf p) pos rSq photonHeap maxPhotons+ | distSq < rSq = let newHeap = insert (GatheredPhoton distSq p) photonHeap+ in Data.Heap.drop (size newHeap - maxPhotons) newHeap -- Discard any excess photons - we get rid of the furthest ones+ | otherwise = photonHeap+ where distSq = pos `distanceSq` (fst . posDir) p++-- Return the contribution of a given photon, including a simple cos term to emulate BRDF plus the cone filter+-- Cone filter is from Realistic Image Synthesis Using Photon Mapping p81+photonContribution :: Double -> SurfaceLocation -> Photon -> Colour+photonContribution kr (pos, (_, _, normal)) photon = power photon Colour.<*> ((Vector.negate normal `sdot3` (snd . posDir) photon) * weight)+ where+ weight = 1 - (pos `distance` (fst . posDir) photon) / (kr + 0.000000001) -- Add on an epsilon to prevent div0 in cone filter++-- Find the overall contribution of a list of photons+-- Radiance estimate algorithm from Realistic Image Synthesis Using Photon Mapping p81+sumPhotonContribution :: Double -> Double -> SurfaceLocation -> [Photon] -> Colour+sumPhotonContribution r k posTanSpace photons = foldl' (\y x -> y + photonContribution (k * r) posTanSpace x) colBlack photons Colour.<*> (1.0 / ((1.0 - 2.0 / (3.0 * k)) * pi * r * r))++-- Look up the resulting irradiance from the photon map at a given point+-- Realistic Image Synthesis Using Photon Mapping, e7.6+irradiance :: PhotonMap -> PhotonMapContext -> Material -> SurfaceLocation -> Colour+irradiance photonMap photonMapContext mat posTanSpace = sumPhotonContribution r k posTanSpace gatheredPhotons * diffuse mat+ where+ r = photonGatherDistance photonMapContext+ maxPhotons+ | directVisualisation photonMapContext = 1+ | otherwise = maxGatherPhotons photonMapContext+ k = coneFilterK photonMapContext+ photonHeap = gatherPhotons (photonMapTree photonMap) (fst posTanSpace) (r * r) Data.Heap.empty maxPhotons+ gatheredPhotons = map (\(GatheredPhoton _ photon) -> photon) (Data.Heap.take maxPhotons photonHeap)
+ app/src/PhotonMap.hs-boot view
@@ -0,0 +1,28 @@+-- HS-boot interface file for PhotonMap to break circular dependencies++module PhotonMap(buildPhotonMap, PhotonMap(photonList), irradiance, PhotonMapContext(PhotonMapContext)) where++import Colour+import Vector+import SceneGraph+import Material+import {-# SOURCE #-} Light++data PhotonMapContext = PhotonMapContext {+ photonGatherDistance :: Double,+ maxGatherPhotons :: Int,+ coneFilterK :: Double,+ directVisualisation :: Bool }++data Photon = Photon { power :: !Colour,+ position :: !Position,+ direction :: !Direction }++data PhotonMapTree = PhotonMapNode { splitAxis :: Int, splitValue :: Double, child0 :: PhotonMapTree, child1 :: PhotonMapTree} + | PhotonMapLeaf { photon :: Photon }++data PhotonMap = PhotonMap { photonList :: [Photon],+ photonMapTree :: PhotonMapTree }++irradiance :: PhotonMap -> PhotonMapContext -> Material -> (Position, TangentSpace) -> Colour+buildPhotonMap :: SceneGraph -> [Light] -> Int -> (PhotonMap, [Light])
+ app/src/Primitive.hs view
@@ -0,0 +1,298 @@+-- Module for general primitives and intersections+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE UnboxedTuples #-}++module Primitive (primitiveBoundingRadius, + primitiveClosestIntersect, + primitiveAnyIntersect,+ primitiveTangentSpace, + Object(Object), + Primitive(Sphere, Plane, TriangleMesh), + primitive, + material, + makeTriangle, + makeQuad, + quadsToTriangles,+ vertPosition, + vertUV, + vertTangentSpace, + transform, + radius,+ triangles, + getCentre, + planeDistance, + primitiveBoundingBox, + objectListBoundingBox, + intersectsBox, + infinite, + boundingBoxValid, + sphereIntersect,+ TangentSpace,+ Vertex) where++import Ray+import Vector+import Material+import Matrix+import BoundingBox+import Misc+import Data.Maybe+import Data.List++-- Triangle object used for triangle meshes+data Vertex = Vertex { vertPosition :: {-# UNPACK #-} !Position, + vertUV :: {-# UNPACK #-} !Position, + vertTangentSpace :: {-# UNPACK #-} !TangentSpace } deriving (Show, Eq)+data Triangle = Triangle { vertices :: ![Vertex], plane :: !Primitive, halfPlanes :: ![Primitive] } deriving (Show, Eq)++-- General object definition+data Object = Object { primitive :: Primitive,+ material :: Material,+ transform :: !Matrix} deriving (Show, Eq)++-- Different kinds of primitives that an object can have+data Primitive = Sphere { radius :: {-# UNPACK #-} !Double }+ | Plane { planeTangentSpace :: {-# UNPACK #-} !TangentSpace, planeDistance :: {-# UNPACK #-} !Double }+ | TriangleMesh { triangles :: [Triangle] } deriving (Show, Eq)++-- Get the centre of an object+getCentre :: Object -> Vector+getCentre !object = getTranslation $ transform object++-- Surface normal for 3 points+surfaceNormal :: Position -> Position -> Position -> Direction+surfaceNormal !v1 !v2 !v3 = (v2 - v1) `cross` (v3 - v1)++-- Make a plane+makePlane :: Position -> Position -> Position -> Primitive+makePlane !v1 !v2 !v3 = Plane (tangent, binormal, normal) (-(v1 `dot3` normal))+ where + !normal = normalise (surfaceNormal v1 v2 v3)+ !tangent = normalise (v2 - v1)+ !binormal = normalise (v3 - v1)++-- -------------------------------------------------------------------------------------------------------------------------------------------------------------------+-- Triangle base functionality++-- Make a triangle+makeTriangle :: Position -> Position -> Position -> Triangle+makeTriangle !v1 !v2 !v3 = Triangle verts newPlane newHalfPlanes+ where newPlane = makePlane v1 v2 v3+ newTanSpace = planeTangentSpace newPlane+ verts = map (\v -> Vertex v zeroVector newTanSpace) [v1, v2, v3]+ edgeVertices = [v1, v2, v3]+ edges = map normalise [v2 - v1, v3 - v2, v1 - v3]+ edgeNormals = map (\edge -> normalise $ thr newTanSpace `cross` edge) edges+ -- TODO - The two vectors passed here are just dummies but they can fairly easily be derived+ newHalfPlanes = zipWith (\edgeNormal edgeVertex -> Plane (Vector 1 0 0 1, Vector 0 1 0 1, edgeNormal) (-(edgeNormal `dot3` edgeVertex))) edgeNormals edgeVertices++makeQuad :: [Position] -> [Triangle]+makeQuad [vert1, vert2, vert3, vert4] = [makeTriangle vert1 vert2 vert3, makeTriangle vert1 vert3 vert4]+makeQuad _ = error "makeQuad: List was invalid size"++-- Turn a list of quad vertices into a triangle list+quadsToTriangles :: [Position] -> [Triangle]+quadsToTriangles positions = quadsToTriangles' positions []+ where+ quadsToTriangles' verts currentList + | length verts >= 4 = quadsToTriangles' (drop 4 verts) (currentList ++ makeQuad (take 4 verts))+ | otherwise = currentList++-- Area of a triangle+triangleArea :: Position -> Position -> Position -> Double+triangleArea !v1 !v2 !v3 = 0.5 * magnitude (surfaceNormal v1 v2 v3)++-- Calculate the barycentric co-ordinates of a point on a triangle+calculateBarycentricCoordinates :: Position -> Triangle -> (Double, Double, Double)+calculateBarycentricCoordinates !pos !triangle = (alpha, beta, gamma)+ where !alpha = triangleArea pos v2 v3 / area+ !beta = triangleArea pos v1 v3 / area+ !gamma = 1 - alpha - beta+ !area = triangleArea v1 v2 v3+ [!v1, !v2, !v3] = map vertPosition (vertices triangle)++-- Distance to a plane+distanceToPlane :: Primitive -> Vector -> Double+distanceToPlane (Plane !(_, _, norm) !dist) !pos = (pos `dot3` norm) + dist+distanceToPlane _ _ = error "distanceToPlane: Unsupported primitive for this function"++-- Use halfplanes to test if a point is inside a triangle+pointInsideTriangle :: Triangle -> Position -> Bool+pointInsideTriangle !tri !point = all (\pln -> distanceToPlane pln point >= 0) (halfPlanes tri)++-- Intersect a ray with a triangle+intersectRayTriangle :: Ray -> Object -> Triangle -> Bool -> (# Bool, Double, Triangle #)+intersectRayTriangle !ray !obj !triangle !doubleSided+ | not doubleSided && direction ray `dot3` (thr . planeTangentSpace . plane) triangle > 0 = (# False, 0, triangle #)+ | otherwise = case primitiveClosestIntersect (plane triangle) ray obj of+ Nothing -> (# False, 0, triangle #)+ Just (!dist', _) -> if pointInsideTriangle triangle (pointAlongRay ray dist')+ then (# True, dist', triangle #)+ else (# False, 0, triangle #)++-- Intersect against a list of triangles+intersectRayTriangleList :: [Triangle] -> Int -> Maybe (Double, Int) -> Ray -> Object -> Maybe (Double, Int)+intersectRayTriangleList !(x:xs) !index !currentResult !currentRay !obj = intersectRayTriangleList xs (index + 1) newResult newRay obj+ where+ (!newRay, !newResult) = case intersectRayTriangle currentRay obj x False of+ (# False, _, _ #) -> (currentRay, currentResult)+ (# True, !dist, _ #) -> (shortenRay currentRay dist, Just (dist, index))+intersectRayTriangleList [] _ !currentResult _ _ = currentResult++-- Intersect against any triangle+intersectRayAnyTriangleList :: [Triangle] -> Int -> Ray -> Object -> Maybe (Double, Int)+intersectRayAnyTriangleList !(x:xs) !index !ray !obj = case intersectRayTriangle ray obj x True of+ (# False, _, _ #) -> intersectRayAnyTriangleList xs (index + 1) ray obj+ (# True, !dist, _ #) -> Just (dist, index)+intersectRayAnyTriangleList [] _ _ _ = Nothing++-- Get the interpolated vertex normal+interpolatedTangentSpace :: Triangle -> Double -> Double -> Double -> TangentSpace+interpolatedTangentSpace !triangle !triAlpha !triBeta !triGamma = (tangent, binormal, normal)+ where [!(tan1, bi1, norm1), !(tan2, bi2, norm2), !(tan3, bi3, norm3)] = map vertTangentSpace (vertices triangle)+ tangent = normalise $ (tan1 <*> triAlpha) + (tan2 <*> triBeta) + (tan3 <*> triGamma)+ binormal = normalise $ (bi1 <*> triAlpha) + (bi2 <*> triBeta) + (bi3 <*> triGamma)+ normal = normalise $ (norm1 <*> triAlpha) + (norm2 <*> triBeta) + (norm3 <*> triGamma)++-- -------------------------------------------------------------------------------------------------------------------------------------------------------------------+-- Family of intersection functions+primitiveClosestIntersect :: Primitive -> Ray -> Object -> Maybe (Double, Int)++-- This function intersects a ray with a sphere, and returns the closest intercept+primitiveClosestIntersect (Sphere !sphereRadius) (Ray !rayOrg !rayDir !rayLen) obj+ | discriminant < 0 = Nothing+ | discriminant == 0 = Just (((-b) ** 0.5), 0)+ | root1 >= 0 && root1 <= rayLen = Just (root1, 0)+ | root2 >= 0 && root2 <= rayLen = Just (root2, 0)+ | otherwise = Nothing+ where + !delta = rayOrg - getCentre obj+ !b = 2 * (delta `dot3` rayDir)+ !c = (delta `dot3` delta) - sphereRadius ** 2+ !discriminant = b ** 2 - 4 * c -- A is 1 because the ray direction is normalised+ !root1 = ((-b) - sqrt discriminant) * 0.5+ !root2 = ((-b) + sqrt discriminant) * 0.5++-- This function intersects a ray with a plane and returns the closest intercept+primitiveClosestIntersect (Plane !(_, _, planeNormal) !planeD) (Ray !rayOrg !rayDir !rayLen) _+ | dirDotNormal == 0 = Nothing+ | intercept >= 0 && intercept <= rayLen = Just (intercept, 0)+ | otherwise = Nothing+ where !dirDotNormal = rayDir `dot3` planeNormal+ !intercept = ((-planeD) - (rayOrg `dot3` planeNormal)) / dirDotNormal++-- Find intersection with a triangle mesh+primitiveClosestIntersect (TriangleMesh !tris) !ray !obj = intersectRayTriangleList tris 0 Nothing ray obj++primitiveAnyIntersect :: Primitive -> Ray -> Object -> Maybe (Double, Int)+primitiveAnyIntersect (TriangleMesh !tris) !ray !obj = intersectRayAnyTriangleList tris 0 ray obj+primitiveAnyIntersect primitive' ray obj = primitiveClosestIntersect primitive' ray obj++-- -------------------------------------------------------------------------------------------------------------------------------------------------------------------+-- Family of normal query functions+primitiveTangentSpace :: Primitive -> Int -> Position -> Object -> TangentSpace++-- Normal for a sphere+primitiveTangentSpace (Sphere !sphereRadius) _ !intersectionPoint obj = (tangent, binormal, normal)+ where+ tangent = Vector 1 0 0 0 -- This is clearly incorrect - fix this later!+ binormal = Vector 0 1 0 0+ normal = (intersectionPoint - getCentre obj) <*> (1 / sphereRadius)++primitiveTangentSpace (Plane !planeNormal _) _ _ _ = planeNormal+primitiveTangentSpace (TriangleMesh !tris) !triId !intersectionPoint _ = interpolatedTangentSpace triangle triAlpha triBeta triGamma+ where !triangle = tris !! triId+ (!triAlpha, !triBeta, !triGamma) = calculateBarycentricCoordinates intersectionPoint triangle++-- -------------------------------------------------------------------------------------------------------------------------------------------------------------------+-- Family of bounding radius functions+primitiveBoundingRadius :: Primitive -> Matrix -> Vector -> Double++primitiveBoundingRadius (Sphere sphereRadius) sphereTransform pos = sphereRadius + (pos `distance` getTranslation sphereTransform)+primitiveBoundingRadius (Plane _ _) _ _ = 0+primitiveBoundingRadius (TriangleMesh tris) _ centre = triangleListRadius 0 tris centre++triangleListRadius :: Double -> [Triangle] -> Vector -> Double+triangleListRadius maximumRadius (tri:tris) centre = triangleListRadius (Prelude.max maximumRadius triangleRadius) tris centre+ where+ radii = map (\v -> centre `distance` (vertPosition v)) (vertices tri)+ triangleRadius = foldl' Prelude.max 0 radii+triangleListRadius maximumRadius [] _ = maximumRadius++-- -------------------------------------------------------------------------------------------------------------------------------------------------------------------+-- Family of bounding box functions++-- Find the bounding box of a primitive+primitiveBoundingBox :: Primitive -> Object -> Maybe AABB+primitiveBoundingBox (Sphere sphereRadius) obj = Just (boxMin, boxMax)+ where+ boxMin = getCentre obj - Vector sphereRadius sphereRadius sphereRadius 0+ boxMax = getCentre obj + Vector sphereRadius sphereRadius sphereRadius 0+primitiveBoundingBox (Plane _ _) _ = Nothing+primitiveBoundingBox (TriangleMesh tris) obj = Just $ triangleListBoundingBox initialInvalidBox (transform obj) tris++-- Bounding box of a list of something+triangleListBoundingBox :: AABB -> Matrix -> [Triangle] -> AABB+triangleListBoundingBox currentBox transformMatrix (tri:tris) = triangleListBoundingBox (boundingBoxUnion currentBox thisTriangleBox) transformMatrix tris+ where+ worldSpaceVertices = map (transformVector transformMatrix . vertPosition) (vertices tri)+ (invalidMin, invalidMax) = initialInvalidBox+ thisTriangleBox = (foldl' Vector.min invalidMin worldSpaceVertices, foldl' Vector.max invalidMax worldSpaceVertices)+triangleListBoundingBox currentBox _ [] = currentBox++objectListBoundingBox :: [Object] -> AABB+objectListBoundingBox = foldr (boundingBoxUnion . (\obj -> fromMaybe initialInvalidBox (primitiveBoundingBox (primitive obj) obj))) initialInvalidBox++-- Does a primitive intersect a box?+-- Could maybe generalise the primitiveClosestIntersect function above via further pattern matching?+intersectsBox :: Primitive -> Matrix -> AABB -> Bool+intersectsBox (Sphere sphereRadius) matrix (boxMin, boxMax) = (centreX + sphereRadius) >= vecX boxMin && (centreX - sphereRadius) <= vecX boxMax &&+ (centreY + sphereRadius) >= vecY boxMin && (centreY - sphereRadius) <= vecY boxMax &&+ (centreZ + sphereRadius) >= vecZ boxMin && (centreZ - sphereRadius) <= vecZ boxMax+ where+ !centre = getTranslation matrix+ !centreX = vecX centre+ !centreY = vecY centre+ !centreZ = vecZ centre++intersectsBox (Plane (_, _, planeNormal) planeD) _ box = signum minDistance /= signum maxDistance+ where+ !minX = selectMinBoxComponent vecX planeNormal box+ !minY = selectMinBoxComponent vecY planeNormal box+ !minZ = selectMinBoxComponent vecZ planeNormal box+ !maxX = selectMaxBoxComponent vecX planeNormal box+ !maxY = selectMaxBoxComponent vecY planeNormal box+ !maxZ = selectMaxBoxComponent vecZ planeNormal box+ !minAgainstPlane = Vector minX minY minZ 1+ !maxAgainstPlane = Vector maxX maxY maxZ 1+ !minDistance = (planeNormal `dot3` minAgainstPlane) + planeD+ !maxDistance = (planeNormal `dot3` maxAgainstPlane) + planeD+ +intersectsBox (TriangleMesh tris) matrix box = boundingBoxOverlaps box triListBox+ where+ triListBox = triangleListBoundingBox initialInvalidBox matrix tris++-- -------------------------------------------------------------------------------------------------------------------------------------------------------------------+-- Aspect querying of objects+infinite :: Primitive -> Bool+infinite (Plane _ _) = True+infinite _ = False++-- -------------------------------------------------------------------------------------------------------------------------------------------------------------------+-- Specialised intersection code for bounding volume hierarchies+sphereIntersect :: Double -> Vector -> Ray -> Maybe Double+sphereIntersect !rad !centre (Ray !rayOrg !rayDir !rayLen)+ | (centre `distanceSq` rayOrg) < (rad * rad) = Just 0 -- Inside the sphere!+ | discriminant < 0 = Nothing+ | discriminant == 0 = Just ((-b) / 2)+ | root1 >= 0 && root1 <= rayLen = Just root1+ | root2 >= 0 && root2 <= rayLen = Just root2+ | otherwise = Nothing+ where + !delta = rayOrg - centre+ !b = 2.0 * (delta `dot3` rayDir)+ !c = (delta `dot3` delta) - rad**2+ !discriminant = b**2 - 4 * c -- A is 1 because the ray direction is normalised+ !root1 = (-b - sqrt discriminant) / 2+ !root2 = (-b + sqrt discriminant) / 2
+ app/src/Primitive.hs-boot view
@@ -0,0 +1,69 @@+-- Module for general primitives and intersections+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE UnboxedTuples #-}++module Primitive (primitiveBoundingRadius, + primitiveClosestIntersect, + primitiveAnyIntersect,+ primitiveTangentSpace, + Object(Object), + Primitive(Sphere, Plane, TriangleMesh), + primitive, + material, + makeTriangle, + makeQuad, + quadsToTriangles,+ vertPosition, + vertUV, + vertTangentSpace, + transform, + radius,+ triangles, + getCentre, + planeDistance, + primitiveBoundingBox, + objectListBoundingBox, + intersectsBox, + infinite, + boundingBoxValid, + sphereIntersect,+ TangentSpace) where++import Ray+import Vector+import Material+import Matrix+import BoundingBox++-- Triangle object used for triangle meshes+data Vertex = Vertex { vertPosition :: {-# UNPACK #-} !Position, + vertUV :: {-# UNPACK #-} !Position, + vertTangentSpace :: {-# UNPACK #-} !TangentSpace }+data Triangle = Triangle { vertices :: ![Vertex], plane :: !Primitive, halfPlanes :: ![Primitive] }++-- General object definition+data Object = Object { primitive :: Primitive,+ material :: Material,+ transform :: !Matrix}++-- Different kinds of primitives that an object can have+data Primitive = Sphere { radius :: {-# UNPACK #-} !Double }+ | Plane { planeTangentSpace :: {-# UNPACK #-} !TangentSpace, planeDistance :: {-# UNPACK #-} !Double }+ | TriangleMesh { triangles :: [Triangle] }++primitiveBoundingRadius :: Primitive -> Matrix -> Vector -> Double+primitiveBoundingBox :: Primitive -> Object -> Maybe AABB++primitiveClosestIntersect :: Primitive -> Ray -> Object -> Maybe (Double, Int)+primitiveAnyIntersect :: Primitive -> Ray -> Object -> Maybe (Double, Int)++primitiveTangentSpace :: Primitive -> Int -> Position -> Object -> TangentSpace+makeTriangle :: Position -> Position -> Position -> Triangle+makeQuad :: [Position] -> [Triangle]+quadsToTriangles :: [Position] -> [Triangle]+getCentre :: Object -> Vector+objectListBoundingBox :: [Object] -> AABB+intersectsBox :: Primitive -> Matrix -> AABB -> Bool+infinite :: Primitive -> Bool+sphereIntersect :: Double -> Vector -> Ray -> Maybe Double
+ app/src/Ray.hs view
@@ -0,0 +1,32 @@+-- Module for handling rays in a raytracer+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module Ray where++import Vector++-- For now, we're sticking to Doubles+data Ray = Ray { origin :: {-# UNPACK #-} !Position, direction :: {-# UNPACK #-} !Direction, rayLength :: {-# UNPACK #-} !Double }++-- Make a ray given the start and end position+rayWithPoints :: Position -> Position -> Ray+rayWithPoints !start !end = Ray start (normalise (end - start)) (end `distance` start)++rayWithDirection :: Position -> Direction -> Double -> Ray+rayWithDirection !start !dir !rayLen = Ray start dir rayLen++rayWithPosDir :: (Position, Direction) -> Double -> Ray+rayWithPosDir !(start, dir) !rayLen = Ray start dir rayLen++-- Given a ray and a distance, produce the point along the ray+pointAlongRay :: Ray -> Double -> Position+pointAlongRay (Ray !org !dir _) !dist = setWTo1 (madd org dir dist)++-- Given some intercept, work out if it is valid, for this ray+validIntercept :: Ray -> Double -> Bool+validIntercept (Ray _ _ !rayLen) !t = t >= 0 && t <= rayLen++-- Make a shorter version of the same ray+shortenRay :: Ray -> Double -> Ray+shortenRay (Ray !org !dir _) !newLength = Ray org dir newLength
+ app/src/RayTrace.hs view
@@ -0,0 +1,254 @@+-- The module where all the tracing actually happens+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module RayTrace (rayTraceImage, findNearestIntersection, findAnyIntersection, GlobalIlluminationFunc) where++import Vector+import {-# SOURCE #-} Light+import Primitive+import Colour+import Ray+import Material+import Matrix+import Misc+import Camera+import Distribution+import SceneGraph+import Control.Parallel.Strategies+import {-# SOURCE #-} PhotonMap (PhotonMap, irradiance)+import IrradianceCache+import Control.Monad.State+import RenderContext++-- Intersect a ray against a sphere tree+intersectSphereTree :: [SphereTreeNode] -> Ray -> Maybe (Object, Double, Int) -> Maybe (Object, Double, Int)+intersectSphereTree !(node:nodes) !ray !currentHit = seq result (intersectSphereTree (newNodeList ++ nodes) newRay thisResult)+ where+ -- Intersect the ray with the bounding volume of this node+ !result = case children node of+ -- If the node has no children, don't bother with it's bounding volume and just check the object (if it has one)+ [] -> case object node of+ Nothing -> error "A node with no children should hold an object"+ Just obj -> case primitiveClosestIntersect (primitive obj) ray obj of+ -- Didn't hit the object. Retain the current hit, and continue with remaining nodes on the list+ Nothing -> (currentHit, [], ray)+ -- We did hit this object. Update the intersection, and continue with remaining nodes on the list+ Just (objHitDistance, objHitId) -> (Just (obj, objHitDistance, objHitId), [], shortenRay ray objHitDistance)+ -- We have children. In this case it makes sense to test our bounding volume+ !nodeChildren -> case sphereIntersect (boundingRadius node) (boundingCentre node) ray of -- (make a sphere centred at the object's transform matrix with given radius)+ -- If we do not find an intersection, we do not update the results and we offer no further nodes to be traversed, thus skipping this subtree+ Nothing -> (currentHit, [], ray)+ -- If we do find an intersection against the bounding volume, then we try again against the actual object (if present)+ Just _ -> case object node of+ Nothing -> (currentHit, nodeChildren, ray) -- No object; just pass to the children+ Just obj -> case primitiveClosestIntersect (primitive obj) ray obj of+ -- Didn't hit the object. Retain the current hit, but offer up the children of the node as we hit the bounding volume+ Nothing -> (currentHit, nodeChildren, ray)+ -- We did hit this object. Update the intersection, and continue with the bounding volume's children+ Just (objHitDistance, objHitId) -> (Just (obj, objHitDistance, objHitId), nodeChildren, shortenRay ray objHitDistance)+ (!thisResult, !newNodeList, !newRay) = result+intersectSphereTree [] _ !currentHit = currentHit++-- Intersect with the list of infinite objects+intersectObjectList :: [Object] -> Ray -> Maybe (Object, Double, Int) -> Maybe (Object, Double, Int)+intersectObjectList !(obj:objs) !ray !currentHit = intersectObjectList objs newRay thisResult+ where+ (!thisResult, !newRay) = case primitiveClosestIntersect (primitive obj) ray obj of+ Nothing -> (currentHit, ray)+ Just (objHitDistance, objHitId) -> (Just (obj, objHitDistance, objHitId), shortenRay ray objHitDistance)+intersectObjectList [] _ !currentHit = currentHit++-- Find the nearest intersection along a line+findNearestIntersection :: SceneGraph -> Ray -> Maybe (Object, Double, Int)+findNearestIntersection sceneGraph' !ray = case intersectObjectList (infiniteObjects sceneGraph') ray Nothing of+ Just (obj, dist, objId) -> intersectSphereTree [root sceneGraph'] (shortenRay ray dist) (Just (obj, dist, objId))+ Nothing -> intersectSphereTree [root sceneGraph'] ray Nothing++-- Intersect a ray against a scene graph. Return first (ie, any) hit without finding the closest+findAnyIntersectionSphereTree :: [SphereTreeNode] -> Ray -> Maybe (Object, Double)+findAnyIntersectionSphereTree (node:nodes) !ray = let sphereIntersectionResult = sphereIntersect (boundingRadius node) (boundingCentre node) ray+ in case sphereIntersectionResult of -- (make a sphere centred at the object's transform matrix with given radius)+ -- If we do not find an intersection, traverse to the rest of the list+ Nothing -> findAnyIntersectionSphereTree nodes ray+ -- If we do find an intersection against the bounding volume, then we try again against the actual object+ Just _ -> case object node of+ Nothing -> findAnyIntersectionSphereTree (nodes ++ children node) ray -- No object here - just offer up the children+ Just obj -> case primitiveAnyIntersect (primitive obj) ray obj of + -- Didn't hit the object. Offer up the children of the scene graph node to continue with (as we did actually hit the bounding volume)+ Nothing -> findAnyIntersectionSphereTree (nodes ++ children node) ray+ -- We did hit this object. Update the intersection, and continue with the bounding volume's children+ Just (objHitDistance, _) -> Just (obj, objHitDistance)+findAnyIntersectionSphereTree [] _ = Nothing++-- Find any intersection against an object list+findAnyIntersectionObjectList :: [Object] -> Ray -> Maybe (Object, Double)+findAnyIntersectionObjectList (obj:objs) !ray = case primitiveAnyIntersect (primitive obj) ray obj of+ -- Didn't hit he object. Retain the current hit, but offer up the children of the node as we hit the bounding volume+ Nothing -> findAnyIntersectionObjectList objs ray+ -- We did hit this object. Update the intersection, and continue with the bounding volume's children+ Just (objHitDistance, _) -> Just (obj, objHitDistance)+findAnyIntersectionObjectList [] _ = Nothing++findAnyIntersection :: SceneGraph -> Ray -> Maybe (Object, Double)+findAnyIntersection sceneGraph' !ray = case findAnyIntersectionObjectList (infiniteObjects sceneGraph') ray of+ Nothing -> findAnyIntersectionSphereTree [root sceneGraph'] ray+ Just x -> Just x++-- Default background colour to return when we can't match anything+defaultColour :: Direction -> Colour+defaultColour _ = colBlue++-- Accumulate the contributions of the lights+lightSurface :: [Light] -> Colour -> RenderContext -> SurfaceLocation -> Material -> Vector -> Colour+lightSurface (x:xs) !acc renderContext !posTanSpace !objMaterial !viewDirection + = let result = acc + applyLight (sceneGraph renderContext) posTanSpace objMaterial viewDirection x+ in seq result (lightSurface xs result renderContext posTanSpace objMaterial viewDirection)+lightSurface [] !acc _ _ _ _ = acc++-- Magic number for the usable radius of an irradaiance cache sample+irrCacheSampleRadius :: Double+irrCacheSampleRadius = 10++-- Abstraction to permit different GI calculations+type GlobalIlluminationFunc = (SurfaceLocation -> IrradianceCache -> Object -> RenderContext -> (Colour, IrradianceCache))++-- Photon map specific GI calculator+photonMapGlobalIllumination :: Maybe PhotonMap -> SurfaceLocation -> IrradianceCache -> Object -> RenderContext -> (Colour, IrradianceCache)+photonMapGlobalIllumination (Just photonMap) !surfaceLocation irrCache obj renderContext = + case renderMode renderContext of+ PhotonMapper -> query irrCache surfaceLocation irradiance'+ _ -> undefined -- Shouldn't hit this path...+ where+ irradiance' x = (irradiance photonMap (photonMapContext renderContext) (material obj) x, irrCacheSampleRadius)+photonMapGlobalIllumination _ _ irrCache _ _ = (colBlack, irrCache)++-- Null GI+nullGI :: SurfaceLocation -> IrradianceCache -> Object -> RenderContext -> (Colour, IrradianceCache)+nullGI _ irrCache _ _ = (colBlack, irrCache)++-- Retrieve the appropriate GI function and calculate the GI at a point+calculateGI :: RenderContext -> Maybe PhotonMap -> GlobalIlluminationFunc+calculateGI renderContext photonMap = case renderMode renderContext of+ PhotonMapper -> photonMapGlobalIllumination photonMap+ _ -> nullGI++-- Perform a full trace of a ray+type RayTraceState = State IrradianceCache Colour+traceRay :: RenderContext -> Maybe PhotonMap -> Ray -> Int -> Direction -> Double -> Double -> RayTraceState++-- Special case for lowest level of recursion (theoretically this should not get hit)+traceRay _ _ _ 0 _ _ _ = error "Should not hit this codepath"++-- Special case for penultimate level - we're not allowed to spawn rays here+traceRay renderContext photonMap !ray 1 !viewDir _ _ = + case findNearestIntersection (sceneGraph renderContext) ray of+ Nothing -> return $! defaultColour (direction ray)+ Just (obj, intersectionDistance, hitId) -> do+ irrCache <- get+ let !intersectionPoint = pointAlongRay ray intersectionDistance+ let !tanSpace = primitiveTangentSpace (primitive obj) hitId intersectionPoint obj+ let !(!surfaceIrradiance, !newIrrCache) = calculateGI renderContext photonMap (intersectionPoint, tanSpace) irrCache obj renderContext+ -- TODO - Need to plug irradiance values into surface shading more correctly+ let resultColour = lightSurface (lights renderContext) surfaceIrradiance renderContext (intersectionPoint, tanSpace) (material obj) viewDir+ put newIrrCache+ return $! resultColour++-- General case+traceRay renderContext photonMap !ray !limit !viewDir !currentIOR !accumulatedReflectivity = + case findNearestIntersection (sceneGraph renderContext) ray of+ Nothing -> return $! defaultColour (direction ray)+ Just (obj, intersectionDistance, hitId) -> do+ -- Evaluate surface-location specific things such as shader results+ let !intersectionPoint = pointAlongRay ray intersectionDistance+ let !tanSpace = primitiveTangentSpace (primitive obj) hitId intersectionPoint obj+ let !normal = thr tanSpace+ let !incoming = Vector.negate $ direction ray+ -- TODO - Need to plug irradiance values into surface shading more correctly++ -- Evaluate result from irradiance cache+ irrCache <- get+ let !(!surfaceIrradiance, !irrCache') = calculateGI renderContext photonMap (intersectionPoint, tanSpace) irrCache obj renderContext+ put irrCache'++ let !surfaceShading = lightSurface (lights renderContext) surfaceIrradiance renderContext (intersectionPoint, tanSpace) (material obj) viewDir++ -- Reflection specific code+ let offsetToExterior = madd intersectionPoint normal surfaceEpsilon+ let reflectionDir = normalise $ reflect incoming normal+ let !shine = reflectivity $ material obj+ let reflectRay = rayWithDirection offsetToExterior reflectionDir (reflectionRayLength renderContext)+ let (reflection, irrCache'') = if shine > 0 && (accumulatedReflectivity * shine) > 0.03 + then runState (traceRay renderContext photonMap reflectRay (limit - 1) viewDir currentIOR (accumulatedReflectivity * shine)) irrCache'+ else (colBlack, irrCache')+ put irrCache''++ -- Refraction specific+ let !eta = if enteringObject incoming normal+ then currentIOR / indexOfRefraction (material obj) + else indexOfRefraction (material obj) / currentIOR+ let refractionDir = normalise $ refract incoming normal eta+ let offsetToInterior = madd intersectionPoint refractionDir surfaceEpsilon+ let !transmittance = transmit $ material obj+ let refractRay = rayWithDirection offsetToInterior refractionDir (refractionRayLength renderContext)+ let (refraction, irrCache''') = if transmittance > 0 + then runState (traceRay renderContext photonMap refractRay (limit - 1) viewDir (indexOfRefraction $ material obj) accumulatedReflectivity) irrCache''+ else (colBlack, irrCache'')+ put irrCache'''++ -- Final colour combine+ return $! (surfaceShading + (reflection Colour.<*> shine) + (refraction Colour.<*> transmittance))+ where+ enteringObject !incoming !normal = incoming `dot3` normal > 0++-- This function converts a pixel co-ordinate to a direction of the ray+makeRayDirection :: Int -> Int -> Camera -> (Int, Int) -> Vector+makeRayDirection !renderWidth !renderHeight !camera (x, y) =+ let x' = (fromIntegral x / fromIntegral renderWidth) * 2.0 - 1.0+ y' = (fromIntegral y / fromIntegral renderHeight) * 2.0 - 1.0+ fov = 0.5 * fieldOfView camera+ fovX = tan (degreesToRadians fov)+ fovY = -tan (degreesToRadians fov)+ aspectRatio = fromIntegral renderWidth / fromIntegral renderHeight+ !dirX = fovX * x'+ !dirY = fovY * (-y') / aspectRatio+ rayDir = normalise (Vector dirX dirY 1 0)+ in normalise $ transformVector (worldToCamera camera) rayDir++-- Trace a list of distributed samples with tail recursion+traceDistributedSample :: RenderContext -> Colour -> [Position] -> Maybe PhotonMap -> (Position, Direction) -> Double -> RayTraceState+traceDistributedSample renderContext !acc (x:xs) photonMap !eyeViewDir !sampleWeighting = + do+ irrCache <- get+ let !dofFocalDistance = depthOfFieldFocalDistance renderContext+ let jitteredRayPosition jitter = fst eyeViewDir + jitter+ let jitteredRayDirection jitter = normalise $ madd jitter (snd eyeViewDir) dofFocalDistance+ let (sampleColour, irrCache') = runState (traceRay renderContext photonMap (rayWithDirection (jitteredRayPosition x) (jitteredRayDirection x) 100000.0) (maximumRayDepth renderContext) (snd eyeViewDir) 1 1) irrCache+ let result = (sampleColour Colour.<*> sampleWeighting) + acc+ put irrCache'+ let (col, irrCache'') = runState (traceDistributedSample renderContext result xs photonMap eyeViewDir sampleWeighting) irrCache'+ put irrCache''+ return $! col+traceDistributedSample _ !acc [] _ _ _ = return $! acc++-- Need to remove hard coded constants of 8 here+-- This traces for a given pixel (x, y)+tracePixel :: RenderContext -> Position -> Maybe PhotonMap -> Direction -> RayTraceState+tracePixel renderContext !eye photonMap !viewDirection = do+ irrCache <- get+ let !distributedPositions = generatePointsOnSphere (numDistribSamples renderContext) (rayOriginDistribution renderContext) 12345+ let (!pixelColour, irrCache') = runState (traceDistributedSample renderContext colBlack distributedPositions photonMap (eye, viewDirection) (1.0 / (fromIntegral . numDistribSamples $ renderContext))) irrCache+ put irrCache'+ return $! pixelColour++-- Generate a list of colours which contains a raytraced image. In parallel+rayTraceImage :: RenderContext -> Camera -> Int -> Int -> Maybe PhotonMap -> [Colour]+rayTraceImage renderContext camera renderWidth renderHeight photonMap = tracePixelPassingState rayDirections irrCache `using` parListChunk 256 rdeepseq+ where !rayDirections = [makeRayDirection renderWidth renderHeight camera (x, y) | y <- [0..(renderHeight - 1)], x <- [0..(renderWidth - 1)]]+ !eyePosition = Camera.position camera+ irrCache = initialiseCache (sceneGraph renderContext)+ -- This function is the equivalent to map, but it passes the ending state of one invocation to the next invocation+ tracePixelPassingState (x:xs) st = result : tracePixelPassingState xs st'+ where+ (!result, !st') = runState (tracePixel renderContext eyePosition photonMap x) st+ tracePixelPassingState [] _ = []
+ app/src/RayTrace.hs-boot view
@@ -0,0 +1,22 @@+module RayTrace (rayTraceImage, findNearestIntersection, findAnyIntersection, RenderContext(RenderContext)) where++import Vector+import Light+import Primitive+import Colour+import Ray+import Material+import Matrix+import Camera+import SceneGraph++data RenderContext = RenderContext {+ numDistribSamples :: Int,+ photonGatherDistance :: Double,+ sceneGraph :: SceneGraph,+ lights :: [Light],+ maximumRayDepth :: Int }++rayTraceImage :: RenderContext -> Camera -> Int -> Int -> PhotonMap -> [Colour]+findNearestIntersection :: SceneGraph -> Ray -> Maybe (Object, Double, Int)+findAnyIntersection :: SceneGraph -> Ray -> Maybe (Object, Double)
+ app/src/RenderContext.hs view
@@ -0,0 +1,23 @@+-- This is a render context, something that describes the general shared variables for rendering++module RenderContext where++import SceneGraph+import {-# SOURCE #-} Light+import {-# SOURCE #-} PhotonMap (PhotonMapContext)++data RenderMode = RayTrace | PhotonMapper | PathTracer deriving (Show)++data RenderContext = RenderContext {+ numDistribSamples :: Int,+ sceneGraph :: SceneGraph,+ lights :: [Light],+ maximumRayDepth :: Int,+ reflectionRayLength :: Double,+ refractionRayLength :: Double,+ photonMapContext :: PhotonMapContext,+ rayOriginDistribution :: Double,+ depthOfFieldFocalDistance :: Double,+ renderMode :: RenderMode,+ useIrradianceCache :: Bool }+
+ app/src/RenderContext.hs-boot view
@@ -0,0 +1,20 @@+-- This is a render context, something that describes the general shared variables for rendering++module RenderContext where++import SceneGraph++data RenderMode = RayTrace | PhotonMapper | PathTracer++data RenderContext = RenderContext {+ numDistribSamples :: Int,+ sceneGraph :: SceneGraph,+ lights :: [Light],+ maximumRayDepth :: Int,+ reflectionRayLength :: Double,+ refractionRayLength :: Double,+ photonMapContext :: PhotonMapContext,+ rayOriginDistribution :: Double,+ depthOfFieldFocalDistance :: Double,+ renderMode :: RenderMode,+ useIrradianceCache :: Bool }
+ app/src/SceneGraph.hs view
@@ -0,0 +1,48 @@+-- The graph structure holding the scene+{-# LANGUAGE MagicHash #-}++module SceneGraph (buildSceneGraph, SphereTreeNode(boundingRadius, boundingCentre, object, children), SceneGraph(root, infiniteObjects, finiteBox)) where++import Primitive+import Vector+import BoundingBox+import Data.List++data SphereTreeNode = SphereTreeNode { object :: Maybe Object, children :: [SphereTreeNode], boundingRadius :: !Double, boundingCentre :: !Vector } deriving (Show)+data SceneGraph = SceneGraph { root :: SphereTreeNode, infiniteObjects :: [Object], finiteBox :: AABB } deriving (Show)++-- Find the mean of a collection of objects+calculateMeanPosition' :: [Object] -> Vector -> Vector+calculateMeanPosition' (obj : objects) acc = calculateMeanPosition' objects acc + getCentre obj+calculateMeanPosition' [] acc = acc++calculateMeanPosition :: [Object] -> Vector+calculateMeanPosition objects = setWTo1 (calculateMeanPosition' objects zeroVector </> len')+ where+ !len' = fromIntegral (length objects)++-- Find the overall bounding radius of a list of objects+calculateBoundingRadius :: [Object] -> Vector -> Double+calculateBoundingRadius objs centre = foldr (Prelude.max . (\obj -> primitiveBoundingRadius (primitive obj) (transform obj) centre)) 0 objs++-- Build up a sphere tree+buildSphereTree :: ([Object] -> [[Object]]) -> [Object] -> SphereTreeNode+buildSphereTree _ (obj : []) = SphereTreeNode (Just obj) [] nodeRadius nodeCentre+ where+ nodeCentre = calculateMeanPosition [obj]+ nodeRadius = calculateBoundingRadius [obj] nodeCentre+buildSphereTree builder (obj:objs)+ | length (obj:objs) == 1 = error "Should have been handled by a different pattern"+ | null (obj:objs) = error "Should not have zero objects"+ | otherwise = SphereTreeNode Nothing nodeChildren nodeRadius nodeCentre+ where+ nodeCentre = calculateMeanPosition (obj:objs)+ nodeRadius = calculateBoundingRadius (obj:objs) nodeCentre+ nodeChildren = map (buildSphereTree builder) (builder (obj:objs))+buildSphereTree _ [] = error "Should not hit this pattern for buildSphereTree" ++-- Build a scene graph+buildSceneGraph :: [Object] -> ([Object] -> [[Object]]) -> SceneGraph+buildSceneGraph objs buildFunction = SceneGraph (buildSphereTree buildFunction nonInfiniteObjs) infiniteObjs (objectListBoundingBox nonInfiniteObjs)+ where+ (infiniteObjs, nonInfiniteObjs) = partition (infinite . primitive) objs
+ app/src/SceneGraphTest.hs view
@@ -0,0 +1,30 @@+module SceneGraphTest where++import SceneGraph+import Primitive+import Material+import Matrix+import Shader+import Colour++testObjs1 = [+ Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100)+ ]++testObjs2 = [+ Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100),+ Object (Sphere 50) (Material colGreen colGreen colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 0 100)+ ]++testObjs3 = [+ Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100),+ Object (Sphere 50) (Material colGreen colGreen colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 0 100),+ Object (Sphere 50) (Material colBlue colBlue colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 200 100)+ ]++testObjs4 = [+ Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100),+ Object (Sphere 50) (Material colGreen colGreen colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 0 100),+ Object (Sphere 50) (Material colBlue colBlue colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 200 100),+ Object (Sphere 50) (Material colYellow colYellow colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 200 100)+ ]
+ app/src/Shader.hs view
@@ -0,0 +1,50 @@+-- Generic shaders to return colour and texture information+{-# LANGUAGE MagicHash #-}++module Shader where++import Vector+import Colour+import Misc++data Shader = CheckedShader { scale :: {-# UNPACK #-} !Vector, colour1 :: {-# UNPACK #-} !Colour, colour2 :: {-# UNPACK #-} !Colour }+ | ShowNormalShader+ | NullShader deriving (Show, Eq)++-- Functions available for each shader+evaluateAmbient :: Shader -> Position -> TangentSpace -> Colour+evaluateDiffuse :: Shader -> Position -> TangentSpace -> Colour+evaluateSpecular :: Shader -> Position -> TangentSpace -> Colour+shadePoint :: Shader -> (Position, Direction) -> (Colour, Colour, Colour) -> Colour++-- Checked shaders+evaluateDiffuse (CheckedShader checkScale checkColour1 checkColour2) position _ = + let scaledPosition = checkScale * position+ scaledX = round (vecX scaledPosition) :: Int+ scaledY = round (vecY scaledPosition) :: Int+ scaledZ = round (vecZ scaledPosition) :: Int+ in if odd scaledX `xor` odd scaledY `xor` odd scaledZ then checkColour1 else checkColour2++-- Normal display+evaluateDiffuse ShowNormalShader _ (_, _, normal) = encodeNormal normal++-- Null shader+evaluateDiffuse NullShader _ _ = colWhite++-- Defaults+evaluateSpecular = evaluateDiffuse+evaluateAmbient = evaluateDiffuse++-- New style shader interface+shadePoint (CheckedShader checkScale checkColour1 checkColour2) (position, _) (ambient, diffuse, specular) = (ambient + diffuse + specular) * checkColour+ where+ scaledPosition = checkScale * position+ scaledX = round (vecX scaledPosition) :: Int+ scaledY = round (vecY scaledPosition) :: Int+ scaledZ = round (vecZ scaledPosition) :: Int+ checkColour = if odd scaledX `xor` odd scaledY `xor` odd scaledZ then checkColour1 else checkColour2++shadePoint ShowNormalShader (_, norm) (_, _, _) = encodeNormal norm++-- Default fallback+shadePoint _ (_, _) (ambient, diffuse, specular) = ambient + diffuse + specular
+ app/src/Tests/BoundingBoxTest.hs view
@@ -0,0 +1,22 @@+module Tests.BoundingBoxTest where++import Test.HUnit+import Vector+import BoundingBox+import GHC.Prim+import GHC.Types++test_InsideBox = TestCase (assertEqual "Inside box" (box `contains` pos) True)+ where+ box = (Vector (-5) (-5) (-5) 1, Vector 5 5 5 1)+ pos = Vector 0 0 0 1++test_OutsideBoxMinX = TestCase (assertEqual "Inside box" (box `contains` pos) False)+ where+ box = (Vector (-5) (-5) (-5) 1, Vector 5 5 5 1)+ pos = Vector (-10) 0 0 1++tests_BoundingBox = TestList [+ TestLabel "InsideBox" test_InsideBox,+ TestLabel "OutsideBoxMinX" test_OutsideBoxMinX+ ]
+ app/src/Tests/ColourTest.hs view
@@ -0,0 +1,84 @@+module Tests.ColourTest where++import Colour+import Test.HUnit++test_Add = TestCase (assertEqual "Colour addition" expected (v1 + v2))+ where+ v1 = Colour 1 2 3 4+ v2 = Colour 10 20 30 40+ expected = Colour 11 22 33 44++test_Sub = TestCase (assertEqual "Colour subtraction" expected (v1 - v2))+ where+ v1 = Colour 10 20 30 40+ v2 = Colour 5 10 15 20+ expected = Colour 5 10 15 20++test_Mul = TestCase (assertEqual "Colour multiplication" expected (v1 * v2))+ where+ v1 = Colour 10 20 30 40+ v2 = Colour 5 10 15 20+ expected = Colour 50 200 450 800++test_Abs = TestCase (assertEqual "Colour abs" expected (abs v1))+ where+ v1 = Colour 10 (-20) 0 (-40)+ expected = Colour 10 20 0 40++test_Signum = TestCase (assertEqual "Colour signum" expected (signum v1))+ where+ v1 = Colour 10 (-20) 0 (-40)+ expected = Colour 1 (-1) 0 (-1)++test_fromInteger = TestCase (assertEqual "Colour fromInteger" expected (signum v1))+ where+ v1 = fromInteger 5+ expected = Colour 5 5 5 5++test_Div = TestCase (assertEqual "Colour division" expected (v1 / v2))+ where+ v1 = Colour 10 20 30 40+ v2 = Colour 5 10 15 20+ expected = Colour 2 2 2 2++test_AddScalar = TestCase (assertEqual "Colour add scalar" expected (v1 <+> k))+ where+ v1 = Colour 10 20 30 40+ k = 1+ expected = Colour 11 21 31 41++test_SubScalar = TestCase (assertEqual "Colour sub scalar" expected (v1 <-> k))+ where+ v1 = Colour 10 20 30 40+ k = 1+ expected = Colour 9 19 29 39++test_Clamp = TestCase (assertEqual "Colour clamp" expected (clamp v1))+ where+ v1 = Colour 2 0 (-5) 0.5+ expected = Colour 1 0 0 0.5++test_Luminance0 = TestCase (assertEqual "Luminance = 0" expected (luminance v1))+ where+ v1 = Colour 0 0 0 0+ expected = 0++test_Luminance1 = TestCase (assertEqual "Luminance = 1" expected (luminance v1))+ where+ v1 = Colour 1 1 1 1+ expected = 0.9999999999999999++tests_Colour = TestList [+ TestLabel "Addition" test_Add,+ TestLabel "Subtraction" test_Sub,+ TestLabel "Multiplication" test_Mul,+ TestLabel "Abs" test_Abs,+ TestLabel "Signum" test_Signum,+ TestLabel "Div" test_Div,+ TestLabel "Add scalar" test_AddScalar,+ TestLabel "Sub scalar" test_SubScalar,+ TestLabel "Clamp" test_Clamp,+ TestLabel "Luminance = 0" test_Luminance0,+ TestLabel "Luminance = 1" test_Luminance1+ ]
+ app/src/Tests/UnitTests.hs view
@@ -0,0 +1,10 @@+-- This defines all unit tests to be executed for this project++import Tests.VectorTest+import Tests.BoundingBoxTest+import Tests.ColourTest+import Test.HUnit++unitTests = [tests_Vector, tests_BoundingBox, tests_Colour]++main = mapM runTestTT unitTests
+ app/src/Tests/VectorTest.hs view
@@ -0,0 +1,144 @@+module Tests.VectorTest where++import Vector+import Test.HUnit++test_Add = TestCase (assertEqual "Vector addition" expectedResult (v1 + v2))+ where+ v1 = Vector 1 2 3 4+ v2 = Vector 10 20 30 40+ expectedResult = Vector 11 22 33 44++test_Sub = TestCase (assertEqual "Vector subtraction" expectedResult (v1 - v2))+ where+ v1 = Vector 10 20 30 40+ v2 = Vector 1 2 3 4+ expectedResult = Vector 9 18 27 36++test_Mul = TestCase (assertEqual "Vector multiplication" expectedResult (v1 * v2))+ where+ v1 = Vector 1 0 2 3+ v2 = Vector 1 10 (-2) 3+ expectedResult = Vector 1 0 (-4) 9++test_Negate = TestCase (assertEqual "Vector negation" expectedResult (Vector.negate v1))+ where+ v1 = Vector 1 (-2) 3 (-4)+ expectedResult = Vector (-1) 2 (-3) 4++test_Abs = TestCase (assertEqual "Vector abs" expectedResult (abs v1))+ where+ v1 = Vector 0 1 (-2) 3+ expectedResult = Vector 0 1 2 3++test_Signum = TestCase (assertEqual "Vector signum" expectedResult (signum v1))+ where+ v1 = Vector 1 (-1) 0 (-2)+ expectedResult = Vector 1 (-1) 0 (-1)++test_Madd = TestCase (assertEqual "Vector madd" expectedResult (madd pos dir k))+ where+ pos = Vector 1 2 3 1+ dir = Vector 0.5 0 1 0+ k = 10+ expectedResult = Vector 6 2 13 1++test_ScalarMul = TestCase (assertEqual "Vector-scalar mul" expectedResult (vectorScalarMul vec k))+ where+ vec = Vector 1 2 (-3) 1+ k = 2+ expectedResult = Vector 2 4 (-6) 2++test_ScalarDiv = TestCase (assertEqual "Vector-scalar div" expectedResult (vec </> k))+ where+ vec = Vector 10 20 (-30) 40+ k = 2+ expectedResult = Vector 5 10 (-15) 20++test_ScalarDot3 = TestCase (assertEqual "dot3" expectedResult (v1 `dot3` v2))+ where+ v1 = Vector 1 2 0 1+ v2 = Vector (-2) 4 (-5) 1+ expectedResult = 6++test_ScalarDot4 = TestCase (assertEqual "dot4" expectedResult (v1 `dot4` v2))+ where+ v1 = Vector 1 2 0 1+ v2 = Vector (-2) 4 (-5) 1+ expectedResult = 7++test_SatScalarDot3 = TestCase (assertEqual "sdot3" expectedResult (v1 `sdot3` v2))+ where+ v1 = Vector 1 2 0 1+ v2 = Vector (-2) 4 (-5) 1+ expectedResult = 1++test_SatScalarDot4 = TestCase (assertEqual "sdot4" expectedResult (v1 `sdot4` v2) )+ where+ v1 = Vector 1 2 0 1+ v2 = Vector (-2) (-4) (-5) 1+ expectedResult = 0++test_Cross = TestCase (assertEqual "cross" expectedResult (v1 `cross` v2))+ where+ v1 = Vector 1 0 0 0+ v2 = Vector 0 1 0 0+ expectedResult = Vector 0 0 1 0++test_Magnitude = TestCase (assertEqual "magnitude" expectedResult (magnitude v1))+ where+ v1 = Vector 3 4 0 0+ expectedResult = 5++test_MagnitudeSq = TestCase (assertEqual "magnitudeSq" expectedResult (magnitudeSq v1))+ where+ v1 = Vector 3 4 0 0+ expectedResult = 25++test_Normalise = TestCase (assertEqual "normalise" expectedResult (normalise v1))+ where+ v1 = Vector 1 (-1) 1 0+ expectedResult = Vector 0.5773502691896258 (-0.5773502691896258) 0.5773502691896258 0++test_Reflect = undefined+test_Refract = undefined++test_LargestAxis = TestCase (assertEqual "largestAxis" expectedResult (largestAxis v1))+ where+ v1 = Vector (-1) 2 (-3) 0+ expectedResult = 2++test_Min = TestCase (assertEqual "min" expectedResult (v1 `Vector.min` v2))+ where+ v1 = Vector (-1) 2 (-3) 8+ v2 = Vector 10 (-20) 50 2+ expectedResult = Vector (-1) (-20) (-3) 2++test_Max = TestCase (assertEqual "max" expectedResult (v1 `Vector.max` v2))+ where+ v1 = Vector (-1) 2 (-3) 8+ v2 = Vector 10 (-20) 50 2+ expectedResult = Vector 10 2 50 8++tests_Vector = TestList [+ TestLabel "Addition" test_Add, + TestLabel "Subtraction" test_Sub,+ TestLabel "Multiplication" test_Mul,+ TestLabel "Negation" test_Negate,+ TestLabel "Abs" test_Abs,+ TestLabel "Signum" test_Signum,+ TestLabel "Madd" test_Madd,+ TestLabel "Vector * scalar" test_ScalarMul,+ TestLabel "Vector / scalar" test_ScalarDiv,+ TestLabel "Dot3" test_ScalarDot3,+ TestLabel "Dot4" test_ScalarDot4,+ TestLabel "Sdot3" test_SatScalarDot3,+ TestLabel "Sdot4" test_SatScalarDot4,+ TestLabel "Cross" test_Cross,+ TestLabel "Magnitude" test_Magnitude,+ TestLabel "MagnitudeSq" test_MagnitudeSq,+ TestLabel "Normalise" test_Normalise,+ TestLabel "LargestAxis" test_LargestAxis,+ TestLabel "Min" test_Min,+ TestLabel "Max" test_Max+ ]
+ app/src/ToneMap.hs view
@@ -0,0 +1,61 @@+-- Tone map an image+{-# LANGUAGE BangPatterns #-}++module ToneMap(toneMapImage, + toneMapIdentity, + toneMapAverageLuminance, + toneMapReinhard, + toneMapHejlBurgessDawson,+ exposeImage,+ imageAverageLogLuminance,+ imageAverageLuminance) where++import Colour++-- x = x+toneMapIdentity :: [Colour] -> [Colour]+toneMapIdentity = map id++-- x = x / avg xs+toneMapAverageLuminance :: [Colour] -> [Colour]+toneMapAverageLuminance xs = map (<*> invAverageBrightness) xs+ where+ invAverageBrightness = 1 / imageAverageLuminance xs++-- Reinhard tone map operator http://filmicgames.com/archives/75+toneMapReinhard :: [Colour] -> [Colour]+toneMapReinhard = map (\(Colour !r !g !b _) -> Colour (r / (r + 1)) (g / (g + 1)) (b / (b + 1)) 1)++-- Hejl-Burgess-Dawson http://filmicgames.com/archives/75+toneMapHejlBurgessDawson :: [Colour] -> [Colour]+toneMapHejlBurgessDawson = map f+ where+ f colour = (x * (x <*> 6.2 <+> 0.5)) / (x * (x <*> 6.2 <+> 1.7) <+> 0.06)+ where+ x = (\x' -> fold max x' 0) (colour <-> 0.004)++-- Apply a tone map operator+toneMapImage :: ([Colour] -> [Colour]) -> [Colour] -> [Colour]+toneMapImage f = f++-- Normal averaging+imageAverageLuminance :: [Colour] -> Double+imageAverageLuminance = imageAverageLuminance' 0 0+ where+ imageAverageLuminance' accLum accCount (x:xs) = imageAverageLuminance' (accLum + luminance x) (accCount + 1) xs+ imageAverageLuminance' accLum 0 [] = accLum+ imageAverageLuminance' accLum accCount [] = accLum / accCount++-- Get the average luminance of a scene, using Reinhard style log-lum averaging to damp down the effect of outlier pixels+imageAverageLogLuminance :: [Colour] -> Double+imageAverageLogLuminance = imageAverageLogLuminance' 0 0+ where+ imageAverageLogLuminance' accLum accCount (x:xs) = imageAverageLogLuminance' (accLum + logLuminance x) (accCount + 1) xs+ imageAverageLogLuminance' accLum 0 [] = exp accLum+ imageAverageLogLuminance' accLum accCount [] = exp (accLum / accCount)++-- Adjust the exposure of an image+exposeImage :: ([Colour] -> Double) -> [Colour] -> Double -> [Colour]+exposeImage f xs exposureScale = map (</> (exposure * exposureScale)) xs+ where+ exposure = f xs
+ app/src/Vector.hs view
@@ -0,0 +1,233 @@+-- Vector library for 3d graphics++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}++module Vector where+import Data.List+import Misc+import GHC.Prim+import GHC.Types+import Control.DeepSeq++data Vector = Vector { vecX :: {-# UNPACK #-} !Double,+ vecY :: {-# UNPACK #-} !Double,+ vecZ :: {-# UNPACK #-} !Double,+ vecW :: {-# UNPACK #-} !Double } deriving (Ord, Eq)+type Position = Vector+type Direction = Vector+type Normal = Direction+type TangentSpace = (Normal, Normal, Normal)+type SurfaceLocation = (Position, TangentSpace)++instance Num Vector where+ {-# SPECIALIZE INLINE (+) :: Vector -> Vector -> Vector #-}+ (Vector !(D# x) !(D# y) !(D# z) !(D# w)) + (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x +## x') (D# $ y +## y') (D# $ z +## z') (D# $ w +## w')+ {-# SPECIALIZE INLINE (-) :: Vector -> Vector -> Vector #-}+ (Vector !(D# x) !(D# y) !(D# z) !(D# w)) - (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x -## x') (D# $ y -## y') (D# $ z -## z') (D# $ w -## w')+ {-# SPECIALIZE INLINE (*) :: Vector -> Vector -> Vector #-}+ (Vector !(D# x) !(D# y) !(D# z) !(D# w)) * (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x *## x') (D# $ y *## y') (D# $ z *## z') (D# $ w *## w')+ abs (Vector !x !y !z !w) = Vector absX absY absZ absW+ where+ !absX = abs x+ !absY = abs y+ !absZ = abs z+ !absW = abs w+ signum (Vector !x !y !z !w) = Vector signumX signumY signumZ signumW+ where+ !signumX = signum x+ !signumY = signum y+ !signumZ = signum z+ !signumW = signum w++ fromInteger x = Vector x' x' x' x'+ where+ !x' = fromInteger x++instance Fractional Vector where+ {-# SPECIALIZE INLINE (/) :: Vector -> Vector -> Vector #-}+ (Vector !(D# x) !(D# y) !(D# z) !(D# w)) / (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x /## x') (D# $ y /## y') (D# $ z /## z') (D# $ w /## w')+ fromRational x = Vector x' x' x' x'+ where+ !x' = fromRational x++instance Show Vector where+ show (Vector !x !y !z !w) = "(" ++ show x ++ ", " ++ show y ++ ", " ++ show z ++ ", " ++ show w ++ ")"++instance NFData Vector where+ rnf (Vector x y z w) = rnf x `seq` rnf y `seq` rnf z `seq` rnf w++tsTangent :: TangentSpace -> Normal+tsTangent (t, _, _) = t++tsBinormal :: TangentSpace -> Normal+tsBinormal (_, b, _) = b++tsNormal :: TangentSpace -> Normal+tsNormal = thr++xaxis :: Vector+xaxis = Vector 1 0 0 0++yaxis :: Vector+yaxis = Vector 0 1 0 0++zaxis :: Vector+zaxis = Vector 0 0 1 0++waxis :: Vector+waxis = Vector 0 0 0 1++zeroVector :: Vector+zeroVector = Vector 0 0 0 0++setWTo1 :: Vector -> Vector+{-# SPECIALIZE INLINE setWTo1 :: Vector -> Vector #-}+setWTo1 (Vector !x !y !z _) = Vector x y z 1++setWTo0 :: Vector -> Vector+{-# SPECIALIZE INLINE setWTo0 :: Vector -> Vector #-}+setWTo0 (Vector !x !y !z _) = Vector x y z 0++restoreOriginalW :: Vector -> Vector -> Vector+{-# SPECIALIZE INLINE restoreOriginalW :: Vector -> Vector -> Vector #-}+restoreOriginalW (Vector _ _ _ !w') (Vector !x !y !z _) = Vector x y z w'++madd :: Position -> Direction -> Double -> Vector+{-# SPECIALIZE INLINE madd :: Vector -> Vector -> Double -> Vector #-}+madd (Vector !(D# x) !(D# y) !(D# z) !(D# w)) (Vector !(D# x') !(D# y') !(D# z') !(D# w')) !(D# scalar) = Vector x'' y'' z'' w''+ where+ x'' = D# (x +## (x' *## scalar))+ y'' = D# (y +## (y' *## scalar))+ z'' = D# (z +## (z' *## scalar))+ w'' = D# (w +## (w' *## scalar))++negate :: Direction -> Direction+{-# SPECIALIZE INLINE Vector.negate :: Vector -> Vector #-}+negate (Vector !x !y !z !w) = Vector (-x) (-y) (-z) (-w)++vectorScalarMul :: Vector -> Double -> Vector+{-# SPECIALIZE INLINE vectorScalarMul :: Vector -> Double -> Vector #-}+(Vector !x !y !z !w) `vectorScalarMul` k = Vector (x * k) (y * k) (z * k) (w * k)++infixl 7 <*>+infixl 7 </>+--infixl 6 <+>+--infixl 6 <->++(</>) :: Vector -> Double -> Vector+{-# SPECIALIZE INLINE (</>) :: Vector -> Double -> Vector #-}+a </> b = a `vectorScalarMul` (1 / b)++(<*>) :: Vector -> Double -> Vector+{-# SPECIALIZE INLINE (<*>) :: Vector -> Double -> Vector #-}+a <*> b = a `vectorScalarMul` b++dot3 :: Vector -> Vector -> Double+{-# SPECIALIZE INLINE dot3 :: Vector -> Vector -> Double #-}+(Vector !(D# x) !(D# y) !(D# z) _) `dot3` (Vector !(D# x') !(D# y') !(D# z') _) = D# $ (x *## x') +## (y *## y') +## (z *## z')++dot4 :: Vector -> Vector -> Double+{-# SPECIALIZE INLINE dot4 :: Vector -> Vector -> Double #-}+(Vector !(D# x) !(D# y) !(D# z) !(D# w)) `dot4` (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = D# $ (x *## x') +## (y *## y') +## (z *## z') +## (w *## w')++sdot3 :: Vector -> Vector -> Double+{-# SPECIALIZE INLINE dot3 :: Vector -> Vector -> Double #-}+(Vector !(D# x) !(D# y) !(D# z) _) `sdot3` (Vector !(D# x') !(D# y') !(D# z') _) = D# $ saturate## ((x *## x') +## (y *## y') +## (z *## z'))++sdot4 :: Vector -> Vector -> Double+{-# SPECIALIZE INLINE sdot4 :: Vector -> Vector -> Double #-}+(Vector !(D# x) !(D# y) !(D# z) !(D# w)) `sdot4` (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = D# $ saturate## ((x *## x') +## (y *## y') +## (z *## z') +## (w *## w'))++cross :: Direction -> Direction -> Direction+{-# SPECIALIZE INLINE cross :: Vector -> Vector -> Vector #-}+(Vector !(D# x1) !(D# y1) !(D# z1) _) `cross` (Vector !(D# x2) !(D# y2) !(D# z2) _) = Vector x y z 0+ where+ !x = D# ((y1 *## z2) -## (y2 *## z1))+ !y = D# ((z1 *## x2) -## (z2 *## x1))+ !z = D# ((x1 *## y2) -## (x2 *## y1))++magnitude :: Vector -> Double+{-# SPECIALIZE INLINE magnitude :: Vector -> Double #-}+magnitude !vec = sqrt (magnitudeSq vec)++magnitudeSq :: Vector -> Double+{-# SPECIALIZE INLINE magnitudeSq :: Vector -> Double #-}+magnitudeSq (Vector !(D# x#) !(D# y#) !(D# z#) _) = D# ((x# *## x#) +## (y# *## y#) +## (z# *## z#))++normalise :: Direction -> Direction+{-# SPECIALIZE INLINE normalise :: Direction -> Direction #-}+normalise !a = setWTo0 (a `vectorScalarMul` (1 / magnitude a))++distance :: Position -> Position -> Double+{-# SPECIALIZE INLINE distance :: Position -> Position -> Double #-}+distance !a !b = magnitude (a - b)++distanceSq :: Position -> Position -> Double+{-# SPECIALIZE INLINE distanceSq :: Position -> Position -> Double #-}+distanceSq !a !b = magnitudeSq (a - b)++reflect :: Direction -> Direction -> Direction+{-# SPECIALIZE INLINE reflect :: Direction -> Direction -> Direction #-}+reflect !incoming !normal = restoreOriginalW incoming $ (normal `vectorScalarMul` (2 * (normal `dot3` incoming))) - incoming++refract :: Direction -> Direction -> Double -> Direction+{-# SPECIALIZE INLINE refract :: Vector -> Vector -> Double -> Vector #-}+refract !incoming !normal !eta+ | cosTheta1 >## 0.0## = setWTo0 $ (l `vectorScalarMul` eta) + (normal `vectorScalarMul` D# (eta# *## cosTheta1 -## cosTheta2))+ | otherwise = setWTo0 $ (l `vectorScalarMul` eta) + (normal `vectorScalarMul` D# (eta# *## cosTheta1 +## cosTheta2))+ where !(D# cosTheta1) = normal `dot3` incoming+ !cosTheta2 = sqrtDouble# (1.0## -## eta# **## 2.0## *## (1.0## -## cosTheta1 **## 2.0##))+ !l = Vector.negate incoming+ !(D# eta#) = eta++largestAxis :: Vector -> Int+largestAxis (Vector !x !y !z _) + | abs x >= abs y && abs x >= abs z = 0+ | abs y >= abs x && abs y >= abs z = 1+ | abs z >= abs x && abs z >= abs y = 2+ | otherwise = error "largestAxis: Undefined case"++nthLargestAxis :: Vector -> Int -> Int+nthLargestAxis (Vector !x !y !z _) order + | order < 3 = snd (sort [(abs x, 0), (abs y, 1), (abs z, 2)] !! order)+ | otherwise = error "nthLargestAXis: Undefined case"++min :: Vector -> Vector -> Vector+{-# SPECIALIZE INLINE Vector.min :: Vector -> Vector -> Vector #-}+min (Vector !x1 !y1 !z1 !w1) (Vector !x2 !y2 !z2 !w2) = Vector x y z w+ where+ !x = Prelude.min x1 x2+ !y = Prelude.min y1 y2+ !z = Prelude.min z1 z2+ !w = Prelude.min w1 w2++max :: Vector -> Vector -> Vector+{-# SPECIALIZE INLINE Vector.max :: Vector -> Vector -> Vector #-}+max (Vector !x1 !y1 !z1 !w1) (Vector !x2 !y2 !z2 !w2) = Vector x y z w+ where+ !x = Prelude.max x1 x2+ !y = Prelude.max y1 y2+ !z = Prelude.max z1 z2+ !w = Prelude.max w1 w2++directionToSpherical :: Direction -> (Double, Double)+directionToSpherical (Vector !x !y !z _) = (theta, phi)+ where+ theta = acos z / pi+ phi = (atan2 y x + pi) / (2.0 * pi)++sphericalToDirection :: Double -> Double -> Direction+sphericalToDirection (D# !theta) (D# !phi) = Vector (D# $ sinDouble# theta *## cosDouble# phi) (D# $ sinDouble# theta *## sinDouble# phi) (D# $ cosDouble# theta) 1++component :: Vector -> Int -> Double+{-# SPECIALIZE INLINE component :: Vector -> Int -> Double #-}+component (Vector !x _ _ _) 0 = x+component (Vector _ !y _ _) 1 = y+component (Vector _ _ !z _) 2 = z+component (Vector _ _ _ !w) 3 = w+component _ _ = error "Invalid component index"++transformDir :: Direction -> TangentSpace -> Direction+{-# SPECIALIZE INLINE transformDir :: Direction -> TangentSpace -> Direction #-}+transformDir (Vector !x !y !z _) !(tangent, binormal, normal) = setWTo0 ((tangent <*> x) + (binormal <*> y) + (normal <*> z))
+ app/src/scripts/build view
@@ -0,0 +1,4 @@+#!/bin/sh+ghc -tmpdir tmp -hidir hi -odir obj -fext-core -fexcess-precision -funbox-strict-fields -threaded -rtsopts -fwarn-missing-signatures -Wall -O2 -fspec-constr -fliberate-case -fstatic-argument-transformation Main.hs -o ../crocodile $1 $2 $3 $4 $5 $6 $7 $8 $9+mv *.hcr core/+scripts/unit_test
+ app/src/scripts/clean view
@@ -0,0 +1,2 @@+#!/bin/sh+rm obj/*.o
+ app/src/scripts/run-hlint view
@@ -0,0 +1,3 @@+#!/bin/sh+~/.cabal/bin/hlint *.hs > hlint-output.txt+less hlint-output.txt
+ app/src/scripts/unit_test view
@@ -0,0 +1,2 @@+#!/bin/sh+runhaskell Tests/UnitTests.hs
+ crocodile.cabal view
@@ -0,0 +1,58 @@++-- Crocodile.cabal auto-generated by cabal init. For additional+-- options, see+-- http://www.haskell.org/cabal/release/cabal-latest/doc/users-guide/authors.html#pkg-descr.+-- The name of the package.+Name: crocodile++-- The package version. See the Haskell package versioning policy+-- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for+-- standards guiding when and how versions should be incremented.+Version: 0.1++-- A short (one-line) description of the package.+Synopsis: An offline renderer supporting ray tracing and photon mapping++-- A longer description of the package.+Description: This is an offline renderer written purely in Haskell, supporting ray tracing and photon mapping with an irradiance cache++-- URL for the project homepage or repository.+Homepage: https://github.com/TomHammersley/HaskellRenderer/++-- The license under which the package is released.+License: GPL-2++-- The file containing the license text.+License-file: LICENSE++-- The package author(s).+Author: Tom Hammersley++-- An email address to which users can send suggestions, bug reports,+-- and patches.+Maintainer: tomhammersley@gmail.com++-- A copyright notice.+-- Copyright: ++Category: Graphics++Build-type: Simple++-- Extra files to be distributed with the package, such as examples or+-- a README.+Extra-source-files: README LICENSE Setup.hs app/src/*.hs app/src/*.hs-boot app/src/Tests/*.hs app/src/scripts/build app/src/scripts/clean app/src/scripts/run-hlint app/src/scripts/unit_test++-- Constraint on the version of Cabal needed to build this package.+Cabal-version: >=1.6++Executable crocodile+ -- Packages needed in order to build this package.+ Build-depends: base >= 4 && < 5, HUnit, heap, bmp, mtl, mersenne-random-pure64, ghc-prim, parallel, bytestring, deepseq++ -- .hs or .lhs file containing the Main module.+ Main-is: Main.hs+ + Hs-Source-Dirs: app/src+ + ghc-options: -O2 -fext-core -fexcess-precision -funfolding-use-threshold=16 -funbox-strict-fields -funfolding-keeness-factor=10 -threaded -rtsopts -fwarn-missing-signatures -Wall -feager-blackholing