Implemented anti-aliasing with permuted halton sequence.

src/algorithm.rs

@@ -0,0 +1,47 @@
+#![allow(dead_code)]
+
+use std;
+
+/// Partitions a slice in-place with the given unary predicate, returning
+/// the index of the first element for which the predicate evaluates
+/// false.
+///
+/// The predicate is executed precisely once on every element in
+/// the slice, and is allowed to modify the elements.
+pub fn partition<T, F>(slc: &mut [T], pred: F) -> usize
+    where F: Fn(&mut T) -> bool
+{
+    // This version uses raw pointers and pointer arithmetic to squeeze more
+    // performance out of the code.
+    unsafe {
+        let mut a = slc.as_mut_ptr();
+        let mut b = a.offset(slc.len() as isize);
+        let start = a as usize;
+
+        loop {
+            loop {
+                if a == b {
+                    return ((a as usize) - start) / std::mem::size_of::<T>();
+                }
+                if !pred(&mut *a) {
+                    break;
+                }
+                a = a.offset(1);
+            }
+
+            loop {
+                b = b.offset(-1);
+                if a == b {
+                    return ((a as usize) - start) / std::mem::size_of::<T>();
+                }
+                if pred(&mut *b) {
+                    break;
+                }
+            }
+
+            std::ptr::swap(a, b);
+
+            a = a.offset(1);
+        }
+    }
+}

src/halton.py

