diff --git a/COPYING b/COPYING new file mode 100644 --- /dev/null +++ b/COPYING @@ -0,0 +1,340 @@ + 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. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 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. + + , 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. + diff --git a/Clr.hs b/Clr.hs new file mode 100644 --- /dev/null +++ b/Clr.hs @@ -0,0 +1,41 @@ +module Clr where + +type CFlt = Double +data Color = Color {r,g,b :: !CFlt} deriving Show + +c_black = Color 0 0 0 +c_white = Color 1 1 1 +c_red = Color 1 0 0 +c_green = Color 0 1 0 +c_blue = Color 0 0 1 + +cadd :: Color -> Color -> Color +cadd (Color r1 g1 b1) (Color r2 g2 b2) = + Color (r1+r2) (g1+g2) (b1+b2) + +cdiv :: Color -> CFlt -> Color +cdiv c1 div = + cscale c1 (1/div) + +cscale :: Color -> CFlt -> Color +cscale (Color r g b) mul = + Color (r * mul) + (g * mul) + (b * mul) + +cmul :: Color -> Color -> Color +cmul (Color r1 g1 b1) (Color r2 g2 b2) = + Color (r1*r2) (g1*g2) (b1*b2) + +cavg :: Color -> Color -> Color +cavg c1 c2 = cscale (cadd c1 c2) 0.5 + +cscaleadd :: Color -> Color -> CFlt -> Color +cscaleadd (Color r1 g1 b1) (Color r2 g2 b2) mul = + Color (r1+(r2*mul)) (g1+(g2*mul)) (b1+(b2*mul)) + +cclamp :: Color -> Color +cclamp (Color r g b) = + Color (if r > 0.0 then r else 0.0) + (if g > 0.0 then g else 0.0) + (if b > 0.0 then b else 0.0) diff --git a/Glome.hs b/Glome.hs new file mode 100644 --- /dev/null +++ b/Glome.hs @@ -0,0 +1,205 @@ +import Vec +import Clr +import Solid +import Trace +import Spd +import TestScene +import Graphics.Rendering.OpenGL +import Graphics.UI.GLUT as GLUT +import System.CPUTime +import Control.Parallel.Strategies +import Data.Time.Clock.POSIX +import System.Exit +import IO +import System +import System.Console.GetOpt +import Data.Maybe( fromMaybe ) + +-- import Debug.Trace +-- import Data.ByteString + +get_color :: Flt -> Flt -> Scene -> Clr.Color +-- +get_color x y scn = + let (Scene sld lights (Camera pos fwd up right) dtex bgcolor) = scn + dir = vnorm $ vadd3 fwd (vscale right (-x)) (vscale up y) + ray = (Ray pos dir) + in + Trace.trace scn ray infinity 2 +--} + +{-- +-- for testing screen-draw overhead +get_color x y scn = + Clr.Color (1*x) (1*y) 0 +--} + +fc :: Flt -> Float +--fc x = if x == x then clamp 0 x 0.5 else error "nan" +fc x = realToFrac x +--fc x = x+delta + +-- double nested for loop +-- for rendering a sub-box of the screen +gen_pixels curx cury stopx stopy maxx maxy scene = + let midx = maxx/2 + midy = maxy/2 + gp x y = + if y>=stopy then + return () + else + if x>=stopx + then + gp curx (y+1) + else + do + let scx = (x-midx) / midx + let scy = (y-midy) / midy + -- let (Clr.Color r g b) = get_color scx (scy*(midy/midx)) scene + let (Clr.Color r g b) = get_color (scx*(midx/midy)) scy scene + currentColor $= Color4 (fc r) (fc g) (fc b) 1 + vertex$Vertex3 scx scy 0 + gp (x+1) y + in gp curx cury + +-- same as above, but non-monadic, instead return list of pixels +gen_pixel_list :: Flt -> Flt -> Flt -> Flt -> Flt -> Flt -> Scene -> [(Flt,Flt,Flt,Flt,Flt)] +gen_pixel_list curx cury stopx stopy maxx maxy scene = + let midx = maxx/2 + midy = maxy/2 + gp x y = + if y>=stopy then + [] + else + if x>=stopx + then + gp curx (y+1) + else + let scx = (x-midx) / midx + scy = (y-midy) / midy + (Clr.Color r g b) = get_color scx (scy*(midy/midx)) scene + --(Clr.Color r g b) = get_color (scx*(midx/midy)) scy scene + in + (scx,scy,r,g,b) : (gp (x+1) y) + in gp curx cury + +-- split screen into little blocks for parallel rendering +{- +gen_blocks maxx maxy block_size scene = + let foo = 1 + bar = 2 + gb x y = + if y>=maxy then + return () + else + if x>=maxx then + gb 0 (y+block_size) + else + do + gen_pixels x y (x+block_size-1) (y+block_size-1) maxx maxy scene + gb (x+block_size) y + in gb 0 0 +-} + +{- +-- this doesn't seem to parallelize +gen_blocks maxx maxy block_size scene = + let xblocks = maxx/block_size + yblocks = maxy/block_size + blocks = Prelude.map (\x -> Prelude.map (\y -> (x*block_size,y*block_size) ) [0..yblocks-1] ) [0..xblocks-1] + in + do + -- mapM_ (\(x,y) -> gen_pixels x y (x+block_size) (y+block_size) maxx maxy scene) (concat blocks) + -- sequence_ $ Prelude.map (\(x,y) -> gen_pixels x y (x+block_size) (y+block_size) maxx maxy scene) (concat blocks) + sequence_ $ (parMap rwhnf) + (\(x,y) -> gen_pixels x y (x+block_size) (y+block_size) maxx maxy scene) + (concat blocks) +-} + +-- as above, but parallel code is pure +-- odd, this seems to be faster, even without multicore +-- rnf is faster than rwhnf +gen_blocks_list maxx maxy block_size scene = + let xblocks = maxx/block_size + yblocks = maxy/block_size + blocks = Prelude.concat $ Prelude.map (\x -> Prelude.map (\y -> (x*block_size,y*block_size) ) [0..yblocks-1] ) [0..xblocks-1] + pixels = map -- (parMap rnf) + (\(x,y) -> gen_pixel_list x y (x+block_size) (y+block_size) maxx maxy scene) + (blocks) + in + do + mapM_ (\pix -> mapM_ (\(x,y,r,g,b) -> do currentColor $= Color4 (fc r) (fc g) (fc b) 1 + vertex$Vertex3 (fc x) (fc y) 0 + ) pix) pixels + + + +-- Haskell opengl tutorial: +-- http://blog.mikael.johanssons.org/archive/2006/09/ +-- opengl-programming-in-haskell-a-tutorial-part-1/ + +-- another tutorial: +-- http://www.tfh-berlin.de/~panitz/hopengl/skript.html + +main :: IO () +main = + do + -- parse arguments + -- http://leiffrenzel.de/papers/commandline-options-in-haskell.html + args <- getArgs + let (flags, nonOpts, msgs) = getOpt RequireOrder options args + print $ "recognized options: " ++ (show (length flags)) + scene <- getscene flags + let sx = 720 :: GLsizei + let sy = 480 :: GLsizei + (name, _) <- getArgsAndInitialize + initialDisplayMode $= [] +-- initialDisplayMode $= [DoubleBuffered] + createWindow name + windowSize $= Size sx sy + displayCallback $= display scene sx sy + keyboardMouseCallback $= Just (keyboard scene) + mainLoop + +display scene sx sy = do + t1 <- getPOSIXTime + clearColor $= Color4 0 0 0 1 + clear [ColorBuffer] + --(Size sx sy) <- GLUT.get windowSize + let sizex = fromIntegral sx + let sizey = fromIntegral sy + -- renderPrimitive Points $ gen_blocks_list 512 512 128 scene + -- renderPrimitive Points $ gen_blocks_list 720 480 80 scene + renderPrimitive Points $ gen_pixels 0 0 sizex sizey sizex sizey scene + swapBuffers + t2 <- getPOSIXTime + -- print ((fromInteger (t2-t1))/1000000000) + print (t2-t1) + +keyboard _ (Char 'q') Down _ _ = + do + exitWith ExitSuccess + +keyboard s (Char 's') Down _ _ = + do + print (show s) + +keyboard _ _ _ _ _ = return () + + +getscene :: [Flag] -> IO Scene +getscene flags = + case flags of + [] -> TestScene.scn + (Filename s:xs) -> do filedes <- openFile s ReadMode + filestring <- (IO.hGetContents filedes) + (scene,s) <- return $ Prelude.head $ reads filestring + return scene + + +data Flag = Filename String | Res Int Int deriving Show + +options :: [OptDescr Flag] +options = + [ Option ['n'] ["filename"] (ReqArg Filename "FILE") "input NFF scene" + ] diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,14 @@ + This program, Glome.hs, is copyright 2008 Jim Snow + + This program is free software; you can redistribute it and/or modify + it under the terms of version 2 of the GNU General Public License as + published by the Free Software Foundation; + + 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 St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,39 @@ +Glome.hs is a haskell port of my ocaml raytracer, glome. +(http://syn.cs.pdx.edu/~jsnow/glome) + +To compile on unix-type systems, execute "./make", to run, execute "./run". +Otherwise, invoke compiler commands manually. + +Update: glome has been converted over to cabal, so you can now invoke + +> runhaskell Setup.lhs configure --prefix=$HOME --user +> runhaskell Setup.lhs build +> runhaskell Setup.lhs install + +(The "make" and "run" scripts should still work.) + +Glome.hs depends on opengl bindings, which come in the standard ghc distribution. +Glome has been tested with ghc 6.8.2. + +Features: +- can load files in NFF format +- handles diffuse illumination, shadows, and reflection +- renders triangles, spheres, cylinders, cones, disks, boxes + and planes as base primitives +- supports csg group, difference, and intersection primitives +- supports transformations (rotate, scale, translate) of + arbitrary geometry +- perlin noise and a few basic textures are implemented +- uses a bounding interval heirarchy acceleration structure + or you can construct a BVH manually with the "Bound" primitive +- multiprocessor support is currently disabled + +Using: to load an NFF scene, run "./Glome -n [filename]". +Otherwise, a default scene is rendered, defined in "TestScene.hs". + +Many of the features are only accessible by programming directly in +Haskell via TestScene.hs, as NFF does not support CSG or textures. + +Refraction and photon mapping are not yet implemented. + +-jim diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,4 @@ +#! /usr/bin/env runhaskell + +> import Distribution.Simple +> main = defaultMain diff --git a/Solid.hs b/Solid.hs new file mode 100644 --- /dev/null +++ b/Solid.hs @@ -0,0 +1,1183 @@ +module Solid where +import Vec +import Clr +import Data.List hiding (group) + +--COMMON DATATYPES AND FUNCTIONS-- +data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show +data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple + +--union of two bounding boxes +bbjoin :: Bbox -> Bbox -> Bbox +bbjoin (Bbox p1a p2a) (Bbox p1b p2b) = + (Bbox (vmin p1a p1b) (vmax p2a p2b)) + +--overlap of two bounding boxes +bboverlap :: Bbox -> Bbox -> Bbox +bboverlap (Bbox p1a p2a) (Bbox p1b p2b) = + (Bbox (vmax p1a p1b) (vmin p2a p2b)) + +--split a bounding box into two +bbsplit :: Bbox -> Int -> Flt -> (Bbox,Bbox) +bbsplit (Bbox p1 p2) axis offset = + if (offset < (va p1 axis)) || (offset > (va p2 axis)) + then error "degenerate bounding box split" + else ((Bbox p1 (vset p2 axis offset)), + (Bbox (vset p1 axis offset) p2)) + +-- generate a bounding box from a list of points +bbpts :: [Vec] -> Bbox +bbpts [] = empty_bbox +bbpts ((Vec x y z):[]) = + Bbox (Vec (x-delta) (y-delta) (z-delta)) + (Vec (x+delta) (y+delta) (z+delta)) + +bbpts ((Vec x y z):pts) = + let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bbpts pts + minx = fmin (x-delta) p1x + miny = fmin (y-delta) p1y + minz = fmin (z-delta) p1z + maxx = fmax (x+delta) p2x + maxy = fmax (y+delta) p2y + maxz = fmax (z+delta) p2z in + Bbox (Vec minx miny minz) (Vec maxx maxy maxz) + +-- surface area, volume +bbsa :: Bbox -> Flt +bbsa (Bbox p1 p2) = + let Vec dx dy dz = vsub p2 p1 + in dx*dy + dx*dz + dy*dz + +bbvol :: Bbox -> Flt +bbvol (Bbox p1 p2) = + let (Vec dx dy dz) = vsub p2 p1 + in dx*dy*dz + +empty_bbox = + Bbox (Vec infinity infinity infinity) + (Vec (-infinity) (-infinity) (-infinity)) + +everything_bbox = + Bbox (Vec (-infinity) (-infinity) (-infinity)) + (Vec infinity infinity infinity) + + +bbclip :: Ray -> Bbox -> Interval +bbclip (Ray (Vec ox oy oz) (Vec dx dy dz)) + (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = + let dxrcp = 1/dx + dyrcp = 1/dy + dzrcp = 1/dz + Interval inx outx = if dx > 0 + then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp) + else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp) + Interval iny outy = if dy > 0 + then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp) + else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp) + Interval inz outz = if dz > 0 + then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp) + else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp) + in + Interval (fmax3 inx iny inz) (fmin3 outx outy outz) + +data Rayint = RayHit { + depth :: !Flt, + pos :: !Vec, + norm :: !Vec, + texture :: !Texture +} | RayMiss deriving Show + +nearest :: Rayint -> Rayint -> Rayint +nearest a RayMiss = a +nearest RayMiss b = b +nearest (RayHit da pa na ta) (RayHit db pb nb tb) = + if da < db + then RayHit da pa na ta + else RayHit db pb nb tb + +furthest :: Rayint -> Rayint -> Rayint +furthest a RayMiss = RayMiss +furthest RayMiss b = RayMiss +furthest (RayHit da pa na ta) (RayHit db pb nb tb) = + if da > db + then RayHit da pa na ta + else RayHit db pb nb tb + +hit :: Rayint -> Bool +hit (RayHit _ _ _ _) = True +hit RayMiss = False + +dist :: Rayint -> Flt +dist RayMiss = infinity +dist (RayHit d _ _ _) = d + + +--LIGHTS-- +data Light = Light {litpos :: !Vec, + litcol :: !Color} deriving Show + +--MATERIALS-- +data Material = Material {clr :: Color, + reflect, refract, ior, + kd, shine :: !Flt} deriving Show +type Texture = Rayint -> Material + +-- this is sort of a no-op; we don't have a +-- good way to show an arbitrary function +showTexture :: Texture -> String +showTexture t = show $ t RayMiss + +instance Show Texture where + show = showTexture + +m_white = (Material c_white 0 0 0 1 2) +t_white ri = m_white + +t_uniform :: Material -> Texture +t_uniform m = \x -> m + +interp :: Flt -> Flt -> Flt -> Flt +interp scale a b = + scale*a + (1-scale)*b + +--not really correct, but we'll go with it for now +m_interp :: Material -> Material -> Flt -> Material +m_interp m1 m2 scale = + let (Material m1c m1refl m1refr m1ior m1kd m1shine) = m1 + (Material m2c m2refl m2refr m2ior m2kd m2shine) = m2 + intp = interp scale + c = cadd (cscale m1c scale) (cscale m2c (1-scale)) + refl = intp m1refl m2refl + refr = intp m1refr m2refr + ior = intp m1ior m2ior + kd = intp m1kd m2kd + shine = intp m1shine m2shine + in (Material c refl refr ior kd shine) + +--SOLID TYPES-- + +data Solid = Sphere {center :: !Vec, + radius, invradius :: !Flt} + | Triangle {v1, v2, v3 :: Vec} + | TriangleNorm {v1, v2, v3, n1, n2, n3 :: Vec} + | Disc !Vec !Vec !Flt -- position, normal, r*r + | Cylinder !Flt !Flt !Flt -- radius height1 height2 + | Cone !Flt !Flt !Flt !Flt -- r clip1 clip2 height + | Plane Vec Flt -- normal, offset from origin + | Box !Bbox + | Group ![Solid] + | Intersection ![Solid] + | Bevel !Solid !Flt + | Bound Solid Solid + | Difference !Solid !Solid + | Bih {bihbb :: !Bbox, bihroot :: !BihNode} + | Instance !Solid !Xfm + | Tex !Solid Texture + | Portal Solid Solid -- if we hit a, intersect with b + | Photon Vec Vec Color Flt -- pos, incident ray, color, radius + | Nothing deriving Show -- conflicts with + -- Nothing :: Maybe a from prelude + +data BihNode = BihLeaf !Solid + | BihBranch {lmax :: !Flt, rmin :: !Flt, ax :: !Int, + l :: BihNode, r :: BihNode} deriving Show + +--CONSTRUCTORS-- +sphere :: Vec -> Flt -> Solid +sphere c r = + Sphere c r (1.0/r) + +triangle :: Vec -> Vec -> Vec -> Solid +triangle v1 v2 v3 = + Triangle v1 v2 v3 + +--simple tesselation +triangles :: [Vec] -> [Solid] +triangles (v1:vs) = + zipWith (\v2 v3 -> triangle v1 v2 v3) vs (tail vs) + +trianglenorm v1 v2 v3 n1 n2 n3 = + -- Triangle v1 v2 v3 + TriangleNorm v1 v2 v3 n1 n2 n3 + +trianglesnorms :: [(Vec,Vec)] -> [Solid] +trianglesnorms (vn1:vns) = + zipWith (\vn2 vn3 -> trianglenorm (fst vn1) (fst vn2) (fst vn3) + (snd vn1) (snd vn2) (snd vn3)) + vns (tail vns) + +box :: Vec -> Vec -> Solid +box p1 p2 = + Box (Bbox p1 p2) + +disc :: Vec -> Vec -> Flt -> Solid +disc pos norm r = + Disc pos norm (r*r) + +cylinder_z :: Flt -> Flt -> Flt -> Solid +cylinder_z r h1 h2 = Cylinder r h1 h2 + +cone_z :: Flt -> Flt -> Flt -> Flt -> Solid +cone_z r h1 h2 height = Cone r h1 h2 height + +-- construct a general cylinder from p1 to p2 with radius r +cylinder :: Vec -> Vec -> Flt -> Solid +cylinder p1 p2 r = + let axis = vsub p2 p1 + len = vlen axis + ax1 = vscale axis (1/len) + (ax2,ax3) = orth ax1 + in Instance (cylinder_z r 0 len) + (compose [ (xyz_to_uvw ax2 ax3 ax1), + (translate p1) ]) + +-- similar for cone +cone :: Vec -> Flt -> Vec -> Flt -> Solid +cone p1 r1 p2 r2 = + if r1 < r2 + then cone p2 r2 p1 r1 + else if r1-r2 < delta + then cylinder p1 p2 r2 + else + let axis = vsub p2 p1 + len = vlen axis + ax1 = vscale axis (1/len) + (ax2,ax3) = orth ax1 + height = (r1*len)/(r1-r2) -- distance to end point + in + Instance (cone_z r1 0 len height) + (compose [ (xyz_to_uvw ax2 ax3 ax1), + (translate p1) ]) + +plane :: Vec -> Vec -> Solid +plane orig norm_ = Plane norm d + where norm = vnorm norm_ + d = vdot orig norm + +-- flatten tree of groups into a single group +flatten_group :: [Solid] -> [Solid] +flatten_group ((Group slds):xs) = + (flatten_group slds) ++ xs +flatten_group ((Solid.Nothing):xs) = xs +flatten_group x = x + +group :: [Solid] -> Solid +group [] = Solid.Nothing +group (sld:[]) = sld +group slds = + Group (flatten_group slds) + +transform :: Solid -> [Xfm] -> Solid +transform (Instance s xfm2) xfm1 = + transform s [compose ([xfm2] ++ xfm1)] + +transform s xfm = + Instance s (compose xfm) + +-- push all the transforms out to the leaves +-- and throw away pre-existing bounding volumes +-- so we can run the bih constructor on the +-- resulting group +flatten_transform :: Solid -> [Solid] +flatten_transform (Group slds) = + flatten_group $ concat (map flatten_transform slds) + +flatten_transform (Instance s xfm) = + case s of + Group slds -> flatten_transform $ group (map (\x -> transform x [xfm]) slds) + Bound sa sb -> flatten_transform (transform sb [xfm]) + Instance sa xfm2 -> flatten_transform (transform s [xfm]) + _ -> [transform s [xfm]] + +flatten_transform (Bound sa sb) = flatten_transform sb + +-- bih construction +build_leaf objs = + BihLeaf (group (map snd objs)) + +max_bih_sa = 0.3 :: Flt + +build_rec :: [(Bbox,Solid)] -> Bbox -> Bbox -> Int -> BihNode +build_rec objs nodebox splitbox depth = + -- if (null objs) || (null $ tail objs) || + -- (null $ tail $ tail objs) + if length objs < 2 + then build_leaf objs + else + let (Bbox nodeboxp1 nodeboxp2) = nodebox + (Bbox splitboxp1 splitboxp2) = splitbox + axis = vmaxaxis (vsub splitboxp2 splitboxp1) + bbmin = va splitboxp1 axis + bbmax = va splitboxp2 axis + candidate = (bbmin + bbmax) * 0.5 + in + if candidate > (va nodeboxp2 axis) then + build_rec objs nodebox + (Bbox splitboxp1 (vset splitboxp2 axis candidate)) + depth + else + if candidate < (va nodeboxp1 axis) then + build_rec objs nodebox ( + Bbox (vset splitboxp1 axis candidate) splitboxp2) + depth + else + -- not sure if this is a big win + let nbsa = bbsa nodebox + (big,small) = partition (\ (bb,_) -> + (bbsa bb) > (nbsa * max_bih_sa)) objs + in + if (not $ null big) && ((length big) < ((length small)*2)) + then (BihBranch (va nodeboxp2 0) (va nodeboxp1 0) 0 + (build_rec big nodebox splitbox (depth+1)) + (build_rec small nodebox splitbox (depth+1)) ) + else + let (l,r) = partition (\((Bbox bbp1 bbp2),_)-> + (((va bbp1 axis)+(va bbp2 axis))*0.5) + < candidate ) objs + lmax = foldl fmax (-infinity) (map (\((Bbox _ p2),_) -> va p2 axis) l) + rmin = foldl fmin infinity (map (\((Bbox p1 _),_) -> va p1 axis) r) + (lsplit,rsplit) = bbsplit splitbox axis candidate + lnb = (Bbox nodeboxp1 (vset nodeboxp2 axis lmax)) + rnb = (Bbox (vset nodeboxp1 axis rmin) nodeboxp2) + in + -- stop if there's no progress being made + if ((null l) && (rmin <= bbmin)) || + ((null r) && (lmax >= bbmax)) + then build_leaf objs + else + (BihBranch (lmax+delta) (rmin-delta) axis + (build_rec l lnb lsplit (depth+1)) + (build_rec r rnb rsplit (depth+1)) ) + +bih :: [Solid] -> Solid +bih [] = Solid.Nothing +-- bih (sld:[]) = sld -- sometimes we'd like to be able to use a + -- single object bih just for its aabb +bih slds = + let objs = map (\x -> ((bound x),x)) (flatten_group slds) + bb = foldl bbjoin empty_bbox (map (\(b,_)->b) objs) + root = build_rec objs bb bb 0 + (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bb + in + if p1x == (-infinity) || p1y == (-infinity) || p1z == (-infinity) || + p2x == infinity || p2y == infinity || p2z == infinity + then + error $ "bih: infinite bounding box " ++ (show objs) + else + (Bih bb root) + +--INTERSECTION TESTS-- +rayint :: Solid -> Ray -> Flt -> Texture -> Rayint + +--Basic Primitives-- +--Triangle-- +-- adaptation of Moller and Trumbore from pbrt page 127 +rayint (Triangle p1 p2 p3) (Ray o dir) dist tex = + let e1 = vsub p2 p1 + e2 = vsub p3 p1 + s1 = vcross dir e2 + divisor = vdot s1 e1 + in + if (divisor == 0) + then RayMiss + else + let invdivisor = 1.0 / divisor + d = vsub o p1 + b1 = (vdot d s1) * invdivisor + in + if (b1 < 0) || (b1 > 1) + then RayMiss + else + let s2 = vcross d e1 + b2 = (vdot dir s2) * invdivisor + in + if (b2 < 0) || (b1 + b2 > 1) + then RayMiss + else + let t = (vdot e2 s2) * invdivisor + in + if (t < 0) || (t > dist) + then RayMiss + else + RayHit t (vscaleadd o dir t) (vnorm (vcross e1 e2)) tex + +rayint (TriangleNorm p1 p2 p3 n1 n2 n3) (Ray o dir) dist tex = + let e1 = vsub p2 p1 + e2 = vsub p3 p1 + s1 = vcross dir e2 + divisor = vdot s1 e1 + in + if (divisor == 0) + then RayMiss + else + let invdivisor = 1.0 / divisor + d = vsub o p1 + b1 = (vdot d s1) * invdivisor + in + if (b1 < 0) || (b1 > 1) + then RayMiss + else + let s2 = vcross d e1 + b2 = (vdot dir s2) * invdivisor + in + if (b2 < 0) || (b1 + b2 > 1) + then RayMiss + else + let t = (vdot e2 s2) * invdivisor + in + if (t < 0) || (t > dist) + then RayMiss + else + let n1scaled = (vscale n1 (1-(b1+b2))) + n2scaled = (vscale n2 b1) + n3scaled = (vscale n3 b2) + norm = vnorm $ vadd3 n1scaled n2scaled n3scaled + in RayHit t (vscaleadd o dir t) norm tex + +--Sphere-- +-- adapted from graphics gems volume 1 +rayint (Sphere center r invr) (Ray e dir) dist t = + let eo = vsub center e + v = vdot eo dir + in + if (dist >= (v - r)) && (v > 0.0) + then + let vsqr = v*v + csqr = vdot eo eo + rsqr = r*r + disc = rsqr - (csqr - vsqr) in + if disc < 0.0 then + RayMiss + else + let d = sqrt disc + hitdist = if (v-d) > 0 then (v-d) else (v+d) + in if (hitdist < 0) || (hitdist > dist) + then RayMiss + else + let p = vscaleadd e dir hitdist + -- n = vscale (vsub p center) invr in + -- n = vsub (vscale p invr) (vscale center invr) in + n = vnorm (vsub p center) + in RayHit hitdist p n t + else + RayMiss + +-- nice and simple +rayint (Disc point norm radius_sqr) r d t = + let (Ray orig dir) = r + dist = plane_int_dist r point norm + in if dist < 0 || dist > d + then RayMiss + else let pos = vscaleadd orig dir dist + offset = vsub pos point + in + if (vdot offset offset) > radius_sqr + then RayMiss + else RayHit dist pos norm t + +-- cylinder aligned to z axis +-- no end caps +-- adapted from pbrt +rayint (Cylinder r h1 h2) (Ray orig dir) d t = + let Vec ox oy oz = orig + Vec dx dy dz = dir + a = dx*dx + dy*dy + b = 2*(dx*ox + dy*oy) + c = ox*ox + oy*oy - r*r + disc = b*b - 4*a*c + in if disc < 0 + then RayMiss + else + let discsqrt = sqrt disc + q = if b < 0 + then (b-discsqrt)*(-0.5) + else (b+discsqrt)*(-0.5) + t0' = q/a + t1' = c/q + t0 = fmin t0' t1' + t1 = fmax t0' t1' + in if t1 < 0 || t0 > d + then RayMiss + else let dist = if t0 < 0 + then t1 + else t0 + in if dist < 0 || dist > d + then RayMiss + else let pos = vscaleadd orig dir dist + Vec posx posy posz = pos + in if posz > h1 && posz < h2 + then RayHit dist pos (Vec (posx/r) (posy/r) 0) t + else if dz > 0 -- ray pointing up from bottom + then if oz < h1 + then rayint (Disc (Vec 0 0 h1) nvz (r*r)) (Ray orig dir) d t + --then rayint_aadisc h1 r (Ray orig dir) d t + else RayMiss + else if oz > h2 + then rayint (Disc (Vec 0 0 h2) vz (r*r)) (Ray orig dir) d t + --rayint_aadisc h2 r (Ray orig dir) d t -- todo: fix normal + else RayMiss + +-- cone centered on z axiz, height of hp, clipped at h1 and h2 +rayint (Cone r clip1 clip2 height) (Ray orig dir) d t = + let Vec ox oy oz = orig + Vec dx dy dz = dir + k' = (r/height) + k = k'*k' + a = dx*dx + dy*dy - k*dz*dz + b = 2*(dx*ox + dy*oy - k*dz*(oz-height)) + c = ox*ox + oy*oy - k*(oz-height)*(oz-height) + disc = b*b - 4*a*c + in if disc < 0 + then RayMiss + else + let discsqrt = sqrt disc + q = if b < 0 + then (b-discsqrt)*(-0.5) + else (b+discsqrt)*(-0.5) + t0' = q/a + t1' = c/q + t0 = fmin t0' t1' + t1 = fmax t0' t1' + in if t1 < 0 || t0 > d + then RayMiss + else let dist = if t0 < 0 + then t1 + else t0 + in if dist < 0 || dist > d + then RayMiss + else + let pos = vscaleadd orig dir dist + Vec posx posy posz = pos + in if posz > clip1 && posz < clip2 + then let invhyp = 1 / (sqrt (height*height + r*r)) + up = r * invhyp + out = height * invhyp + r_ = sqrt (posx*posx + posy*posy) + correction = (out)/(r_) + in RayHit dist pos (Vec (posx*correction) (posy*correction) up) t + else + if dz > 0 -- ray pointing up from bottom + then if oz < clip1 + then rayint (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d t + else RayMiss + else if oz > clip2 + then let r2 = r*(1-((clip2-clip1)/(height))) + in rayint (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d t + --rayint_aadisc clip2 r2 (Ray orig dir) d t + else RayMiss + -- then rayint_aadisc clip1 r (Ray orig dir) d t + -- else RayMiss -- rayint_aadisc clip2 + -- (r*((clip2-clip1)/height)) + -- (Ray orig dir) d t -- todo: fix normal + +--Plane-- +rayint (Plane norm offset) (Ray orig dir) d t = + let hit = -(((vdot norm orig)-offset) / (vdot norm dir)) + in if hit < 0 || hit > d + then RayMiss + else RayHit hit (vscaleadd orig dir hit) norm t + +--Box-- +-- this could be optimized a bit more +rayint (Box (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z))) r d t = + let (Ray orig dir) = r + (Vec ox oy oz) = orig + (Vec dx dy dz) = dir + dxrcp = 1/dx + dyrcp = 1/dy + dzrcp = 1/dz + Interval inx outx = if dx > 0 + then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp) + else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp) + Interval iny outy = if dy > 0 + then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp) + else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp) + Interval inz outz = if dz > 0 + then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp) + else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp) + lastin = (fmax3 inx iny inz) + firstout = (fmin3 outx outy outz) + in if lastin > firstout || firstout < 0 || lastin > d + then RayMiss + else + let n = if inx == lastin + then if dx > 0 then nvx else vx + else if iny == lastin + then if dy > 0 then nvy else vy + else if dz > 0 then nvz else vz + norm = if lastin > 0 then n else vinvert n + hitdepth = fmax 0 lastin + in + RayHit hitdepth (vscaleadd orig dir hitdepth) norm t + +--Composite objects-- +--Instance-- +-- transforming the distance is a little awkward +-- the normal shouldn't have to be re-normalized, should it? +rayint (Instance sld xfm) (Ray orig dir) d t = + let newdir = invxfm_vec xfm dir + neworig = invxfm_point xfm orig + lenscale = vlen newdir + invlenscale = 1/lenscale + in + case (rayint_check sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale) t) of + RayMiss -> RayMiss + RayHit depth pos n tex -> RayHit (depth*invlenscale) (xfm_point xfm pos) (vnorm (invxfm_norm xfm n)) tex + +--Group-- +-- rayint (Group lst) r d t = foldl nearest RayMiss (map (\x -> rayint x r d t) lst) +rayint (Group xs) r d t = + let rig [] = RayMiss + rig (x:xs) = nearest (rayint_check x r d t) (rig xs) + in rig xs + +--Difference-- +-- csg of object a - object b -- +rayint (Difference sa sb) r d t = + let dif = Difference sa sb + Ray orig dir = r + ria = rayint sa r d t + in + case ria of + RayMiss -> RayMiss + RayHit ad ap an at -> + if inside sb orig + then + case rayint sb r d t of + RayMiss -> RayMiss + RayHit bd bp bn bt -> + if bd < ad + then if inside sa bp + then RayHit bd bp (vinvert bn) bt + else rayint_advance dif r d t bd + else rayint_advance dif r d t bd + else + if inside sb ap + then rayint_advance dif r d t ad + else RayHit ad ap an at + + +--Intersection-- +-- fixme: there's some numerical instability near edges +rayint (Intersection slds) r d t = + let (Ray orig dir) = r + in + if null slds || d < 0 + then RayMiss + else + let s = head slds + in case tail slds of + [] -> rayint_check s r d t + ss -> if inside s orig + then case rayint s r d t of + RayMiss -> rayint (Intersection ss) r d t + RayHit sd sp sn st -> + case rayint_check (Intersection ss) r sd t of + RayMiss -> rayint_advance (Intersection slds) + r d t sd + hit -> hit + else case rayint s r d t of + RayMiss -> RayMiss + RayHit sd sp sn st -> + if inside (Intersection ss) sp + then RayHit sd sp sn st + else rayint_advance (Intersection slds) + r d t sd + + +rayint (Bound sa sb) r d t = + let (Ray orig _) = r + in if inside sa orig || shadow sa r d + then rayint sb r d t + else RayMiss + +--Bih-- +rayint (Bih bb root) r d t = + let Ray orig dir = r + dir_rcp = vrcp dir + Interval near far = bbclip r bb + traverse (BihLeaf s) near far = rayint_check s r (fmin d far) t + traverse (BihBranch lsplit rsplit axis l r) near far = + let dirr = va dir_rcp axis + o = va orig axis + dl = (lsplit - o) * dirr + dr = (rsplit - o) * dirr + in + if near > far + then RayMiss + else + -- this is ugly and verbose, + -- but it does what it needs to + if dirr > 0 + then + (nearest + (if near < dl + then traverse l near (fmin dl far) + else RayMiss) + (if dr < far + then traverse r (fmax dr near) far + else RayMiss)) + else + (nearest + (if near < dr + then traverse r near (fmin dr far) + else RayMiss) + (if dl < far + then traverse l (fmax dl near) far + else RayMiss)) --} + + in + traverse root near far + +-- anything that hits it gets teleported +-- one-way door to another world +-- this is dangerous and doesn't really work right +rayint (Portal port world) r d t = + case rayint port r d t of + RayMiss -> RayMiss + RayHit depth _ _ _ -> + rayint_advance world r d t depth + +--Tex-- +-- this is a little odd; rather than associate +-- a texture with each primitive, we use a +-- container object; everything inside has that +-- texture, unless it's overridden by a nested +-- texture +rayint (Tex s tex) r d t = rayint_check s r d tex + +--Nothing-- +rayint (Solid.Nothing) _ _ _ = RayMiss + +-- default case: miss +-- rayint _ _ _ _ = RayMiss + +-- various specialized ray intersections, used as helper functions + +-- used by cylinder / cone code +-- broken, do not use +rayint_aadisc :: Flt -> Flt-> Ray -> Flt -> Texture -> Rayint +rayint_aadisc height radius r d t = + let Ray orig dir = r + -- Vec _ _ oz = orig + -- Vec _ _ dz = dir + -- dist = (height-oz)/dz + dist = plane_int_dist r (Vec 0 0 height) vx + in if dist <= 0 || dist >= d || isNaN dist + then RayMiss + else let pos = vscaleadd orig dir dist + (Vec px py _) = pos + in + if (px*px + py+py) > radius*radius + then RayMiss + else RayHit dist pos vz t + +-- move ray forward, intersect, fix result +-- useful in csg +rayint_advance :: Solid -> Ray -> Flt -> Texture -> Flt -> Rayint +rayint_advance s r d t adv = + let a = adv+delta + in + case (rayint s (ray_move r a) (d-a) t) of + RayMiss -> RayMiss + RayHit depth pos norm tex -> RayHit (depth+a) pos norm tex + +-- check results of a ray-intersection test +rayint_check s r d t = + let Ray orig dir = r in + case rayint s r d t of + RayMiss -> RayMiss + RayHit depth pos norm tex -> + if depth < 0 + then error $ "rayint depth < 0 " ++ (show depth) ++ " " ++ (show s) + else + if depth > d + then error $ "rayint depth (" ++ (show depth) ++ ") > d (" ++ (show d) ++ ") " ++ (show s) + else + if not $ veq pos (vscaleadd orig dir depth) + then error $ "rayint position doesn't match depth" ++ + (show pos) ++ (show $ vscaleadd orig dir depth) ++ (show depth) ++ (show s) + else + if (vlen norm) < 1-delta + then error "normal too short" + else + if (vlen norm) > 1+delta + then error $ "normal too long " ++ (show norm) ++ " " ++ (show s) + else RayHit depth pos norm tex + + +--SHADOW-- +shadow :: Solid -> Ray -> Flt -> Bool + +--Sphere-- +shadow (Sphere center r invr) (Ray e dir) dist = + let eo = vsub center e + v = vdot eo dir + in + if (dist >= (v - r)) && (v > 0.0) + then + let vsqr = v*v + csqr = vdot eo eo + rsqr = r*r + disc = rsqr - (csqr - vsqr) in + if disc < 0.0 then + False + else + let d = sqrt disc + hitdist = if (v-d) > 0 then (v-d) else (v+d) + in if (hitdist < 0) || (hitdist > dist) + then False + else True + else + False + + +shadow (Triangle p1 p2 p3) (Ray o dir) dist = + let e1 = vsub p2 p1 + e2 = vsub p3 p1 + s1 = vcross dir e2 + divisor = vdot s1 e1 + in + if (divisor == 0) + then False + else + let invdivisor = 1.0 / divisor + d = vsub o p1 + b1 = (vdot d s1) * invdivisor + in + if (b1 < 0) || (b1 > 1) + then False + else + let s2 = vcross d e1 + b2 = (vdot dir s2) * invdivisor + in + if (b2 < 0) || (b1 + b2 > 1) + then False + else + let t = (vdot e2 s2) * invdivisor + in + if (t < 0) || (t > dist) + then False + else True + +shadow (TriangleNorm p1 p2 p3 n1 n2 n3) r d = + shadow (Triangle p1 p2 p3) r d + +shadow (Disc point norm radius_sqr) r d = + let (Ray orig dir) = r + dist = plane_int_dist r point norm + in if dist < 0 || dist > d + then False + else let pos = vscaleadd orig dir dist + offset = vsub pos point + in + if (vdot offset offset) > radius_sqr + then False + else True + +shadow (Box box) r d = + let Interval near far = bbclip r box + in + if (near > far) || far <= 0 || far > d + then False + else True + +-- cone centered on z axiz, height of hp, clipped at h1 and h2 +shadow (Cone r clip1 clip2 height) (Ray orig dir) d = + let Vec ox oy oz = orig + Vec dx dy dz = dir + k' = (r/height) + k = k'*k' + a = dx*dx + dy*dy - k*dz*dz + b = 2*(dx*ox + dy*oy - k*dz*(oz-height)) + c = ox*ox + oy*oy - k*(oz-height)*(oz-height) + disc = b*b - 4*a*c + in if disc < 0 + then False + else + let discsqrt = sqrt disc + q = if b < 0 + then (b-discsqrt)*(-0.5) + else (b+discsqrt)*(-0.5) + t0' = q/a + t1' = c/q + t0 = fmin t0' t1' + t1 = fmax t0' t1' + in if t1 < 0 || t0 > d + then False + else let dist = if t0 < 0 + then t1 + else t0 + in if dist < 0 || dist > d + then False + else + let pos = vscaleadd orig dir dist + Vec posx posy posz = pos + in if posz > clip1 && posz < clip2 + then True + else + if dz > 0 -- ray pointing up from bottom + then if oz < clip1 + then shadow (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d + else False + else if oz > clip2 + then let r2 = r*(1-((clip2-clip1)/(height))) + in shadow (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d + else False + +shadow (Instance sld xfm) (Ray orig dir) d = + let newdir = invxfm_vec xfm dir + neworig = invxfm_point xfm orig + lenscale = vlen newdir + invlenscale = 1/lenscale + in + shadow sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale) + +shadow (Tex s t) r d = shadow s r d + +shadow (Group xs) r d = + let sg [] = False + sg (x:xs) = (shadow x r d) || (sg xs) + in sg xs + +shadow (Bound sa sb) r d = + let (Ray orig _ ) = r + in if inside sa orig || shadow sa r d + then shadow sb r d + else False + +shadow (Bih bb root) r d = + let (Ray orig dir) = r + dir_rcp = vrcp dir + Interval near far = bbclip r bb + traverse (BihLeaf s) near far = shadow s r (fmin d far) + traverse (BihBranch lsplit rsplit axis l r) near far = + let dirr = va dir_rcp axis + o = va orig axis + dl = (lsplit - o) * dirr + dr = (rsplit - o) * dirr + in + if near > far + then False + else + if dirr > 0 + then + ((if near < dl + then traverse l near (fmin dl far) + else False) + || + (if dr < far + then traverse r (fmax dr near) far + else False)) + else + ((if near < dr + then traverse r near (fmin dr far) + else False) + || + (if dl < far + then traverse l (fmax dl near) far + else False)) + + in traverse root near far + +-- default shadow test, in case an optimized test +-- isn't implemented yet, we use the regular +-- trace function; we have to be careful with +-- container objects, though; using the slow +-- test on the container means doing the slow +-- test against everything it contains +shadow s r d = + case (rayint s r d t_white) of + RayHit _ _ _ _ -> True + RayMiss -> False + +--INSIDE-- +inside :: Solid -> Vec -> Bool +inside (Sphere center r invr) pt = + let offset = vsub center pt + in (vdot offset offset) < r*r + +inside (Group slds) pt = + foldl' (||) False (map (\x -> inside x pt) slds) + +inside (Difference sa sb) pt = + (inside sa pt) && (not $ inside sb pt) + +-- note: inside is True for an empty intersection. +-- this is actually the preferred semantics in +-- some cases, strange as it may seem. +inside (Intersection slds) pt = + foldl' (&&) True (map (\x -> inside x pt) slds) + +inside (Instance s xfm) pt = + inside s (xfm_point xfm pt) + +inside (Plane norm offset) pt = + let onplane = (vscale norm offset) + newvec = vsub onplane pt + in vdot newvec norm > 0 + +inside (Box (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2))) (Vec x y z) = + x > x1 && x < x2 && y > y1 && y < y2 && z > z1 && z < z2 + +inside (Tex s t) pt = inside s pt + +inside (Bih (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2)) root) pt = + let (Vec x y z) = pt + traverse (BihLeaf s) = inside s pt + traverse (BihBranch lsplit rsplit axis l r) = + let o = va pt axis + in (if o < lsplit + then (traverse l) + else False) + || + (if o > rsplit + then (traverse r) + else False) + in + (x > x1) && (x < x2) && (y > y1) && (y < y2) && (z > z1) && (z < z2) && (traverse root) + +inside (Cylinder r h1 h2) (Vec x y z) = + z > h1 && z < h2 && x*x + y*y < r*r + +inside (Cone rbase h1 h2 height) (Vec x y z) = + let r = rbase*(1-(((z-h1)/height))) + in z > h1 && z < h2 && x*x + y*y < r*r + +inside (Bound sa sb) pt = inside sa pt && inside sb pt + +inside _ _ = False + +-- return distance to surface, positive if inside, negative if outside +-- this isn't used for anything in particular yet +power :: Solid -> Vec -> Flt + +power (Sphere center r invr) pt = + let offset = vsub center pt + in r - (vlen offset) + +power (Group slds) pt = + foldl' (max) (-infinity) (map (\x -> power x pt) slds) + +-- not accurate +power (Instance s xfm) pt = + power s (xfm_point xfm pt) + +power (Plane norm offset) pt = + let onplane = (vscale norm offset) + newvec = vsub onplane pt + in vdot newvec norm + + +--BOUND-- +bound :: Solid -> Bbox +bound (Sphere center r invr) = + let offset = (vec r r r) in + (Bbox (vsub center offset) (vadd center offset)) + +bound (Triangle (Vec v1x v1y v1z) + (Vec v2x v2y v2z) + (Vec v3x v3y v3z)) = + Bbox + (Vec ((fmin (fmin v1x v2x) v3x) - delta) + ((fmin (fmin v1y v2y) v3y) - delta) + ((fmin (fmin v1z v2z) v3z) - delta) ) + + (Vec ((fmax (fmax v1x v2x) v3x) + delta) + ((fmax (fmax v1y v2y) v3y) + delta) + ((fmax (fmax v1z v2z) v3z) + delta) ) + +bound (TriangleNorm v1 v2 v3 n1 n2 n3) = + bound (Triangle v1 v2 v3) + +--dangerous to use inside a bih +bound (Plane norm offset) = everything_bbox + +bound (Box box) = box + +-- this could be a tighter fit +bound (Disc pos norm rsqr) = + bound (sphere pos (sqrt rsqr)) + +bound (Cylinder r h1 h2) = + Bbox (Vec (-r) (-r) h1) (Vec r r h2) + +bound (Cone r h1 h2 height) = + Bbox (Vec (-r) (-r) h1) (Vec r r h2) + +bound (Group slds) = + foldl' bbjoin empty_bbox (map bound slds) + +bound (Difference sa sb) = bound sa + + +bound (Intersection slds) = + if null slds + then empty_bbox + else foldl' bboverlap everything_bbox (map bound slds) + + + +-- the reason the following doesn't work: +-- bound (Bound sa sb) = bboverlap (bound sa) (bound sb) +-- is that the ray has to hit solid a, and it might not +-- if the bound is smaller than sa. + +-- update: this is probably no longer true, as I added an "inside" +-- test to the bounds check as well + +bound (Bound sa sb) = bound sa + +-- not optimal, but it does the job +bound (Instance sld xfm) = + let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bound sld + pxfm = xfm_point xfm + in + bbpts [(pxfm (Vec x y z)) | x <- [p1x,p2x], + y <- [p1y,p2y], + z <- [p1z,p2z]] + + +bound (Bih bb root) = bb + +bound (Portal port world) = bound port + +bound (Tex s t) = bound s + +bound (Photon p indir clr r) = bound (sphere p r) + +bound Solid.Nothing = empty_bbox + +-- no default case: we want an exception +-- if there is no match + + +--CAMERA-- +data Camera = Camera {campos, fwd, up, right :: !Vec} + deriving Show + +default_cam = (Camera (vec 0.0 0.0 (-3.0)) + (vec 0.0 0.0 1.0) + (vec 0.0 1.0 0.0) + (vec 1.0 0.0 0.0) ) + +camera :: Vec -> Vec -> Vec -> Flt -> Camera +camera pos at up angle = + let fwd = vnorm $ vsub at pos + right = vnorm $ vcross up fwd + up_ = vnorm $ vcross fwd right + cam_scale = tan ((pi/180)*(angle/2)) + in + Camera pos fwd + (vscale up_ cam_scale) + (vscale right cam_scale) + +--SCENE-- +data Scene = Scene {sld :: !Solid, + lights :: ![Light], + cam :: !Camera, + dtex :: !Texture, + bground :: !Color} deriving Show + +default_scene = (Scene (sphere (vec 0.0 0.0 0.0) 1.0) + [] default_cam t_white c_white) diff --git a/SolidTexture.hs b/SolidTexture.hs new file mode 100644 --- /dev/null +++ b/SolidTexture.hs @@ -0,0 +1,127 @@ +{- +module SolidTexture (square_wave, triangle_wave, sine_wave, + stripe, noise, turbulence + ) where -} + +module SolidTexture where +import Clr +import Vec +import Data.Array.IArray + +-- INTERPOLATION FUNCTIONS -- +square_wave :: Flt -> Flt +square_wave x = + let offset = x - (fromIntegral (floor x)) + in if offset < 0.5 then 0 else 1 + +triangle_wave :: Flt -> Flt +triangle_wave x = + let offset = x - (fromIntegral (floor x)) + in if offset < 0.5 + then (offset*2) + else (2-(offset*2)) + +sine_wave :: Flt -> Flt +sine_wave x = (sin (x*2*pi))*0.5 + 0.5 + + +lump_wave :: Flt -> Flt +lump_wave x = 1 - x*x*x + +-- SCALAR TEXTURE FUNCTIONS -- + +-- These are simple solid texture functions that take a +-- point as argument and return a number 0 < n < 1 + +stripe :: Vec -> (Flt -> Flt) -> (Vec -> Flt) +stripe axis interp = + let len = vlen axis + in + (\pos -> let offset = vdot pos axis + in interp offset) + + +-- PERLIN NOISE -- + +-- (-6 t^5 + 15 t^4 - 10t^3 +1) +-- "realistic ray tracing 2nd edition" inconsistent +-- on whether it should be t^5 or t^6, +-- but t^5 works and t^6 doesn't. +omega :: Flt -> Flt +omega t_ = + let t = abs t_ + tsqr = t*t + tcube = tsqr*t + in (-6)*tcube*tsqr + 15*tcube*t - 10*tcube + 1 + +-- questionably random +phi :: Array Int Int +phi = listArray (0,11) [3,0,2,7,4,1,5,11,8,10,9,6] + +grad :: Array Int Vec +grad = listArray (0,11) + $ filter (\x -> let l = vlen x in l < 1.5 && l > 1.1) + [Vec x y z | x <- [(-1),0,1], + y <- [(-1),0,1], + z <- [(-1),0,1]] + +gamma :: Int -> Int -> Int -> Vec +gamma i j k = + let a = phi!(mod (abs k) 12) + b = phi!(mod (abs (j+a)) 12) + c = phi!(mod (abs (i+b)) 12) + in grad!c + +knot :: Int -> Int -> Int -> Vec -> Flt +knot i j k v = + let Vec x y z = v + in (omega x) * (omega y) * (omega z) * (vdot (gamma i j k) v) + +intGamma :: Int -> Int -> Int +intGamma i j = + let a = phi!(mod (abs j) 16) + b = phi!(mod (abs (i+a)) 16) + in b + +turbulence :: Vec -> Int -> Flt +turbulence p 1 = abs(noise(p)) +turbulence p n = + let newp = vscale p 0.5 + t = abs (noise p) + in t + (0.5 * (turbulence newp (n-1))) + +noise :: Vec -> Flt +noise (Vec x y z) = + let i = floor x + j = floor y + k = floor z + u = x-(fromIntegral i) + v = y-(fromIntegral j) + w = z-(fromIntegral k) + in knot i j k (Vec u v w) + + knot (i+1) j k (Vec (u-1) v w) + + knot i (j+1) k (Vec u (v-1) w) + + knot i j (k+1) (Vec u v (w-1)) + + knot (i+1) (j+1) k (Vec (u-1) (v-1) w) + + knot (i+1) j (k+1) (Vec (u-1) v (w-1)) + + knot i (j+1) (k+1) (Vec u (v-1) (w-1)) + + knot (i+1) (j+1) (k+1) (Vec (u-1) (v-1) (w-1)) + +perlin :: Vec -> Flt +perlin v = + let p = ((noise v)+1)*0.5 + in if p > 1 + then error $ "perlin noise error, 1 < " ++ (show p) + else if p < 0 + then error $ "perlin noise error, 0 > " ++ (show p) + else p + +--untested +perlin_turb :: Vec -> Int -> Flt +perlin_turb v l = + let p = turbulence v l + in if p > 1 + then error $ "perlin turbulence error, 1 < " ++ (show p) + else if p < 0 + then error $ "perlin turbulence error, 0 > " ++ (show p) + else p diff --git a/Spd.hs b/Spd.hs new file mode 100644 --- /dev/null +++ b/Spd.hs @@ -0,0 +1,246 @@ +module Spd where +import Vec +import Clr +import Solid + +-- NFF file format description: +-- http://tog.acm.org/resources/SPD/NFF.TXT + +-- this would probably be much shorter if I used scanf instead of lex + +lexignore s = + let t = lex s + in + case t of + [] -> [] -- newline + [("",s1)] -> lexcr s1 + [(_,s1)] -> lexignore s1 + +lexcr s = + let t = lex s + in + case t of + [(s1,[])] -> t + [("",s1)] -> lexcr s1 + [("#",s1)] -> lexignore s1 + _ -> t + +data BgColor = BgColor(Color) + +readsSpdVec :: ReadS Vec +readsSpdVec s = [((Vec x y z),s3) | (x,s1) <- reads s :: [(Flt,String)], + (y,s2) <- reads s1 :: [(Flt,String)], + (z,s3) <- reads s2 :: [(Flt,String)] ] +instance Read Vec where + readsPrec _ = readsSpdVec + +readsSpdVecNorm :: ReadS (Vec,Vec) +readsSpdVecNorm s = [(((Vec x y z),(Vec nx ny nz)),s6) + | (x,s1) <- reads s :: [(Flt,String)], + (y,s2) <- reads s1 :: [(Flt,String)], + (z,s3) <- reads s2 :: [(Flt,String)], + (nx,s4) <- reads s3 :: [(Flt,String)], + (ny,s5) <- reads s4 :: [(Flt,String)], + (nz,s6) <- reads s5 :: [(Flt,String)] ] + +instance Read (Vec,Vec) where + readsPrec _ = readsSpdVecNorm + + +-- if this seems intuitive, there's something wrong with you +readsSpdVecs :: ReadS [Vec] +readsSpdVecs s = + let parses = reads s :: [(Vec,String)] + in + if null parses + then [([],s)] + else + let (v,rest) = head parses + (vs,returns) = head (readsSpdVecs rest) + in [((v:vs),returns)] + +instance Read [Vec] where + readsPrec _ = readsSpdVecs + +readsSpdVecsNorms :: ReadS [(Vec,Vec)] +readsSpdVecsNorms s = + let parses = readsSpdVecNorm s :: [((Vec,Vec),String)] + in + if null parses + then [([],s)] + else + let (v,rest) = head parses + (vs,returns) = head (readsSpdVecsNorms rest) + in [((v:vs),returns)] + +instance Read [(Vec,Vec)] where + readsPrec _ = readsSpdVecsNorms + +{- "v" + "from" Fx Fy Fz + "at" Ax Ay Az + "up" Ux Uy Uz + "angle" angle + "hither" hither + "resolution" xres yres -} +readsSpdCam :: ReadS Camera +readsSpdCam s = [ (camera from at up angle,s14) | ("v", s1) <- lexcr s, + ("from", s2) <- lexcr s1, + (from, s3) <- reads s2 :: [(Vec,String)], + ("at", s4) <- lexcr s3, + (at, s5) <- reads s4 :: [(Vec,String)], + ("up", s6) <- lexcr s5, + (up, s7) <- reads s6 :: [(Vec,String)], + ("angle", s8) <- lexcr s7, + (angle, s9) <- reads s8 :: [(Flt,String)], + ("hither", s10)<- lexcr s9, + (_,s11) <- lexcr s10, + ("resolution", s12) <- lexcr s11, + (_, s13) <- lexcr s12, + (_, s14) <- lexcr s13 ] +instance Read Camera where + readsPrec _ = readsSpdCam + + + +readsSpdClr :: ReadS Color +readsSpdClr s = [((Color r g b), s3) | (r, s1) <- reads s :: [(Flt,String)], + (g, s2) <- reads s1 :: [(Flt,String)], + (b, s3) <- reads s2 :: [(Flt,String)] ] +instance Read Color where + readsPrec _ = readsSpdClr + + +-- "b" R G B +readsSpdBackground :: ReadS BgColor +readsSpdBackground s = [((BgColor clr), s2) | ("b", s1) <- lexcr s, + (clr, s2) <- reads s1 :: [(Color,String)] ] +instance Read BgColor where + readsPrec _ = readsSpdBackground + + +-- "l" X Y Z [R G B] +readsSpdLight :: ReadS Light +readsSpdLight s = [((Light pos clr),s3) | ("l", s1) <- lexcr s, + (pos, s2) <- reads s1 :: [(Vec,String)], + (clr, s3) <- reads s2 :: [(Color,String)] ] + ++ + [((Light pos (Color 1 1 1)),s2) | ("l", s1) <- lexcr s, + (pos, s2) <- reads s1 :: [(Vec,String)] ] +instance Read Light where + readsPrec _ = readsSpdLight + +-- "f" red green blue Kd Ks Shine T index_of_refraction +-- data Material = Material {clr :: Color, reflect, refract, ior, kd, ks, shine :: Flt} +readsSpdFill :: ReadS Texture +readsSpdFill s = [(\ri->Material clr ks (1-trans) ior kd shine, s7) | ("f", s1) <- lexcr s, + (clr, s2) <- reads s1 :: [(Color,String)], + (kd, s3) <- reads s2 :: [(Flt,String)], + (ks, s4) <- reads s3 :: [(Flt,String)], + (shine, s5) <- reads s4 :: [(Flt,String)], + (trans, s6) <- reads s5 :: [(Flt,String)], + (ior, s7) <- reads s6 :: [(Flt,String)] ] + +instance Read (Rayint -> Material) where + readsPrec _ = readsSpdFill + + +-- "s" center.x center.y center.z radius + +-- "c" +-- base.x base.y base.z base_radius +-- apex.x apex.y apex.z apex_radius + +-- "p" total_vertices +-- vert1.x vert1.y vert1.z +-- [etc. for total_vertices vertices] + +readsSpdSolid :: ReadS Solid +readsSpdSolid s = [((sphere center radius),s3) | ("s", s1) <- lexcr s, + (center,s2) <- reads s1 :: [(Vec,String)], + (radius,s3) <- reads s2 :: [(Flt,String)] ] + ++ + [((cone e1 r1 e2 r2),s5) | ("c",s1) <- lexcr s, + (e1,s2) <- reads s1 :: [(Vec,String)], + (r1,s3) <- reads s2 :: [(Flt,String)], + (e2,s4) <- reads s3 :: [(Vec,String)], + (r2,s5) <- reads s4 :: [(Flt,String)] ] + ++ + [(group (triangles verts),s3) | ("p",s1) <- lexcr s, + (n,s2) <- reads s1 :: [(Int,String)], + (verts,s3) <- readsSpdVecs s2 :: [([Vec],String)] ] + ++ + [(group (trianglesnorms (vns)),s3) | ("pp",s1) <- lexcr s, + (n,s2) <- reads s1 :: [(Int,String)], + (vns,s3) <- readsSpdVecsNorms s2 :: [([(Vec,Vec)],String)] ] + {- ++ + [(Tex(Solid.Nothing,t),s1) | (t,s1) <- reads s :: [(Texture,String)]] -} + + +-- instance Read Solid where +-- readsPrec _ = readsSpdSolid + + +-- same as readSpdVecs, just different types +readsSpdPrims :: ReadS [Solid] +readsSpdPrims s = + let parses = readsSpdSolid s :: [(Solid,String)] + in + if null parses + then [([],s)] + else + let (v,rest) = head parses + (vs,returns) = head (readsSpdPrims rest) + in [((v:vs),returns)] + +instance Read [Solid] where + readsPrec _ = readsSpdPrims + + +readsSpdTextureGroup :: ReadS Solid +readsSpdTextureGroup s = + [((Tex (bih prims) t),s2) | (t,s1) <- reads s :: [(Texture,String)], + (prims,s2) <- readsSpdPrims s1 :: [([Solid],String)] ] + +instance Read Solid where + readsPrec _ = readsSpdTextureGroup + +accum_rss :: [Camera] -> [Light] -> [Solid] -> [BgColor] -> String -> ([Camera],[Light],[Solid],[BgColor],String) +accum_rss cams lights prims background s = + if null s + then (cams,lights,prims,background,s) + else + let cam = reads s :: [(Camera,String)] + sld = reads s :: [(Solid,String)] + lit = reads s :: [(Light,String)] + bgc = reads s :: [(BgColor,String)] + in + if not $ null cam + then + let (c1,s1) = head cam + in accum_rss (c1:cams) lights prims background s1 + else + if not $ null sld + then + let (s2,s1) = head sld + in accum_rss cams lights (s2:prims) background s1 + else + if not $ null lit + then + let (l1,s1) = head lit + in accum_rss cams (l1:lights) prims background s1 + else + if not $ null bgc + then + let (b1,s1) = head bgc + in accum_rss cams lights prims (b1:background) s1 + else + (cams,lights,prims,background,s) + +readsSpdScene :: ReadS Scene +readsSpdScene s = + let ((cam:cams),lights,prims,(BgColor(bgc):bgcs),s1) = accum_rss [] [] [] [] s + in [((Scene (bih prims) lights cam t_white bgc),s1)] + +instance Read Scene where + readsPrec _ = readsSpdScene diff --git a/TestScene.hs b/TestScene.hs new file mode 100644 --- /dev/null +++ b/TestScene.hs @@ -0,0 +1,168 @@ +module TestScene where +import Vec +import Solid +import Clr +import Data.List hiding (group) +import SolidTexture +import System.Random + +{- +lits = [Light {litpos = Vec {x = -3.0, y = 1.7, z = 5.0}, + litcol = Color 0.9 1 1 }, + Light {litpos = Vec {x = 1.1, y = -4.2, z = 4.0}, + litcol = Color 1 0.9 1 }, + Light {litpos = Vec {x = 4.3, y = 3.0, z = 2.3}, + litcol = Color 1 1 0.9}] +-} + +lits = [ Light (Vec (-100) 70 (140)) (cscale (Color 1 0.8 0.8) 1500) + , Light (Vec (-3) 5 8) (Color 1.5 2 2) + ] + +icosahedron pos r = + let gr = ((1+(sqrt 5))/2) -- golden ratio, 1.618033988749895 + n11 = [(-r),r] + ngrgr = [(-gr)*r,gr*r] + grrcp = [((-r)/gr),(r/gr)] + points = [Vec x y z | x <- n11, + y <- n11, + z <- n11 ] ++ + [Vec 0 y z | y <- grrcp, + z <- ngrgr ] ++ + [Vec x y 0 | x <- grrcp, + y <- ngrgr] ++ + [Vec x 0 z | x <- ngrgr, + z <- grrcp] + pln x = (Plane (vnorm x) (r+(vdot (vnorm x) pos))) + in + Intersection ((sphere pos (1.26*r)):(map pln points)) + +dodecahedron pos r = + let gr = (1+(sqrt 5))/2 -- golden ratio, 1.618033988749895 + n11 = [(-r),r] + ngrgr = [(-gr)*r,gr*r] + points = [Vec 0 y z | y <- n11, z <- ngrgr] ++ + [Vec x 0 z | z <- n11, x <- ngrgr] ++ + [Vec x y 0 | x <- n11, y <- ngrgr] + pln x = (Plane (vnorm x) (r+(vdot (vnorm x) pos))) + in + Intersection ((sphere pos (1.26*r)):(map pln points)) + +cust_cam = camera (vec (-2) (4.3) (15)) (vec 0 2 0) (vec 0 1 0) 45 + + +spiral = [ ((Vec ((sin (rot n))*n) + ((cos (rot n))*n) + (n-3)), (n/15)) | n <- [0, 0.01..6]] + + +coil = bih (zipWith (\ (p1,r1) (p2,r2) -> (Solid.group [(cone p1 r1 p2 r2), + (sphere p1 r1)] )) + spiral + (tail spiral)) + +-- we branch once per year +-- not really a plausible oak, but it's getting there +oak :: Flt -> StdGen -> Solid +oak age rng = + if age < 0 + then Solid.Nothing + else + let year :: Int = floor age + season = age-(fromIntegral year) + thickness = 0.03 + minbranch = deg 10 + maxbranch = deg 25 + tree 0 r = Solid.Nothing + tree 1 r = -- cone (Vec 0 0 0) thickness (Vec 0 season 0) 0 + Tex (sphere (Vec 0 0 0) season) (t_matte (Color 0.2 1 0.4)) + tree n r = let nf = fromIntegral n + height = nf + (rng1,rng2) = split r + (rng3,rng4) = split rng1 + (r1,rng5) = randomR (0,0.5) rng4 + (r2,rng6) = randomR (minbranch,maxbranch) rng5 + (r3,rng7) = randomR (0.8,0.95) rng6 + seglen = 0.5 + r1 + branchang = r2 + scaling = r3 + -- we make our own manual bounding heirarchy + -- (bih doesn't know what to do with heirachies + -- of transformed objects) + in Bound (sphere (Vec 0 (height/2) 0) (height/2)) + (group [ cone (Vec 0 0 0) (thickness*nf) (Vec 0 seglen 0) (thickness*(nf-1)*scaling) + , transform (tree (n-1) rng2) [(scale (Vec scaling scaling scaling)), + (rotate (Vec 0 0 1) branchang), + (rotate (Vec 0 1 0) (deg 30)), + (translate (Vec 0 seglen 0))] + , transform (tree (n-1) rng3) [(scale (Vec scaling scaling scaling)), + (rotate (Vec 0 0 1) (-branchang)), + (rotate (Vec 0 1 0) (deg 30)), + (translate (Vec 0 seglen 0))] + ]) + in Tex (bih (flatten_transform (tree year rng))) (t_matte (Color 0.8 0.5 0.4)) + +sphereint = Intersection [ (sphere (Vec (-1) 0 0) 2), + (sphere (Vec 1 0 0) 2), + (sphere (Vec 0 (-1) 0) 2), + (sphere (Vec 0 1 0) 2) ] + +-- scene with everything +-- can't have the plane inside a bih structure because it's infinite +geom = group [ Tex (plane (Vec 0 0 0) (Vec 0 1 0)) (t_matte (Color 0 0.8 0.3)) + , bih [ Tex (dodecahedron (Vec (-6) 3 0) 1) t_stripe + , Tex (icosahedron (Vec 4 1.5 3) 1.5) t_mottled + --,(oak 11.6 (mkStdGen 42)) + , transform (oak 11.6 (mkStdGen 42)) [ scale (Vec 1.5 1.5 1.5)] + , Tex (transform (coil) [ scale (Vec (1/3) (1/3) (1/3)) + , rotate (Vec 0 1 0) (deg 65) + , translate (Vec (-3.5) 1 (5)) + ]) t_mirror --} + , cone (Vec (-6) 0 0) 1 (Vec (-6) 3 0) 0 + , Tex (Difference (sphere (Vec 0 (-4) 5) 4.7) (sphere (Vec 1.5 (1.5) 5.2) 1.6)) t_mirror + , transform (Tex sphereint (t_matte (Color 0.5 0 1))) [ scale (Vec 0.6 0.6 0.6), + translate (Vec (-5.2) 1 5)] + ] + ] + +geom1 = transform (oak 7.2 (mkStdGen 42)) [ scale (Vec 1.5 1.5 1.5)] + +-- color reflect refract ior kd shine +m_shiny_white :: Material +m_shiny_white = (Material c_white 0.3 0 0 0.7 10) + +t_shiny_white :: Texture +t_shiny_white ri = m_shiny_white + +m_dull_gray :: Material +m_dull_gray = (Material (Color 0.4 0.3 0.35) 0 0 0 0.2 1) + +t_mottled (RayHit _ pos norm _) = + --let scale = (stripe (Vec 1 1 1) sine_wave) pos + let scale = perlin (vscale pos 3) + in if scale < 0 then error "foo" + else if scale > 1 + then error "bar" + else m_interp m_shiny_white m_dull_gray scale + +--shouldn't happen +t_mottled RayMiss = m_shiny_white + +t_stripe (RayHit _ pos norm _) = + let scale = (stripe (Vec 2 4 3) sine_wave) pos + in if scale < 0 then error "foo" + else if scale > 1 + then error "bar" + else m_interp m_shiny_white m_dull_gray scale + +--shouldn't happen +t_stripe RayMiss = m_shiny_white + +t_matte c = + (\ri -> (Material c 0 0 0 1 2)) + +t_mirror = + (\ri -> (Material (Color 0.8 0.8 1) 1 0 0 0.2 100)) + +scn :: IO Scene +scn = return (Scene geom lits cust_cam (t_matte (Color 0.8 0.5 0.4)) (Color 0.4 0.5 0.8)) diff --git a/Trace.hs b/Trace.hs new file mode 100644 --- /dev/null +++ b/Trace.hs @@ -0,0 +1,106 @@ +module Trace where +import Vec +import Clr +import Solid +import Data.List +import Control.Concurrent.MVar +import System.IO.Unsafe + +{- + We put lighting code in this file because it needs to be + mutually recursive with the trace function, for refraction + and reflection. + -} + +{- +simple_shade :: Rayint -> [Light] -> Solid -> Color -> Color +simple_shade ri lights s bg = + case ri of + (RayHit d p n t) -> + let (Material clr refl refr ior kd shine) = t ri + in cscale clr (vdot n (Vec 0.0 1.0 0.0)) + (RayMiss) -> bg +-} + +debug_norm_shade :: Rayint -> Ray -> Scene -> Int -> Int -> Color +debug_norm_shade ri (Ray o indir) scn recurs debug = + case ri of + RayHit d p (Vec nx ny nz) t -> (Color (abs $ nx/2) (abs $ ny/2) (abs $ nz/2)) + RayMiss -> bground scn + +-- handles diffuse light, shadows, and reflection +-- todo: specular highlights, refraction +shade :: Rayint -> Ray -> Scene -> Int -> Int -> Color +shade ri (Ray o indir) scn recurs !debug = + case ri of + (RayHit d p n t) -> + let (Material clr refl refr ior kd shine) = t ri + s = sld scn + lits = lights scn + direct = foldl' cadd c_black + (map (\ (Light lp lc) -> + let eyedir = vinvert indir + lvec = vsub lp p + llen = vlen lvec + ldir = vscale lvec (1.0/llen) + halfangle = bisect ldir eyedir + ldotn = fmax 0 $ vdot ldir n + -- blinn = fmax 0 ((vdot halfangle n)**(shine*3)) + blinn = fmax 0 $ ((vdot halfangle n) ** shine) * ldotn + blinn_correct = if isNaN blinn then 0 else blinn + -- indotn = fmax 0 $ vdot eyedir n + intensity = 5.0 / (llen*llen) + --intensity = 0.2 + in + --if blinn /= blinn + --then error $ "nan " ++ (show (vdot halfangle n)) ++ " " ++ + -- (show $ shine*3) ++ " " ++ (show blinn) + --else + if not $ shadow s (Ray (vscaleadd p n delta) ldir) (llen-(2*delta)) + then + cadd + -- diffuse + --c_black + (cmul clr $ cscale lc $ ldotn * intensity) + -- blinn/torrance-sparrow highlight (pbrt p 440) + (cscale lc $ blinn_correct * intensity) + -- c_black + else + c_black) lits) + reflect = + if (refl > delta) && (recurs > 0) + then let outdir = Vec.reflect indir n + in cscale (trace scn + (Ray (vscaleadd p outdir delta) outdir) + infinity (recurs-1) ) refl + else c_black + refract = + if (refr > delta) && (recurs > 0) + then c_black + else c_black + in + cadd direct $ cadd reflect refract + + (RayMiss) -> bground scn + + +trace :: Scene -> Ray -> Flt -> Int -> Color + +trace scn ray depth recurs = + let (Scene sld lights cam dtex bgcolor) = scn + in shade (rayint_check sld ray depth dtex) ray scn recurs 0 + + +{- experiments in MVar usage +trace scn ray depth recurs = + let (Scene sld lights cam dtex bgcolor) = scn + in + unsafePerformIO $ + do + debug_start <- (readMVar bihctr) + let ri = rayint sld ray depth dtex + debug_end <- (readMVar bihctr) + let bihhits = debug_end - debug_start + print bihhits + return $ shade ri ray scn recurs bihhits +-} diff --git a/Vec.hs b/Vec.hs new file mode 100644 --- /dev/null +++ b/Vec.hs @@ -0,0 +1,464 @@ +{-# OPTIONS_GHC -fexcess-precision #-} +-- make sure it gets enabled +-- (doesn't seem to help much, though) + +module Vec where + +-- Performance is pretty similar with Floats or Doubles +-- best performance seems to be doubles with -fvia-C +type Flt = Double + +-- This is unnecessary, because haskell has sane mod +-- semantics, unlike ocaml and c. +sane_mod :: Int -> Int -> Int +sane_mod a b = + let modres = mod a b in + if modres < 0 + then modres + b + else modres + +-- maybe this is defined somewhere? +infinity :: Flt +--infinity = 1.0 / 0.0 +infinity = 1000000.0 + +-- convert from degrees to native angle format (radians) +deg :: Flt -> Flt +deg x = (x*3.1415926535897)/180 + +-- convert from radians (noop) +rad :: Flt -> Flt +rad x = x + +-- convert from rotations +rot :: Flt -> Flt +rot x = x*3.1415926535897*2 + +{-abs_sub :: Flt -> Flt -> Flt +abs_sub a b = + if a > 0.0 + then + if b < a + then a-b + else 0.0 + else + if b < (-a) + then a+b + else 0.0 +-} + +clamp :: Flt -> Flt -> Flt -> Flt +clamp !min !x !max + | x < min = min + | x > max = max + | otherwise = x + +-- delta = 0.00001 :: Flt +delta = 0.0001 :: Flt + +-- non-polymorphic versions; this speeds +-- things up in ocaml, not sure about haskell +fmin :: Flt -> Flt -> Flt +fmin !a !b = if a > b then b else a + +fmax :: Flt -> Flt -> Flt +fmax !a !b = if a > b then a else b + +fmin3 :: Flt -> Flt -> Flt -> Flt +fmin3 !a !b !c = if a > b + then if b > c + then c + else b + else if a > c + then c + else a + +fmax3 :: Flt -> Flt -> Flt -> Flt +fmax3 !a !b !c = if a > b + then if a > c + then a + else c + else if b > c + then b + else c + +-- true if a and b are "almost" equal +about_equal :: Flt -> Flt -> Bool +about_equal !a !b = + (abs $ a - b) < (delta*10) + +data Vec = Vec {x, y, z :: !Flt} deriving Show +data Ray = Ray {origin, dir :: !Vec} deriving Show +--data Plane = Plane {norm :: !Vec, offset :: !Flt} deriving Show + +vec !x !y !z = (Vec x y z) +vzero = Vec 0.0 0.0 0.0 + +-- for when we need a unit vector, but we +-- don't care where it points +vunit = vx + +-- unit axis vectors +vx = Vec 1 0 0 +vy = Vec 0 1 0 +vz = Vec 0 0 1 +nvx = Vec (-1) 0 0 +nvy = Vec 0 (-1) 0 +nvz = Vec 0 0 (-1) + +va :: Vec -> Int -> Flt +va (Vec x y z) !i = + case i of + 0 -> x + 1 -> y + 2 -> z + +vset :: Vec -> Int -> Flt -> Vec +vset (Vec x y z) i f = + case i of + 0 -> Vec f y z + 1 -> Vec x f z + 2 -> Vec x y f + +vdot :: Vec -> Vec -> Flt +vdot !(Vec !x1 !y1 !z1) !(Vec !x2 !y2 !z2) = + (x1*x2)+(y1*y2)+(z1*z2) + +vcross :: Vec -> Vec -> Vec +vcross !(Vec !x1 !y1 !z1) !(Vec !x2 !y2 !z2) = + Vec + ((y1 * z2) - (z1 * y2)) + ((z1 * x2) - (x1 * z2)) + ((x1 * y2) - (y1 * x2)) + +vmap :: (Flt -> Flt) -> Vec -> Vec +vmap f !v1 = + Vec (f (x v1)) (f (y v1)) (f (z v1)) + +vmap2 :: (Flt -> Flt -> Flt) -> Vec -> Vec -> Vec +vmap2 f v1 v2 = + Vec (f (x v1) (x v2)) + (f (y v1) (y v2)) + (f (z v1) (z v2)) + +vinvert :: Vec -> Vec +vinvert (Vec !x1 !y1 !z1) = + Vec (-x1) (-y1) (-z1) + +vlensqr :: Vec -> Flt +vlensqr !v1 = vdot v1 v1 + +vlen :: Vec -> Flt +vlen !v1 = sqrt (vdot v1 v1) + +vadd :: Vec -> Vec -> Vec +vadd !(Vec !x1 !y1 !z1) !(Vec !x2 !y2 !z2) = + Vec (x1 + x2) + (y1 + y2) + (z1 + z2) + +vadd3 :: Vec -> Vec -> Vec -> Vec +vadd3 !(Vec !x1 !y1 !z1) !(Vec !x2 !y2 !z2) !(Vec !x3 !y3 !z3) = + Vec (x1 + x2 + x3) + (y1 + y2 + y3) + (z1 + z2 + z3) + +vsub :: Vec -> Vec -> Vec +vsub !(Vec !x1 !y1 !z1) !(Vec !x2 !y2 !z2) = + Vec (x1 - x2) + (y1 - y2) + (z1 - z2) + +vinc :: Vec -> Flt -> Vec +vinc v1 n = + Vec ((x v1) + n) + ((y v1) + n) + ((z v1) + n) + +vdec :: Vec -> Flt -> Vec +vdec v1 n = + Vec ((x v1) - n) + ((y v1) - n) + ((z v1) - n) + +vmax :: Vec -> Vec -> Vec +vmax (Vec x1 y1 z1) (Vec x2 y2 z2) = + Vec (fmax x1 x2) + (fmax y1 y2) + (fmax z1 z2) + +vmin :: Vec -> Vec -> Vec +vmin (Vec x1 y1 z1) (Vec x2 y2 z2) = + Vec (fmin x1 x2) + (fmin y1 y2) + (fmin z1 z2) + +vmaxaxis :: Vec -> Int +vmaxaxis !(Vec !x !y !z) = + if (x > y) + then if (x > z) + then 0 + else 2 + else if (y > z) + then 1 + else 2 + +vscale :: Vec -> Flt -> Vec +vscale !v1 !fac = + Vec ((x v1) * fac) + ((y v1) * fac) + ((z v1) * fac) + +vscaleadd :: Vec -> Vec -> Flt -> Vec +vscaleadd !v1 !v2 fac = + Vec ((x v1) + ((x v2) * fac)) + ((y v1) + ((y v2) * fac)) + ((z v1) + ((z v2) * fac)) + +vnorm :: Vec -> Vec +vnorm (Vec x1 y1 z1) = + let len = 1.0 / (sqrt ((x1*x1)+(y1*y1)+(z1*z1))) in + Vec (x1*len) (y1*len) (z1*len) + +assert_norm :: Vec -> Vec +assert_norm v = + let l = vdot v v + in if l > (1+delta) + then error "vector too long" + else if l < (1-delta) + then error "vector too short" + else v + +bisect :: Vec -> Vec -> Vec +bisect v1 v2 = vnorm (vadd v1 v2) + +vdist :: Vec -> Vec -> Flt +vdist v1 v2 = + let d = vsub v2 v1 in vlen d + +reflect :: Vec -> Vec -> Vec +reflect v norm = + -- vadd v $ vscale norm $ (-2) * (vdot v norm) + vscaleadd v norm $ (-2) * (vdot v norm) + +vrcp :: Vec -> Vec +vrcp (Vec x y z) = + Vec (1/x) (1/y) (1/z) + +-- test for equality +veq :: Vec -> Vec -> Bool +veq (Vec ax ay az) (Vec bx by bz) = + (about_equal ax bx) && (about_equal ay by) && (about_equal az bz) + +-- translate a ray's origin in ray's direction by d amount +ray_move :: Ray -> Flt -> Ray +ray_move (Ray orig dir) d = + (Ray (vscaleadd orig dir d) dir) + +-- find orthogonal vectors +orth :: Vec -> (Vec,Vec) +orth v1 = + if about_equal (vdot v1 v1) 1 + then + let x = (Vec 1 0 0) + y = (Vec 0 1 0) + dvx = vdot v1 x + v2 = if dvx < 0.8 && dvx > (-0.8) -- don't want to cross with a + then vnorm $ vcross v1 x -- vector that's too similar + else vnorm $ vcross v1 y + v3 = vcross v1 v2 + in (v2,v3) + else error $ "orth: unnormalized vector" ++ (show v1) + +-- intersect a ray with a plane +-- defined by a point and a normal +-- (ray need not be normalized) +plane_int :: Ray -> Vec -> Vec -> Vec +plane_int (Ray orig dir) p norm = + let newo = vsub orig p + dist = -(vdot norm newo) / (vdot norm dir) + in vscaleadd orig dir dist + +plane_int_dist :: Ray -> Vec -> Vec -> Flt +plane_int_dist (Ray orig dir) p norm = + let newo = vsub orig p + in -(vdot norm newo) / (vdot norm dir) + +-- find intersection with plane +-- from graphics gems -- an efficient ray-polygon intersection +-- it seems that the ray need not be normalized +-- let plane_intersect ray (n,d) = +-- let t = -.((d +. (vdot n ray.origin)) /. (vdot n ray.dir)) +-- in vadd ray.origin (vscale ray.dir t) + + +-- TRANSFORMATIONS -- + +data Matrix = Matrix !Flt !Flt !Flt !Flt + !Flt !Flt !Flt !Flt + !Flt !Flt !Flt !Flt deriving Show +data Xfm = Xfm !Matrix !Matrix deriving Show + +ident_matrix = (Matrix 1 0 0 0 0 1 0 0 0 0 1 0) +ident_xfm = Xfm ident_matrix ident_matrix + +mat_mult :: Matrix -> Matrix -> Matrix +mat_mult (Matrix a00 a01 a02 a03 a10 a11 a12 a13 a20 a21 a22 a23) + (Matrix b00 b01 b02 b03 b10 b11 b12 b13 b20 b21 b22 b23) = + Matrix + (a00*b00 + a01*b10 + a02*b20) + (a00*b01 + a01*b11 + a02*b21) + (a00*b02 + a01*b12 + a02*b22) + (a00*b03 + a01*b13 + a02*b23 + a03) + + (a10*b00 + a11*b10 + a12*b20) + (a10*b01 + a11*b11 + a12*b21) + (a10*b02 + a11*b12 + a12*b22) + (a10*b03 + a11*b13 + a12*b23 + a13) + + (a20*b00 + a21*b10 + a22*b20) + (a20*b01 + a21*b11 + a22*b21) + (a20*b02 + a21*b12 + a22*b22) + (a20*b03 + a21*b13 + a22*b23 + a23) + +xfm_mult :: Xfm -> Xfm -> Xfm +xfm_mult (Xfm a inva) (Xfm b invb) = + Xfm (mat_mult a b) (mat_mult invb inva) + +-- UTILITY FUNCTIONS -- + +-- If we multiply two transformation matricies, we get +-- a transformation matrix equivalent to applying the +-- second then the first to a vector. + +-- By reversing the list, we apply transforms in expected order. +compose :: [Xfm] -> Xfm +compose xfms = check_xfm $ foldr xfm_mult ident_xfm (reverse xfms) + +check_xfm :: Xfm -> Xfm +check_xfm (Xfm m i) = + let (Matrix m00 m01 m02 m03 + m10 m11 m12 m13 + m20 m21 m22 m23) = mat_mult m i + ae = about_equal + in + if ae m00 1 && ae m01 0 && ae m02 0 && ae m03 0 && + ae m10 0 && ae m11 1 && ae m12 0 && ae m13 0 && + ae m20 0 && ae m21 0 && ae m22 1 && ae m23 0 + then (Xfm m i) + else error $ "corrupt matrix " ++ (show (Xfm m i)) + +-- point is treated as (x y z 1) +xfm_point :: Xfm -> Vec -> Vec +xfm_point (Xfm (Matrix m00 m01 m02 m03 + m10 m11 m12 m13 + m20 m21 m22 m23) inv) + (Vec x y z) = + Vec (m00*x + m01*y + m02*z + m03) + (m10*x + m11*y + m12*z + m13) + (m20*x + m21*y + m22*z + m23) + +invxfm_point :: Xfm -> Vec -> Vec +invxfm_point (Xfm fwd (Matrix i00 i01 i02 i03 + i10 i11 i12 i13 + i20 i21 i22 i23)) + (Vec x y z) = + Vec (i00*x + i01*y + i02*z + i03) + (i10*x + i11*y + i12*z + i13) + (i20*x + i21*y + i22*z + i23) + +-- vector is treated as (x y z 0) +xfm_vec :: Xfm -> Vec -> Vec +xfm_vec (Xfm (Matrix m00 m01 m02 m03 + m10 m11 m12 m13 + m20 m21 m22 m23) inv) + (Vec x y z) = + Vec (m00*x + m01*y + m02*z) + (m10*x + m11*y + m12*z) + (m20*x + m21*y + m22*z) + +invxfm_vec :: Xfm -> Vec -> Vec +invxfm_vec (Xfm fwd (Matrix i00 i01 i02 i03 + i10 i11 i12 i13 + i20 i21 i22 i23)) + (Vec x y z) = + Vec (i00*x + i01*y + i02*z) + (i10*x + i11*y + i12*z) + (i20*x + i21*y + i22*z) + +-- this one is tricky +-- we transform by the inverse transpose +invxfm_norm :: Xfm -> Vec -> Vec +invxfm_norm (Xfm fwd (Matrix i00 i01 i02 i03 + i10 i11 i12 i13 + i20 i21 i22 i23)) + (Vec x y z) = + Vec (i00*x + i10*y + i20*z) + (i01*x + i11*y + i21*z) + (i02*x + i12*y + i22*z) + +xfm_ray :: Xfm -> Ray -> Ray +xfm_ray xfm (Ray orig dir) = + Ray (xfm_point xfm orig) (vnorm (xfm_vec xfm dir)) + +invxfm_ray xfm (Ray orig dir) = + Ray (invxfm_point xfm orig) (vnorm (invxfm_vec xfm dir)) + +-- BASIC TRANSFORMS -- +translate :: Vec -> Xfm +translate (Vec x y z) = + check_xfm $ Xfm (Matrix 1 0 0 x 0 1 0 y 0 0 1 z) + (Matrix 1 0 0 (-x) 0 1 0 (-y) 0 0 1 (-z)) + +scale :: Vec -> Xfm +scale (Vec x y z) = + check_xfm $ Xfm (Matrix x 0 0 0 0 y 0 0 0 0 z 0) + (Matrix (1/x) 0 0 0 0 (1/y) 0 0 0 0 (1/z) 0) + +rotate :: Vec -> Flt -> Xfm +rotate (Vec x y z) angle = + let s = sin angle + c = cos angle + + m00 = ((x*x)+((1-(x*x))*c)) + m01 = (((x*y)*(1-c))-(z*s)) + m02 = ((x*z*(1-c))+(y*s)) + + m10 = (((x*y)*(1-c))+(z*s)) + m11 = ((y*y)+((1-(y*y))*c)) + m12 = ((y*z*(1-c))-(x*s)) + + m20 = ((x*z*(1-c))-(y*s)) + m21 = ((y*z*(1-c))+(x*s)) + m22 = ((z*z)+((1-(z*z))*c)) + in + check_xfm $ Xfm (Matrix m00 m01 m02 0 m10 m11 m12 0 m20 m21 m22 0) + (Matrix m00 m10 m20 0 m01 m11 m21 0 m02 m12 m22 0) + +-- convert canonical coordinates to uvw coordinates +xyz_to_uvw :: Vec -> Vec -> Vec -> Xfm +xyz_to_uvw u v w = + let Vec ux uy uz = u + Vec vx vy vz = v + Vec wx wy wz = w + in if (vdot u u) `about_equal` 1 + then + if (vdot v v) `about_equal` 1 + then + if (vdot w w) `about_equal` 1 + then + if ((vdot u v) `about_equal` 0) && + ((vdot u w) `about_equal` 0) && + ((vdot v w) `about_equal` 0) + then + check_xfm $ Xfm (Matrix ux vx wx 0 uy vy wy 0 uz vz wz 0) + (Matrix ux uy uz 0 vx vy vz 0 wx wy wz 0) + else error "vectors aren't orthogonal" + else error $ "unnormalized w " ++ (show w) + else error $ "unnormalized v " ++ (show v) + else error $ "unnormalized u " ++ (show u) + +uvw_to_xyz :: Vec -> Vec -> Vec -> Xfm +uvw_to_xyz (Vec ux uy uz) (Vec vx vy vz) (Vec wx wy wz) = + check_xfm $ Xfm (Matrix ux uy uz 0 vx vy vz 0 wx wy wz 0) + (Matrix ux vx wx 0 uy vy wy 0 uz vz wz 0) diff --git a/balls3.spd b/balls3.spd new file mode 100644 --- /dev/null +++ b/balls3.spd @@ -0,0 +1,838 @@ +b 0.078 0.361 0.753 +v +from 2.1 1.3 1.7 +at 0 0 0 +up 0 0 1 +angle 45 +hither 0.01 +resolution 512 512 +l 4 3 2 +l 1 -4 4 +l -3 1 5 +f 1 0.75 0.33 0.8 0 100000 0 1 +p 4 +12 12 -0.5 +-12 12 -0.5 +-12 -12 -0.5 +12 -12 -0.5 +f 1 0.9 0.7 0.5 0.5 3.0827 0 1 +s 0 0 0 0.5 +s 0.272166 0.272166 0.544331 0.166667 +s 0.420314 0.420314 0.618405 0.0555556 +s 0.470715 0.470715 0.598245 0.0185185 +s 0.461623 0.409245 0.557924 0.0185185 +s 0.409245 0.461623 0.557924 0.0185185 +s 0.429405 0.481784 0.658726 0.0185185 +s 0.367935 0.472692 0.618405 0.0185185 +s 0.379004 0.431383 0.678886 0.0185185 +s 0.481784 0.429405 0.658726 0.0185185 +s 0.431383 0.379004 0.678886 0.0185185 +s 0.472692 0.367935 0.618405 0.0185185 +s 0.461844 0.304709 0.43322 0.0555556 +s 0.492085 0.33495 0.372739 0.0185185 +s 0.424345 0.305171 0.369341 0.0185185 +s 0.435193 0.368397 0.406378 0.0185185 +s 0.529584 0.334488 0.436618 0.0185185 +s 0.472692 0.367935 0.470257 0.0185185 +s 0.499343 0.304247 0.497099 0.0185185 +s 0.518736 0.271262 0.399581 0.0185185 +s 0.488495 0.241021 0.460062 0.0185185 +s 0.450996 0.241483 0.396183 0.0185185 +s 0.304709 0.461844 0.43322 0.0555556 +s 0.33495 0.492085 0.372739 0.0185185 +s 0.368397 0.435193 0.406378 0.0185185 +s 0.305171 0.424345 0.369341 0.0185185 +s 0.271262 0.518736 0.399581 0.0185185 +s 0.241483 0.450996 0.396183 0.0185185 +s 0.241021 0.488495 0.460062 0.0185185 +s 0.334488 0.529584 0.436618 0.0185185 +s 0.304247 0.499343 0.497099 0.0185185 +s 0.367935 0.472692 0.470257 0.0185185 +s 0.230635 0.38777 0.729516 0.0555556 +s 0.2506 0.446614 0.769837 0.0185185 +s 0.301839 0.407906 0.732914 0.0185185 +s 0.253236 0.449775 0.695877 0.0185185 +s 0.179397 0.426478 0.766439 0.0185185 +s 0.182032 0.429639 0.692479 0.0185185 +s 0.159431 0.367634 0.726118 0.0185185 +s 0.227999 0.384609 0.803476 0.0185185 +s 0.208034 0.325765 0.763155 0.0185185 +s 0.279238 0.345901 0.766553 0.0185185 +s 0.115031 0.4293 0.544331 0.0555556 +s 0.102505 0.502308 0.544331 0.0185185 +s 0.160392 0.474661 0.581368 0.0185185 +s 0.160392 0.474661 0.507294 0.0185185 +s 0.0571437 0.456947 0.507294 0.0185185 +s 0.115031 0.4293 0.470257 0.0185185 +s 0.0696698 0.383939 0.507294 0.0185185 +s 0.0571437 0.456947 0.581368 0.0185185 +s 0.0696698 0.383939 0.581368 0.0185185 +s 0.115031 0.4293 0.618405 0.0185185 +s 0.082487 0.239622 0.655442 0.0555556 +s 0.0438957 0.258053 0.715923 0.0185185 +s 0.117687 0.252557 0.719322 0.0185185 +s 0.0863845 0.308551 0.682285 0.0185185 +s 0.00869528 0.245118 0.652044 0.0185185 +s 0.0511841 0.295616 0.618405 0.0185185 +s 0.0472866 0.226687 0.591563 0.0185185 +s 0.0399982 0.189123 0.689081 0.0185185 +s 0.0785895 0.170692 0.6286 0.0185185 +s 0.11379 0.183628 0.692479 0.0185185 +s 0.38777 0.230635 0.729516 0.0555556 +s 0.446614 0.2506 0.769837 0.0185185 +s 0.449775 0.253236 0.695877 0.0185185 +s 0.407906 0.301839 0.732914 0.0185185 +s 0.384609 0.227999 0.803476 0.0185185 +s 0.345901 0.279238 0.766553 0.0185185 +s 0.325765 0.208034 0.763155 0.0185185 +s 0.426478 0.179397 0.766439 0.0185185 +s 0.367634 0.159431 0.726118 0.0185185 +s 0.429639 0.182032 0.692479 0.0185185 +s 0.239622 0.082487 0.655442 0.0555556 +s 0.258053 0.0438957 0.715923 0.0185185 +s 0.308551 0.0863845 0.682285 0.0185185 +s 0.252557 0.117687 0.719322 0.0185185 +s 0.189123 0.0399982 0.689081 0.0185185 +s 0.183628 0.11379 0.692479 0.0185185 +s 0.170692 0.0785895 0.6286 0.0185185 +s 0.245118 0.00869528 0.652044 0.0185185 +s 0.226687 0.0472866 0.591563 0.0185185 +s 0.295616 0.0511841 0.618405 0.0185185 +s 0.4293 0.115031 0.544331 0.0555556 +s 0.502308 0.102505 0.544331 0.0185185 +s 0.474661 0.160392 0.507294 0.0185185 +s 0.474661 0.160392 0.581368 0.0185185 +s 0.456947 0.0571437 0.581368 0.0185185 +s 0.4293 0.115031 0.618405 0.0185185 +s 0.383939 0.0696698 0.581368 0.0185185 +s 0.456947 0.0571437 0.507294 0.0185185 +s 0.383939 0.0696698 0.507294 0.0185185 +s 0.4293 0.115031 0.470257 0.0185185 +s 0.643951 0.172546 1.11022e-16 0.166667 +s 0.802608 0.281471 -0.111111 0.0555556 +s 0.824035 0.30566 -0.177765 0.0185185 +s 0.787796 0.241352 -0.171592 0.0185185 +s 0.752156 0.305221 -0.15987 0.0185185 +s 0.838847 0.34578 -0.117284 0.0185185 +s 0.766968 0.345341 -0.099389 0.0185185 +s 0.817421 0.321591 -0.0506299 0.0185185 +s 0.874487 0.28191 -0.129006 0.0185185 +s 0.853061 0.257721 -0.062352 0.0185185 +s 0.838248 0.217602 -0.122833 0.0185185 +s 0.643951 0.172546 -0.222222 0.0555556 +s 0.61371 0.202787 -0.282703 0.0185185 +s 0.5724 0.191718 -0.222222 0.0185185 +s 0.624779 0.244096 -0.222222 0.0185185 +s 0.68526 0.183615 -0.282703 0.0185185 +s 0.696329 0.224924 -0.222222 0.0185185 +s 0.715501 0.153374 -0.222222 0.0185185 +s 0.632882 0.131237 -0.282703 0.0185185 +s 0.663122 0.100996 -0.222222 0.0185185 +s 0.591572 0.120168 -0.222222 0.0185185 +s 0.594141 0.358439 -0.111111 0.0555556 +s 0.619127 0.408291 -0.15987 0.0185185 +s 0.665691 0.37761 -0.111111 0.0185185 +s 0.638217 0.337042 -0.166667 0.0185185 +s 0.547576 0.389119 -0.15987 0.0185185 +s 0.566667 0.317871 -0.166667 0.0185185 +s 0.522591 0.339267 -0.111111 0.0185185 +s 0.57505 0.429687 -0.104315 0.0185185 +s 0.550064 0.379835 -0.0555556 0.0185185 +s 0.621614 0.399007 -0.0555556 0.0185185 +s 0.802608 0.281471 0.111111 0.0555556 +s 0.858698 0.329459 0.104938 0.0185185 +s 0.845371 0.280879 0.0506299 0.0185185 +s 0.798572 0.337088 0.062352 0.0185185 +s 0.815936 0.330051 0.165419 0.0185185 +s 0.75581 0.33768 0.122833 0.0185185 +s 0.759846 0.282063 0.171592 0.0185185 +s 0.862735 0.273842 0.153697 0.0185185 +s 0.806645 0.225855 0.15987 0.0185185 +s 0.849407 0.225263 0.099389 0.0185185 +s 0.594141 0.358439 0.111111 0.0555556 +s 0.613592 0.428945 0.122833 0.0185185 +s 0.665691 0.37761 0.111111 0.0185185 +s 0.621614 0.399007 0.0555556 0.0185185 +s 0.542042 0.409774 0.122833 0.0185185 +s 0.550064 0.379835 0.0555556 0.0185185 +s 0.522591 0.339267 0.111111 0.0185185 +s 0.586119 0.388377 0.178389 0.0185185 +s 0.566667 0.317871 0.166667 0.0185185 +s 0.638217 0.337042 0.166667 0.0185185 +s 0.643951 0.172546 0.222222 0.0555556 +s 0.674191 0.202787 0.282703 0.0185185 +s 0.715501 0.191718 0.222222 0.0185185 +s 0.663122 0.244096 0.222222 0.0185185 +s 0.602641 0.183615 0.282703 0.0185185 +s 0.591572 0.224924 0.222222 0.0185185 +s 0.5724 0.153374 0.222222 0.0185185 +s 0.655019 0.131237 0.282703 0.0185185 +s 0.624779 0.100996 0.222222 0.0185185 +s 0.696329 0.120168 0.222222 0.0185185 +s 0.852418 0.0955788 1.89979e-16 0.0555556 +s 0.922609 0.11107 -0.0178949 0.0185185 +s 0.867231 0.135698 -0.0604812 0.0185185 +s 0.877966 0.164775 0.00679642 0.0185185 +s 0.907797 0.0709499 0.0425863 0.0185185 +s 0.863153 0.124655 0.0672777 0.0185185 +s 0.837606 0.0554592 0.0604812 0.0185185 +s 0.897062 0.0418734 -0.0246914 0.0185185 +s 0.826871 0.0263827 -0.00679642 0.0185185 +s 0.841683 0.0665023 -0.0672777 0.0185185 +s 0.69376 -0.0133465 0.111111 0.0555556 +s 0.740325 -0.0440268 0.15987 0.0185185 +s 0.76531 0.0058253 0.111111 0.0185185 +s 0.721234 0.0272215 0.166667 0.0185185 +s 0.668775 -0.0631985 0.15987 0.0185185 +s 0.649684 0.00804971 0.166667 0.0185185 +s 0.62221 -0.0325183 0.111111 0.0185185 +s 0.712851 -0.0845947 0.104315 0.0185185 +s 0.666287 -0.0539145 0.0555556 0.0185185 +s 0.737837 -0.0347427 0.0555556 0.0185185 +s 0.69376 -0.0133465 -0.111111 0.0555556 +s 0.745859 -0.0646815 -0.122833 0.0185185 +s 0.76531 0.0058253 -0.111111 0.0185185 +s 0.737837 -0.0347427 -0.0555556 0.0185185 +s 0.674309 -0.0838533 -0.122833 0.0185185 +s 0.666287 -0.0539145 -0.0555556 0.0185185 +s 0.62221 -0.0325183 -0.111111 0.0185185 +s 0.701782 -0.0432853 -0.178389 0.0185185 +s 0.649684 0.00804971 -0.166667 0.0185185 +s 0.721234 0.0272215 -0.166667 0.0185185 +s 0.172546 0.643951 1.11022e-16 0.166667 +s 0.281471 0.802608 -0.111111 0.0555556 +s 0.30566 0.824035 -0.177765 0.0185185 +s 0.305221 0.752156 -0.15987 0.0185185 +s 0.241352 0.787796 -0.171592 0.0185185 +s 0.28191 0.874487 -0.129006 0.0185185 +s 0.217602 0.838248 -0.122833 0.0185185 +s 0.257721 0.853061 -0.062352 0.0185185 +s 0.34578 0.838847 -0.117284 0.0185185 +s 0.321591 0.817421 -0.0506299 0.0185185 +s 0.345341 0.766968 -0.099389 0.0185185 +s 0.358439 0.594141 -0.111111 0.0555556 +s 0.408291 0.619127 -0.15987 0.0185185 +s 0.337042 0.638217 -0.166667 0.0185185 +s 0.37761 0.665691 -0.111111 0.0185185 +s 0.429687 0.57505 -0.104315 0.0185185 +s 0.399007 0.621614 -0.0555556 0.0185185 +s 0.379835 0.550064 -0.0555556 0.0185185 +s 0.389119 0.547576 -0.15987 0.0185185 +s 0.339267 0.522591 -0.111111 0.0185185 +s 0.317871 0.566667 -0.166667 0.0185185 +s 0.172546 0.643951 -0.222222 0.0555556 +s 0.142305 0.674191 -0.282703 0.0185185 +s 0.100996 0.663122 -0.222222 0.0185185 +s 0.153374 0.715501 -0.222222 0.0185185 +s 0.213855 0.655019 -0.282703 0.0185185 +s 0.224924 0.696329 -0.222222 0.0185185 +s 0.244096 0.624779 -0.222222 0.0185185 +s 0.161477 0.602641 -0.282703 0.0185185 +s 0.191718 0.5724 -0.222222 0.0185185 +s 0.120168 0.591572 -0.222222 0.0185185 +s 0.0955788 0.852418 1.03629e-16 0.0555556 +s 0.11107 0.922609 -0.0178949 0.0185185 +s 0.164775 0.877966 0.00679642 0.0185185 +s 0.135698 0.867231 -0.0604812 0.0185185 +s 0.0418734 0.897062 -0.0246914 0.0185185 +s 0.0665023 0.841683 -0.0672777 0.0185185 +s 0.0263827 0.826871 -0.00679642 0.0185185 +s 0.0709499 0.907797 0.0425863 0.0185185 +s 0.0554592 0.837606 0.0604812 0.0185185 +s 0.124655 0.863153 0.0672777 0.0185185 +s -0.0133465 0.69376 -0.111111 0.0555556 +s -0.0646815 0.745859 -0.122833 0.0185185 +s -0.0347427 0.737837 -0.0555556 0.0185185 +s 0.0058253 0.76531 -0.111111 0.0185185 +s -0.0432853 0.701782 -0.178389 0.0185185 +s 0.0272215 0.721234 -0.166667 0.0185185 +s 0.00804971 0.649684 -0.166667 0.0185185 +s -0.0838533 0.674309 -0.122833 0.0185185 +s -0.0325183 0.62221 -0.111111 0.0185185 +s -0.0539145 0.666287 -0.0555556 0.0185185 +s -0.0133465 0.69376 0.111111 0.0555556 +s -0.0440268 0.740325 0.15987 0.0185185 +s 0.0272215 0.721234 0.166667 0.0185185 +s 0.0058253 0.76531 0.111111 0.0185185 +s -0.0845947 0.712851 0.104315 0.0185185 +s -0.0347427 0.737837 0.0555556 0.0185185 +s -0.0539145 0.666287 0.0555556 0.0185185 +s -0.0631985 0.668775 0.15987 0.0185185 +s -0.0325183 0.62221 0.111111 0.0185185 +s 0.00804971 0.649684 0.166667 0.0185185 +s 0.281471 0.802608 0.111111 0.0555556 +s 0.329459 0.858698 0.104938 0.0185185 +s 0.337088 0.798572 0.062352 0.0185185 +s 0.280879 0.845371 0.0506299 0.0185185 +s 0.273842 0.862735 0.153697 0.0185185 +s 0.225263 0.849407 0.099389 0.0185185 +s 0.225855 0.806645 0.15987 0.0185185 +s 0.330051 0.815936 0.165419 0.0185185 +s 0.282063 0.759846 0.171592 0.0185185 +s 0.33768 0.75581 0.122833 0.0185185 +s 0.172546 0.643951 0.222222 0.0555556 +s 0.202787 0.674191 0.282703 0.0185185 +s 0.244096 0.663122 0.222222 0.0185185 +s 0.191718 0.715501 0.222222 0.0185185 +s 0.131237 0.655019 0.282703 0.0185185 +s 0.120168 0.696329 0.222222 0.0185185 +s 0.100996 0.624779 0.222222 0.0185185 +s 0.183615 0.602641 0.282703 0.0185185 +s 0.153374 0.5724 0.222222 0.0185185 +s 0.224924 0.591572 0.222222 0.0185185 +s 0.358439 0.594141 0.111111 0.0555556 +s 0.428945 0.613592 0.122833 0.0185185 +s 0.399007 0.621614 0.0555556 0.0185185 +s 0.37761 0.665691 0.111111 0.0185185 +s 0.388377 0.586119 0.178389 0.0185185 +s 0.337042 0.638217 0.166667 0.0185185 +s 0.317871 0.566667 0.166667 0.0185185 +s 0.409774 0.542042 0.122833 0.0185185 +s 0.339267 0.522591 0.111111 0.0185185 +s 0.379835 0.550064 0.0555556 0.0185185 +s -0.371785 0.0996195 0.544331 0.166667 +s -0.393621 0.220501 0.729516 0.0555556 +s -0.368601 0.279642 0.766439 0.0185185 +s -0.321889 0.238665 0.726118 0.0185185 +s -0.372464 0.281062 0.692479 0.0185185 +s -0.440333 0.261479 0.769837 0.0185185 +s -0.444196 0.262898 0.695877 0.0185185 +s -0.465353 0.202338 0.732914 0.0185185 +s -0.389758 0.219082 0.803476 0.0185185 +s -0.414778 0.15994 0.766553 0.0185185 +s -0.343046 0.178104 0.763155 0.0185185 +s -0.191247 0.166275 0.655442 0.0555556 +s -0.130089 0.20793 0.652044 0.0185185 +s -0.154295 0.172673 0.591563 0.0185185 +s -0.192135 0.230419 0.618405 0.0185185 +s -0.167041 0.201532 0.715923 0.0185185 +s -0.229087 0.224021 0.682285 0.0185185 +s -0.228199 0.159877 0.719322 0.0185185 +s -0.129201 0.143787 0.689081 0.0185185 +s -0.190359 0.102131 0.692479 0.0185185 +s -0.153407 0.108529 0.6286 0.0185185 +s -0.31427 0.31427 0.544331 0.0555556 +s -0.277961 0.367156 0.507294 0.0185185 +s -0.252306 0.297666 0.507294 0.0185185 +s -0.31427 0.31427 0.470257 0.0185185 +s -0.339925 0.383759 0.544331 0.0185185 +s -0.376234 0.330873 0.507294 0.0185185 +s -0.376234 0.330873 0.581368 0.0185185 +s -0.277961 0.367156 0.581368 0.0185185 +s -0.31427 0.31427 0.618405 0.0185185 +s -0.252306 0.297666 0.581368 0.0185185 +s -0.574159 0.153845 0.618405 0.0555556 +s -0.612768 0.202534 0.658726 0.0185185 +s -0.543919 0.184086 0.678886 0.0185185 +s -0.554987 0.225396 0.618405 0.0185185 +s -0.643008 0.172294 0.598245 0.0185185 +s -0.585228 0.195155 0.557924 0.0185185 +s -0.6044 0.123605 0.557924 0.0185185 +s -0.631939 0.130984 0.658726 0.0185185 +s -0.593331 0.0822954 0.618405 0.0185185 +s -0.56309 0.112536 0.678886 0.0185185 +s -0.494808 0.247614 0.43322 0.0555556 +s -0.494287 0.313607 0.399581 0.0185185 +s -0.452978 0.302539 0.460062 0.0185185 +s -0.434629 0.269833 0.396183 0.0185185 +s -0.536117 0.258683 0.372739 0.0185185 +s -0.476459 0.214908 0.369341 0.0185185 +s -0.536638 0.19269 0.406378 0.0185185 +s -0.554467 0.291389 0.436618 0.0185185 +s -0.554987 0.225396 0.470257 0.0185185 +s -0.513157 0.28032 0.497099 0.0185185 +s -0.552323 0.0329639 0.43322 0.0555556 +s -0.625877 0.0248832 0.436618 0.0185185 +s -0.584567 0.0138144 0.497099 0.0185185 +s -0.593331 0.0822954 0.470257 0.0185185 +s -0.593633 0.0440327 0.372739 0.0185185 +s -0.561087 0.101445 0.406378 0.0185185 +s -0.520079 0.0521134 0.369341 0.0185185 +s -0.584869 -0.0244483 0.399581 0.0185185 +s -0.511316 -0.0163676 0.396183 0.0185185 +s -0.54356 -0.0355172 0.460062 0.0185185 +s -0.451136 0.0058509 0.729516 0.0555556 +s -0.447081 0.0051487 0.803476 0.0185185 +s -0.386138 0.0172804 0.763155 0.0185185 +s -0.439178 0.0688765 0.766553 0.0185185 +s -0.512079 -0.00628079 0.769837 0.0185185 +s -0.504176 0.0574471 0.732914 0.0185185 +s -0.516134 -0.00557859 0.695877 0.0185185 +s -0.459039 -0.0578769 0.766439 0.0185185 +s -0.463094 -0.0571747 0.692479 0.0185185 +s -0.398096 -0.0457452 0.726118 0.0185185 +s -0.4293 -0.115031 0.544331 0.0555556 +s -0.424299 -0.178985 0.581368 0.0185185 +s -0.367336 -0.131634 0.581368 0.0185185 +s -0.4293 -0.115031 0.618405 0.0185185 +s -0.486264 -0.162382 0.544331 0.0185185 +s -0.491265 -0.0984274 0.581368 0.0185185 +s -0.491265 -0.0984274 0.507294 0.0185185 +s -0.424299 -0.178985 0.507294 0.0185185 +s -0.4293 -0.115031 0.470257 0.0185185 +s -0.367336 -0.131634 0.507294 0.0185185 +s -0.248762 -0.0483751 0.655442 0.0555556 +s -0.183785 -0.0599222 0.689081 0.0185185 +s -0.187119 -0.0172857 0.6286 0.0185185 +s -0.215921 0.00673113 0.692479 0.0185185 +s -0.245428 -0.0910116 0.715923 0.0185185 +s -0.277564 -0.0243582 0.719322 0.0185185 +s -0.310405 -0.0794645 0.682285 0.0185185 +s -0.216626 -0.115028 0.652044 0.0185185 +s -0.281603 -0.103481 0.618405 0.0185185 +s -0.21996 -0.0723919 0.591563 0.0185185 +s -0.471405 0.471405 1.11022e-16 0.166667 +s -0.508983 0.690426 8.51251e-17 0.0555556 +s -0.484794 0.755941 -0.0246914 0.0185185 +s -0.436283 0.7029 -0.00679642 0.0185185 +s -0.478434 0.695668 -0.0672777 0.0185185 +s -0.557494 0.743468 -0.0178949 0.0185185 +s -0.551134 0.683194 -0.0604812 0.0185185 +s -0.581682 0.677953 0.00679642 0.0185185 +s -0.515343 0.7507 0.0425863 0.0185185 +s -0.539531 0.685185 0.0672777 0.0185185 +s -0.466832 0.697658 0.0604812 0.0185185 +s -0.335322 0.607487 0.111111 0.0555556 +s -0.283164 0.659645 0.104315 0.0185185 +s -0.286452 0.603979 0.0555556 0.0185185 +s -0.33883 0.656357 0.0555556 0.0185185 +s -0.332034 0.663153 0.15987 0.0185185 +s -0.3877 0.659866 0.111111 0.0185185 +s -0.384191 0.610996 0.166667 0.0185185 +s -0.279656 0.610775 0.15987 0.0185185 +s -0.331813 0.558618 0.166667 0.0185185 +s -0.282943 0.555109 0.111111 0.0185185 +s -0.335322 0.607487 -0.111111 0.0555556 +s -0.313405 0.629404 -0.178389 0.0185185 +s -0.331813 0.558618 -0.166667 0.0185185 +s -0.384191 0.610996 -0.166667 0.0185185 +s -0.316914 0.678274 -0.122833 0.0185185 +s -0.3877 0.659866 -0.111111 0.0185185 +s -0.33883 0.656357 -0.0555556 0.0185185 +s -0.264535 0.625895 -0.122833 0.0185185 +s -0.286452 0.603979 -0.0555556 0.0185185 +s -0.282943 0.555109 -0.111111 0.0185185 +s -0.645066 0.554344 -0.111111 0.0555556 +s -0.681385 0.616373 -0.129006 0.0185185 +s -0.649723 0.609912 -0.062352 0.0185185 +s -0.607573 0.617144 -0.122833 0.0185185 +s -0.676727 0.560805 -0.177765 0.0185185 +s -0.602915 0.561576 -0.171592 0.0185185 +s -0.640408 0.498776 -0.15987 0.0185185 +s -0.718878 0.553573 -0.117284 0.0185185 +s -0.682558 0.491544 -0.099389 0.0185185 +s -0.687216 0.547112 -0.0506299 0.0185185 +s -0.471405 0.471405 -0.222222 0.0555556 +s -0.501645 0.501645 -0.282703 0.0185185 +s -0.542955 0.490576 -0.222222 0.0185185 +s -0.490576 0.542955 -0.222222 0.0185185 +s -0.430095 0.482473 -0.282703 0.0185185 +s -0.419026 0.523783 -0.222222 0.0185185 +s -0.399854 0.452233 -0.222222 0.0185185 +s -0.482473 0.430095 -0.282703 0.0185185 +s -0.452233 0.399854 -0.222222 0.0185185 +s -0.523783 0.419026 -0.222222 0.0185185 +s -0.607487 0.335322 -0.111111 0.0555556 +s -0.659645 0.283164 -0.104315 0.0185185 +s -0.603979 0.286452 -0.0555556 0.0185185 +s -0.656357 0.33883 -0.0555556 0.0185185 +s -0.663153 0.332034 -0.15987 0.0185185 +s -0.659866 0.3877 -0.111111 0.0185185 +s -0.610996 0.384191 -0.166667 0.0185185 +s -0.610775 0.279656 -0.15987 0.0185185 +s -0.558618 0.331813 -0.166667 0.0185185 +s -0.555109 0.282943 -0.111111 0.0185185 +s -0.645066 0.554344 0.111111 0.0555556 +s -0.668521 0.610229 0.153697 0.0185185 +s -0.598918 0.585648 0.15987 0.0185185 +s -0.619787 0.622977 0.099389 0.0185185 +s -0.714669 0.578925 0.104938 0.0185185 +s -0.665934 0.591673 0.0506299 0.0185185 +s -0.691213 0.52304 0.062352 0.0185185 +s -0.6938 0.541596 0.165419 0.0185185 +s -0.670344 0.48571 0.122833 0.0185185 +s -0.624197 0.517014 0.171592 0.0185185 +s -0.607487 0.335322 0.111111 0.0555556 +s -0.629404 0.313405 0.178389 0.0185185 +s -0.558618 0.331813 0.166667 0.0185185 +s -0.610996 0.384191 0.166667 0.0185185 +s -0.678274 0.316914 0.122833 0.0185185 +s -0.659866 0.3877 0.111111 0.0185185 +s -0.656357 0.33883 0.0555556 0.0185185 +s -0.625895 0.264535 0.122833 0.0185185 +s -0.603979 0.286452 0.0555556 0.0185185 +s -0.555109 0.282943 0.111111 0.0185185 +s -0.471405 0.471405 0.222222 0.0555556 +s -0.441164 0.501645 0.282703 0.0185185 +s -0.399854 0.490576 0.222222 0.0185185 +s -0.452233 0.542955 0.222222 0.0185185 +s -0.512714 0.482473 0.282703 0.0185185 +s -0.523783 0.523783 0.222222 0.0185185 +s -0.542955 0.452233 0.222222 0.0185185 +s -0.460336 0.430095 0.282703 0.0185185 +s -0.490576 0.399854 0.222222 0.0185185 +s -0.419026 0.419026 0.222222 0.0185185 +s -0.643951 -0.172546 1.11022e-16 0.166667 +s -0.835815 -0.157543 0.111111 0.0555556 +s -0.871646 -0.122136 0.165419 0.0185185 +s -0.799077 -0.135649 0.171592 0.0185185 +s -0.82339 -0.0854653 0.122833 0.0185185 +s -0.908384 -0.14403 0.104938 0.0185185 +s -0.860128 -0.107359 0.062352 0.0185185 +s -0.872552 -0.179437 0.0506299 0.0185185 +s -0.884071 -0.194213 0.153697 0.0185185 +s -0.84824 -0.229621 0.099389 0.0185185 +s -0.811502 -0.207727 0.15987 0.0185185 +s -0.643951 -0.172546 0.222222 0.0555556 +s -0.61371 -0.142305 0.282703 0.0185185 +s -0.5724 -0.153374 0.222222 0.0185185 +s -0.624779 -0.100996 0.222222 0.0185185 +s -0.68526 -0.161477 0.282703 0.0185185 +s -0.696329 -0.120168 0.222222 0.0185185 +s -0.715501 -0.191718 0.222222 0.0185185 +s -0.632882 -0.213855 0.282703 0.0185185 +s -0.663122 -0.244096 0.222222 0.0185185 +s -0.591572 -0.224924 0.222222 0.0185185 +s -0.69376 0.0133465 0.111111 0.0555556 +s -0.674309 0.0838533 0.122833 0.0185185 +s -0.62221 0.0325183 0.111111 0.0185185 +s -0.666287 0.0539145 0.0555556 0.0185185 +s -0.745859 0.0646815 0.122833 0.0185185 +s -0.737837 0.0347427 0.0555556 0.0185185 +s -0.76531 -0.0058253 0.111111 0.0185185 +s -0.701782 0.0432853 0.178389 0.0185185 +s -0.721234 -0.0272215 0.166667 0.0185185 +s -0.649684 -0.00804971 0.166667 0.0185185 +s -0.835815 -0.157543 -0.111111 0.0555556 +s -0.899353 -0.119969 -0.117284 0.0185185 +s -0.868703 -0.130205 -0.0506299 0.0185185 +s -0.836885 -0.0844101 -0.099389 0.0185185 +s -0.866465 -0.147308 -0.177765 0.0185185 +s -0.803997 -0.111748 -0.15987 0.0185185 +s -0.802927 -0.184881 -0.171592 0.0185185 +s -0.898283 -0.193102 -0.129006 0.0185185 +s -0.834745 -0.230676 -0.122833 0.0185185 +s -0.867633 -0.203337 -0.062352 0.0185185 +s -0.69376 0.0133465 -0.111111 0.0555556 +s -0.668775 0.0631985 -0.15987 0.0185185 +s -0.62221 0.0325183 -0.111111 0.0185185 +s -0.649684 -0.00804971 -0.166667 0.0185185 +s -0.740325 0.0440268 -0.15987 0.0185185 +s -0.721234 -0.0272215 -0.166667 0.0185185 +s -0.76531 -0.0058253 -0.111111 0.0185185 +s -0.712851 0.0845947 -0.104315 0.0185185 +s -0.737837 0.0347427 -0.0555556 0.0185185 +s -0.666287 0.0539145 -0.0555556 0.0185185 +s -0.643951 -0.172546 -0.222222 0.0555556 +s -0.674191 -0.202787 -0.282703 0.0185185 +s -0.715501 -0.191718 -0.222222 0.0185185 +s -0.663122 -0.244096 -0.222222 0.0185185 +s -0.602641 -0.183615 -0.282703 0.0185185 +s -0.591572 -0.224924 -0.222222 0.0185185 +s -0.5724 -0.153374 -0.222222 0.0185185 +s -0.655019 -0.131237 -0.282703 0.0185185 +s -0.624779 -0.100996 -0.222222 0.0185185 +s -0.696329 -0.120168 -0.222222 0.0185185 +s -0.786005 -0.343435 7.89572e-17 0.0555556 +s -0.82165 -0.392454 0.0425863 0.0185185 +s -0.753118 -0.370774 0.0604812 0.0185185 +s -0.80984 -0.323622 0.0672777 0.0185185 +s -0.854538 -0.365116 -0.0178949 0.0185185 +s -0.842728 -0.296284 0.00679642 0.0185185 +s -0.818893 -0.316097 -0.0604812 0.0185185 +s -0.797815 -0.412267 -0.0246914 0.0185185 +s -0.76217 -0.363249 -0.0672777 0.0185185 +s -0.729282 -0.390587 -0.00679642 0.0185185 +s -0.594141 -0.358439 -0.111111 0.0555556 +s -0.542042 -0.409774 -0.122833 0.0185185 +s -0.522591 -0.339267 -0.111111 0.0185185 +s -0.550064 -0.379835 -0.0555556 0.0185185 +s -0.613592 -0.428945 -0.122833 0.0185185 +s -0.621614 -0.399007 -0.0555556 0.0185185 +s -0.665691 -0.37761 -0.111111 0.0185185 +s -0.586119 -0.388377 -0.178389 0.0185185 +s -0.638217 -0.337042 -0.166667 0.0185185 +s -0.566667 -0.317871 -0.166667 0.0185185 +s -0.594141 -0.358439 0.111111 0.0555556 +s -0.547576 -0.389119 0.15987 0.0185185 +s -0.522591 -0.339267 0.111111 0.0185185 +s -0.566667 -0.317871 0.166667 0.0185185 +s -0.619127 -0.408291 0.15987 0.0185185 +s -0.638217 -0.337042 0.166667 0.0185185 +s -0.665691 -0.37761 0.111111 0.0185185 +s -0.57505 -0.429687 0.104315 0.0185185 +s -0.621614 -0.399007 0.0555556 0.0185185 +s -0.550064 -0.379835 0.0555556 0.0185185 +s 0.0996195 -0.371785 0.544331 0.166667 +s 0.220501 -0.393621 0.729516 0.0555556 +s 0.279642 -0.368601 0.766439 0.0185185 +s 0.281062 -0.372464 0.692479 0.0185185 +s 0.238665 -0.321889 0.726118 0.0185185 +s 0.219082 -0.389758 0.803476 0.0185185 +s 0.178104 -0.343046 0.763155 0.0185185 +s 0.15994 -0.414778 0.766553 0.0185185 +s 0.261479 -0.440333 0.769837 0.0185185 +s 0.202338 -0.465353 0.732914 0.0185185 +s 0.262898 -0.444196 0.695877 0.0185185 +s 0.31427 -0.31427 0.544331 0.0555556 +s 0.367156 -0.277961 0.507294 0.0185185 +s 0.31427 -0.31427 0.470257 0.0185185 +s 0.297666 -0.252306 0.507294 0.0185185 +s 0.367156 -0.277961 0.581368 0.0185185 +s 0.297666 -0.252306 0.581368 0.0185185 +s 0.31427 -0.31427 0.618405 0.0185185 +s 0.383759 -0.339925 0.544331 0.0185185 +s 0.330873 -0.376234 0.581368 0.0185185 +s 0.330873 -0.376234 0.507294 0.0185185 +s 0.166275 -0.191247 0.655442 0.0555556 +s 0.20793 -0.130089 0.652044 0.0185185 +s 0.230419 -0.192135 0.618405 0.0185185 +s 0.172673 -0.154295 0.591563 0.0185185 +s 0.143787 -0.129201 0.689081 0.0185185 +s 0.108529 -0.153407 0.6286 0.0185185 +s 0.102131 -0.190359 0.692479 0.0185185 +s 0.201532 -0.167041 0.715923 0.0185185 +s 0.159877 -0.228199 0.719322 0.0185185 +s 0.224021 -0.229087 0.682285 0.0185185 +s 0.0058509 -0.451136 0.729516 0.0555556 +s 0.0051487 -0.447081 0.803476 0.0185185 +s 0.0688765 -0.439178 0.766553 0.0185185 +s 0.0172804 -0.386138 0.763155 0.0185185 +s -0.0578769 -0.459039 0.766439 0.0185185 +s -0.0457452 -0.398096 0.726118 0.0185185 +s -0.0571747 -0.463094 0.692479 0.0185185 +s -0.00628079 -0.512079 0.769837 0.0185185 +s -0.00557859 -0.516134 0.695877 0.0185185 +s 0.0574471 -0.504176 0.732914 0.0185185 +s -0.0483751 -0.248762 0.655442 0.0555556 +s -0.0599222 -0.183785 0.689081 0.0185185 +s 0.00673113 -0.215921 0.692479 0.0185185 +s -0.0172857 -0.187119 0.6286 0.0185185 +s -0.115028 -0.216626 0.652044 0.0185185 +s -0.0723919 -0.21996 0.591563 0.0185185 +s -0.103481 -0.281603 0.618405 0.0185185 +s -0.0910116 -0.245428 0.715923 0.0185185 +s -0.0794645 -0.310405 0.682285 0.0185185 +s -0.0243582 -0.277564 0.719322 0.0185185 +s -0.115031 -0.4293 0.544331 0.0555556 +s -0.178985 -0.424299 0.581368 0.0185185 +s -0.115031 -0.4293 0.618405 0.0185185 +s -0.131634 -0.367336 0.581368 0.0185185 +s -0.178985 -0.424299 0.507294 0.0185185 +s -0.131634 -0.367336 0.507294 0.0185185 +s -0.115031 -0.4293 0.470257 0.0185185 +s -0.162382 -0.486264 0.544331 0.0185185 +s -0.0984274 -0.491265 0.507294 0.0185185 +s -0.0984274 -0.491265 0.581368 0.0185185 +s 0.153845 -0.574159 0.618405 0.0555556 +s 0.202534 -0.612768 0.658726 0.0185185 +s 0.225396 -0.554987 0.618405 0.0185185 +s 0.184086 -0.543919 0.678886 0.0185185 +s 0.130984 -0.631939 0.658726 0.0185185 +s 0.112536 -0.56309 0.678886 0.0185185 +s 0.0822954 -0.593331 0.618405 0.0185185 +s 0.172294 -0.643008 0.598245 0.0185185 +s 0.123605 -0.6044 0.557924 0.0185185 +s 0.195155 -0.585228 0.557924 0.0185185 +s 0.0329639 -0.552323 0.43322 0.0555556 +s 0.0248832 -0.625877 0.436618 0.0185185 +s 0.0822954 -0.593331 0.470257 0.0185185 +s 0.0138144 -0.584567 0.497099 0.0185185 +s -0.0244483 -0.584869 0.399581 0.0185185 +s -0.0355172 -0.54356 0.460062 0.0185185 +s -0.0163676 -0.511316 0.396183 0.0185185 +s 0.0440327 -0.593633 0.372739 0.0185185 +s 0.0521134 -0.520079 0.369341 0.0185185 +s 0.101445 -0.561087 0.406378 0.0185185 +s 0.247614 -0.494808 0.43322 0.0555556 +s 0.313607 -0.494287 0.399581 0.0185185 +s 0.269833 -0.434629 0.396183 0.0185185 +s 0.302539 -0.452978 0.460062 0.0185185 +s 0.291389 -0.554467 0.436618 0.0185185 +s 0.28032 -0.513157 0.497099 0.0185185 +s 0.225396 -0.554987 0.470257 0.0185185 +s 0.258683 -0.536117 0.372739 0.0185185 +s 0.19269 -0.536638 0.406378 0.0185185 +s 0.214908 -0.476459 0.369341 0.0185185 +s -0.172546 -0.643951 1.11022e-16 0.166667 +s -0.157543 -0.835815 0.111111 0.0555556 +s -0.122136 -0.871646 0.165419 0.0185185 +s -0.0854653 -0.82339 0.122833 0.0185185 +s -0.135649 -0.799077 0.171592 0.0185185 +s -0.194213 -0.884071 0.153697 0.0185185 +s -0.207727 -0.811502 0.15987 0.0185185 +s -0.229621 -0.84824 0.099389 0.0185185 +s -0.14403 -0.908384 0.104938 0.0185185 +s -0.179437 -0.872552 0.0506299 0.0185185 +s -0.107359 -0.860128 0.062352 0.0185185 +s 0.0133465 -0.69376 0.111111 0.0555556 +s 0.0838533 -0.674309 0.122833 0.0185185 +s 0.0539145 -0.666287 0.0555556 0.0185185 +s 0.0325183 -0.62221 0.111111 0.0185185 +s 0.0432853 -0.701782 0.178389 0.0185185 +s -0.00804971 -0.649684 0.166667 0.0185185 +s -0.0272215 -0.721234 0.166667 0.0185185 +s 0.0646815 -0.745859 0.122833 0.0185185 +s -0.0058253 -0.76531 0.111111 0.0185185 +s 0.0347427 -0.737837 0.0555556 0.0185185 +s -0.172546 -0.643951 0.222222 0.0555556 +s -0.142305 -0.61371 0.282703 0.0185185 +s -0.100996 -0.624779 0.222222 0.0185185 +s -0.153374 -0.5724 0.222222 0.0185185 +s -0.213855 -0.632882 0.282703 0.0185185 +s -0.224924 -0.591572 0.222222 0.0185185 +s -0.244096 -0.663122 0.222222 0.0185185 +s -0.161477 -0.68526 0.282703 0.0185185 +s -0.191718 -0.715501 0.222222 0.0185185 +s -0.120168 -0.696329 0.222222 0.0185185 +s -0.343435 -0.786005 7.27893e-17 0.0555556 +s -0.392454 -0.82165 0.0425863 0.0185185 +s -0.323622 -0.80984 0.0672777 0.0185185 +s -0.370774 -0.753118 0.0604812 0.0185185 +s -0.412267 -0.797815 -0.0246914 0.0185185 +s -0.390587 -0.729282 -0.00679642 0.0185185 +s -0.363249 -0.76217 -0.0672777 0.0185185 +s -0.365116 -0.854538 -0.0178949 0.0185185 +s -0.316097 -0.818893 -0.0604812 0.0185185 +s -0.296284 -0.842728 0.00679642 0.0185185 +s -0.358439 -0.594141 0.111111 0.0555556 +s -0.389119 -0.547576 0.15987 0.0185185 +s -0.317871 -0.566667 0.166667 0.0185185 +s -0.339267 -0.522591 0.111111 0.0185185 +s -0.429687 -0.57505 0.104315 0.0185185 +s -0.379835 -0.550064 0.0555556 0.0185185 +s -0.399007 -0.621614 0.0555556 0.0185185 +s -0.408291 -0.619127 0.15987 0.0185185 +s -0.37761 -0.665691 0.111111 0.0185185 +s -0.337042 -0.638217 0.166667 0.0185185 +s -0.358439 -0.594141 -0.111111 0.0555556 +s -0.409774 -0.542042 -0.122833 0.0185185 +s -0.379835 -0.550064 -0.0555556 0.0185185 +s -0.339267 -0.522591 -0.111111 0.0185185 +s -0.388377 -0.586119 -0.178389 0.0185185 +s -0.317871 -0.566667 -0.166667 0.0185185 +s -0.337042 -0.638217 -0.166667 0.0185185 +s -0.428945 -0.613592 -0.122833 0.0185185 +s -0.37761 -0.665691 -0.111111 0.0185185 +s -0.399007 -0.621614 -0.0555556 0.0185185 +s -0.157543 -0.835815 -0.111111 0.0555556 +s -0.119969 -0.899353 -0.117284 0.0185185 +s -0.0844101 -0.836885 -0.099389 0.0185185 +s -0.130205 -0.868703 -0.0506299 0.0185185 +s -0.193102 -0.898283 -0.129006 0.0185185 +s -0.203337 -0.867633 -0.062352 0.0185185 +s -0.230676 -0.834745 -0.122833 0.0185185 +s -0.147308 -0.866465 -0.177765 0.0185185 +s -0.184881 -0.802927 -0.171592 0.0185185 +s -0.111748 -0.803997 -0.15987 0.0185185 +s -0.172546 -0.643951 -0.222222 0.0555556 +s -0.202787 -0.61371 -0.282703 0.0185185 +s -0.244096 -0.624779 -0.222222 0.0185185 +s -0.191718 -0.5724 -0.222222 0.0185185 +s -0.131237 -0.632882 -0.282703 0.0185185 +s -0.120168 -0.591572 -0.222222 0.0185185 +s -0.100996 -0.663122 -0.222222 0.0185185 +s -0.183615 -0.68526 -0.282703 0.0185185 +s -0.153374 -0.715501 -0.222222 0.0185185 +s -0.224924 -0.696329 -0.222222 0.0185185 +s 0.0133465 -0.69376 -0.111111 0.0555556 +s 0.0631985 -0.668775 -0.15987 0.0185185 +s -0.00804971 -0.649684 -0.166667 0.0185185 +s 0.0325183 -0.62221 -0.111111 0.0185185 +s 0.0845947 -0.712851 -0.104315 0.0185185 +s 0.0539145 -0.666287 -0.0555556 0.0185185 +s 0.0347427 -0.737837 -0.0555556 0.0185185 +s 0.0440268 -0.740325 -0.15987 0.0185185 +s -0.0058253 -0.76531 -0.111111 0.0185185 +s -0.0272215 -0.721234 -0.166667 0.0185185 +s 0.471405 -0.471405 1.11022e-16 0.166667 +s 0.690426 -0.508983 1.83812e-16 0.0555556 +s 0.755941 -0.484794 -0.0246914 0.0185185 +s 0.695668 -0.478434 -0.0672777 0.0185185 +s 0.7029 -0.436283 -0.00679642 0.0185185 +s 0.7507 -0.515343 0.0425863 0.0185185 +s 0.697658 -0.466832 0.0604812 0.0185185 +s 0.685185 -0.539531 0.0672777 0.0185185 +s 0.743468 -0.557494 -0.0178949 0.0185185 +s 0.677953 -0.581682 0.00679642 0.0185185 +s 0.683194 -0.551134 -0.0604812 0.0185185 +s 0.607487 -0.335322 -0.111111 0.0555556 +s 0.629404 -0.313405 -0.178389 0.0185185 +s 0.610996 -0.384191 -0.166667 0.0185185 +s 0.558618 -0.331813 -0.166667 0.0185185 +s 0.625895 -0.264535 -0.122833 0.0185185 +s 0.555109 -0.282943 -0.111111 0.0185185 +s 0.603979 -0.286452 -0.0555556 0.0185185 +s 0.678274 -0.316914 -0.122833 0.0185185 +s 0.656357 -0.33883 -0.0555556 0.0185185 +s 0.659866 -0.3877 -0.111111 0.0185185 +s 0.607487 -0.335322 0.111111 0.0555556 +s 0.659645 -0.283164 0.104315 0.0185185 +s 0.656357 -0.33883 0.0555556 0.0185185 +s 0.603979 -0.286452 0.0555556 0.0185185 +s 0.610775 -0.279656 0.15987 0.0185185 +s 0.555109 -0.282943 0.111111 0.0185185 +s 0.558618 -0.331813 0.166667 0.0185185 +s 0.663153 -0.332034 0.15987 0.0185185 +s 0.610996 -0.384191 0.166667 0.0185185 +s 0.659866 -0.3877 0.111111 0.0185185 +s 0.554344 -0.645066 0.111111 0.0555556 +s 0.610229 -0.668521 0.153697 0.0185185 +s 0.622977 -0.619787 0.099389 0.0185185 +s 0.585648 -0.598918 0.15987 0.0185185 +s 0.541596 -0.6938 0.165419 0.0185185 +s 0.517014 -0.624197 0.171592 0.0185185 +s 0.48571 -0.670344 0.122833 0.0185185 +s 0.578925 -0.714669 0.104938 0.0185185 +s 0.52304 -0.691213 0.062352 0.0185185 +s 0.591673 -0.665934 0.0506299 0.0185185 +s 0.471405 -0.471405 0.222222 0.0555556 +s 0.501645 -0.441164 0.282703 0.0185185 +s 0.542955 -0.452233 0.222222 0.0185185 +s 0.490576 -0.399854 0.222222 0.0185185 +s 0.430095 -0.460336 0.282703 0.0185185 +s 0.419026 -0.419026 0.222222 0.0185185 +s 0.399854 -0.490576 0.222222 0.0185185 +s 0.482473 -0.512714 0.282703 0.0185185 +s 0.452233 -0.542955 0.222222 0.0185185 +s 0.523783 -0.523783 0.222222 0.0185185 +s 0.335322 -0.607487 0.111111 0.0555556 +s 0.313405 -0.629404 0.178389 0.0185185 +s 0.384191 -0.610996 0.166667 0.0185185 +s 0.331813 -0.558618 0.166667 0.0185185 +s 0.264535 -0.625895 0.122833 0.0185185 +s 0.282943 -0.555109 0.111111 0.0185185 +s 0.286452 -0.603979 0.0555556 0.0185185 +s 0.316914 -0.678274 0.122833 0.0185185 +s 0.33883 -0.656357 0.0555556 0.0185185 +s 0.3877 -0.659866 0.111111 0.0185185 +s 0.554344 -0.645066 -0.111111 0.0555556 +s 0.616373 -0.681385 -0.129006 0.0185185 +s 0.617144 -0.607573 -0.122833 0.0185185 +s 0.609912 -0.649723 -0.062352 0.0185185 +s 0.553573 -0.718878 -0.117284 0.0185185 +s 0.547112 -0.687216 -0.0506299 0.0185185 +s 0.491544 -0.682558 -0.099389 0.0185185 +s 0.560805 -0.676727 -0.177765 0.0185185 +s 0.498776 -0.640408 -0.15987 0.0185185 +s 0.561576 -0.602915 -0.171592 0.0185185 +s 0.335322 -0.607487 -0.111111 0.0555556 +s 0.283164 -0.659645 -0.104315 0.0185185 +s 0.33883 -0.656357 -0.0555556 0.0185185 +s 0.286452 -0.603979 -0.0555556 0.0185185 +s 0.279656 -0.610775 -0.15987 0.0185185 +s 0.282943 -0.555109 -0.111111 0.0185185 +s 0.331813 -0.558618 -0.166667 0.0185185 +s 0.332034 -0.663153 -0.15987 0.0185185 +s 0.384191 -0.610996 -0.166667 0.0185185 +s 0.3877 -0.659866 -0.111111 0.0185185 +s 0.471405 -0.471405 -0.222222 0.0555556 +s 0.441164 -0.441164 -0.282703 0.0185185 +s 0.399854 -0.452233 -0.222222 0.0185185 +s 0.452233 -0.399854 -0.222222 0.0185185 +s 0.512714 -0.460336 -0.282703 0.0185185 +s 0.523783 -0.419026 -0.222222 0.0185185 +s 0.542955 -0.490576 -0.222222 0.0185185 +s 0.460336 -0.512714 -0.282703 0.0185185 +s 0.490576 -0.542955 -0.222222 0.0185185 +s 0.419026 -0.523783 -0.222222 0.0185185 diff --git a/glome-hs.cabal b/glome-hs.cabal new file mode 100644 --- /dev/null +++ b/glome-hs.cabal @@ -0,0 +1,22 @@ +Name: glome-hs +Version: 0.4.1 +Synopsis: ray tracer +Description: Ray Tracer capable of rendering a variety of primitives, + with support for CSG (difference and intersection of solids), + BIH-based acceleration structure, and ability to load NFF + format files. +License: GPL +License-file: LICENSE +Author: Jim Snow +Maintainer: Jim Snow +Copyright: Copyright 2008 Jim Snow +Homepage: http://syn.cs.pdx.edu/~jsnow/glome +Stability: experimental +Category: graphics +Build-Depends: base,haskell98,time,parallel,GLUT,OpenGL,random,array +build-type: Simple +Executable: glome +ghc-options: -fglasgow-exts -funbox-strict-fields +extensions: BangPatterns +Main-is: Glome.hs + diff --git a/make b/make new file mode 100644 --- /dev/null +++ b/make @@ -0,0 +1,9 @@ +#!/bin/bash +#rm *.o +#ghc -O3 --make Glome.hs +#ghc Glome.hs --make -O2 -threaded -fasm -optc-march=athlon64 -XFlexibleInstances -XTypeSynonymInstances +#ghc Glome.hs --make -O2 -fasm -fglasgow-exts -funbox-strict-fields -fbang-patterns -fexcess-precision -optc-ffast-math -optc-O2 -optc-mfpmath=sse -optc-msse2 +#ghc Glome.hs --make -fasm -fglasgow-exts -funbox-strict-fields -fbang-patterns -fexcess-precision -prof -auto + +runhaskell Setup.lhs configure --prefix=$HOME --user +runhaskell Setup.lhs build diff --git a/run b/run new file mode 100644 --- /dev/null +++ b/run @@ -0,0 +1,4 @@ +#!/bin/bash +#./Glome +RTS -N2 -sstderr -RTS +#./Glome +RTS +./dist/build/glome/glome