@@ -0,0 +1,190 @@
+#!/usr/bin/env python
+
+# Copyright (c) 2012 Leonhard Gruenschloss (leonhard@gruenschloss.org)
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights to
+# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+# of the Software, and to permit persons to whom the Software is furnished to do
+# so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+# Adapted to generate Rust instead of C by Nathan Vegdahl
+# Generate Rust code for evaluating Halton points with Faure-permutations for different bases.
+
+# How many components to generate.
+num_dimensions = 256
+
+# Check primality. Not optimized, since it's not performance-critical.
+def is_prime(p):
+    for i in range(2, p):
+        if not p % i:
+            return False
+    return True
+
+# Init prime number array.
+primes = []
+candidate = 1
+for i in range(num_dimensions):
+    while (True):
+        candidate += 1
+        if (is_prime(candidate)):
+            break;
+    primes.append(candidate)
+
+# Compute the Faure digit permutation for 0, ..., b - 1.
+def get_faure_permutation(b):
+    if b < 2:
+        return (0,)
+
+    elif b == 2:
+        return (0, 1)
+
+    elif b & 1: # odd
+        c = (b - 1) / 2
+
+        def faure_odd(i):
+            if i == c:
+                return c
+
+            f = faure[b - 1][i - int(i > c)]
+            return f + int(f >= c)
+
+        return tuple((faure_odd(i) for i in range(b)))
+
+    else: # even
+        c = b / 2
+
+        def faure_even(i):
+            if i < c:
+                return 2 * faure[c][i]
+            else:
+                return 2 * faure[c][i - c] + 1
+
+        return tuple((faure_even(i) for i in range(b)))
+
+# Init Faure permutations.
+faure = []
+for b in range(primes[-1] + 1):
+    faure.append(get_faure_permutation(b))
+
+# Compute the radical inverse with Faure permutations.
+def invert(base, index, digits):
+    result = 0
+    for i in range(digits):
+        index, remainder = divmod(index, base)
+        result = result * base + faure[base][remainder]
+    return result
+
+# Print the beginning bits of the file
+print '''#![allow(dead_code)]
+// Copyright (c) 2012 Leonhard Gruenschloss (leonhard@gruenschloss.org)
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights to
+// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+// of the Software, and to permit persons to whom the Software is furnished to do
+// so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+// This file is automatically generated.
+
+// Compute points of the Halton sequence with with Faure-permutations for different bases.
+
+pub const MAX_DIMENSION: u32 = %d;
+''' % num_dimensions
+
+# Print the sampling function
+print '''
+pub fn sample(dimension: u32, index: u32) -> f32 {
+    match dimension {
+'''
+
+for i in range(num_dimensions):
+    print '        %d => halton%d(index),' % (i, primes[i])
+
+print '''
+        _ => panic!("Exceeded max dimensions."),
+    }
+}
+'''
+
+# Print the special-cased first dimension
+print '''
+// Special case: radical inverse in base 2, with direct bit reversal.
+#[inline(always)]
+fn halton2(mut index: u32) -> f32 {
+    index = (index << 16) | (index >> 16);
+    index = ((index & 0x00ff00ff) << 8) | ((index & 0xff00ff00) >> 8);
+    index = ((index & 0x0f0f0f0f) << 4) | ((index & 0xf0f0f0f0) >> 4);
+    index = ((index & 0x33333333) << 2) | ((index & 0xcccccccc) >> 2);
+    index = ((index & 0x55555555) << 1) | ((index & 0xaaaaaaaa) >> 1);
+    return (index as f32) * (1.0 / ((1u64 << 32) as f32));
+}
+'''
+
+for i in range(1, num_dimensions): # Skip base 2.
+    base = primes[i]
+
+    # Based on the permutation table size, we process multiple digits at once.
+    digits = 1
+    pow_base = base
+    while pow_base * base <= 500: # Maximum permutation table size.
+        pow_base *= base
+        digits += 1
+
+    max_power = pow_base
+    powers = []
+    while max_power * pow_base < (1 << 32): # 32-bit unsigned precision
+        powers.append(max_power)
+        max_power *= pow_base
+
+    # Build the permutation table.
+    perm = []
+    for j in range(pow_base):
+        perm.append(invert(base, j, digits))
+        
+    power = max_power / pow_base
+    print '''
+
+#[inline(always)]
+fn halton%d(index: u32) -> f32 {
+    const PERM%d: [u16; %d] = [%s];
+    ''' % (base, base, len(perm), ', '.join(str(k) for k in perm))
+
+    print '''    return (unsafe{*PERM%d.get_unchecked((index %% %d) as usize)} as u32 * %d +''' % \
+        (base, pow_base, power)
+
+    # Advance to next set of digits.
+    div = 1
+    while power / pow_base > 1:
+        div *= pow_base
+        power /= pow_base
+        print '            unsafe{*PERM%d.get_unchecked(((index / %d) %% %d) as usize)} as u32 * %d +' % (base, div, pow_base, power)
+
+    print '''            unsafe{*PERM%d.get_unchecked(((index / %d) %% %d) as usize)} as u32) as f32 * (0.999999940395355224609375f32 / (%du32 as f32)); // Results in [0,1).
+}
+''' % (base, div * pow_base, pow_base, max_power)
+

src/main.rs

@@ -2,6 +2,7 @@ extern crate rustc_serialize;
 extern crate docopt;
 
 mod math;
+mod algorithm;
 mod lerp;
 mod float4;
 mod ray;
@@ -9,6 +10,7 @@ mod bbox;
 mod data_tree;
 mod image;
 mod triangle;
+mod halton;
 
 use std::path::Path;
 
@@ -16,7 +18,7 @@ use docopt::Docopt;
 
 use image::Image;
 use data_tree::DataTree;
-use math::{Point, Vector};
+use math::{Point, Vector, fast_logit};
 use ray::Ray;
 
 // ----------------------------------------------------------------
@@ -45,6 +47,18 @@ struct Args {
 
 // ----------------------------------------------------------------
 
+fn hash_u32(n: u32, seed: u32) -> u32 {
+    let mut hash = n;
+
+    for _ in 0..3 {
+        hash = hash.wrapping_mul(1936502639);
+        hash ^= hash.wrapping_shr(16);
+        hash = hash.wrapping_add(seed);
+    }
+
+    return hash;
+}
+
 fn main() {
     // Parse command line arguments.
     let args: Args = Docopt::new(USAGE.replace("<VERSION>", VERSION))
@@ -64,14 +78,31 @@ fn main() {
     let mut img = Image::new(512, 512);
     for y in 0..img.height() {
         for x in 0..img.width() {
-            let ray = Ray::new(Point::new(x as f32, y as f32, 0.0),
+            let mut r = 0.0;
-                               Vector::new(0.0, 0.0, 1.0));
+            let mut g = 0.0;
-            if let Some((_, u, v)) = triangle::intersect_ray(&ray, (p1, p2, p3)) {
+            let mut b = 0.0;
-                let r = (u * 255.0) as u8;
+            let offset = hash_u32(((x as u32) << 16) ^ (y as u32), 0);
-                let g = (v * 255.0) as u8;
+            const SAMPLES: usize = 16;
-                let b = ((1.0 - u - v) * 255.0).max(0.0) as u8;
+            for si in 0..SAMPLES {
-                img.set(x, y, (r, g, b));
+                let ray = Ray::new(Point::new(x as f32 +
+                                              fast_logit(halton::sample(0, offset + si as u32),
+                                                         1.5),
+                                              y as f32 +
+                                              fast_logit(halton::sample(3, offset + si as u32),
+                                                         1.5),
+                                              0.0),
+                                   Vector::new(0.0, 0.0, 1.0));
+                if let Some((_, u, v)) = triangle::intersect_ray(&ray, (p1, p2, p3)) {
+                    r += u;
+                    g += v;
+                    b += (1.0 - u - v).max(0.0);
+                }
             }
+            r *= 255.0 / SAMPLES as f32;
+            g *= 255.0 / SAMPLES as f32;
+            b *= 255.0 / SAMPLES as f32;
+
+            img.set(x, y, (r as u8, g as u8, b as u8));
         }
     }
     let _ = img.write_binary_ppm(Path::new(&args.arg_imgpath));

src/math/mod.rs

@@ -30,3 +30,24 @@ pub trait CrossProduct
 pub fn cross<T: CrossProduct>(a: T, b: T) -> T {
     a.cross(b)
 }
+
+// Adapted from from http://fastapprox.googlecode.com
+pub fn fast_ln(x: f32) -> f32 {
+    use std::mem::transmute_copy;
+
+    let mut y = unsafe { transmute_copy::<f32, u32>(&x) as f32 };
+    y *= 8.2629582881927490e-8;
+    return y - 87.989971088;
+}
+
+pub fn fast_logit(p: f32, width: f32) -> f32 {
+    let n = 0.001 + (p * 0.998);
+
+    fast_ln((n / (1.0 - n))) * width * (0.6266 / 4.0)
+}
+
+pub fn logit(p: f32, width: f32) -> f32 {
+    let n = 0.001 + (p * 0.998);
+
+    (n / (1.0 - n)).ln() * width * (0.6266 / 4.0)
+}

src/math/point.rs

@@ -25,6 +25,20 @@ impl Point {
     pub fn norm(&self) -> Point {
         Point { co: self.co / self.co[3] }
     }
+
+    pub fn min(&self, other: Point) -> Point {
+        let n1 = self.norm();
+        let n2 = other.norm();
+
+        Point { co: n1.co.v_min(n2.co) }
+    }
+
+    pub fn max(&self, other: Point) -> Point {
+        let n1 = self.norm();
+        let n2 = other.norm();
+
+        Point { co: n1.co.v_max(n2.co) }
+    }
 }
 
 

src/triangle.rs

@@ -29,6 +29,6 @@ pub fn intersect_ray(ray: &Ray, tri: (Point, Point, Point)) -> Option<(f32, f32,
         let t = dot(edge2, qvec) * inv_det;
         return Some((t, u, v));
     } else {
-        return None;  // Ray is parallell to the plane of the triangle
+        return None;
     }
 }

todo.txt

@@ -0,0 +1,2 @@
+- Implement a basic camera.
+- Implement a BVH for lists of triangles.