Java port of the postcard raytracer: http://fabiensanglard.net/postcard_pathtracer/index.html

Pathtracer.java

package org.example.pathtracer;

import java.io.*;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.nio.file.Paths;
import java.util.Random;

import static org.example.pathtracer.Vec.*;

class Vec {
    public float x, y, z;

    public Vec(float a, float b, float c) {
        x = a;
        y = b;
        z = c;
    }

    public static Vec add(Vec q, Vec r) {
        return new Vec(q.x + r.x, q.y + r.y, q.z + r.z);
    }

    public static Vec sub(Vec q, Vec r) {
        return new Vec(q.x - r.x, q.y - r.y, q.z - r.z);
    }

    public static Vec mul(Vec q, Vec r) {
        return new Vec(q.x * r.x, q.y * r.y, q.z * r.z);
    }

    public static Vec mul(Vec q, float r) {
        return mul(q, new Vec(r, r, r));
    }

    public static float dot(Vec q, Vec r) {
        return q.x * r.x + q.y * r.y + q.z * r.z;
    }

    public static Vec invSqrt(Vec q) {
        return mul(q, (1.0f / (float) Math.sqrt(dot(q, q))));
    }

    public Vec copy() {
        return new Vec(this.x, this.y, this.z);
    }

    @Override
    public String toString() {
        return x + ", " + y + ", " + z;
    }
}

public class Pathtracer {

    public static Vec Vec(float a) {
        return new Vec(a, a, a);
    }

    public static Vec Vec(float a, float b) {
        return new Vec(a, b, 0);
    }

    public static Vec Vec(float a, float b, float c) {
        return new Vec(a, b, c);
    }

    private static float min(float l, float r) {
        return Math.min(l, r);
    }

    private static Random random = new Random();

    private static float randomVal() {
        return random.nextFloat();
    }

    private static float fmodf(float x, float y) {
        return x % y; // according to https://stackoverflow.com/a/2690516
    }

    private static float fabsf(float x) {
        return Math.abs(x);
    }

    private static float sqrtf(float x) {
        return (float) Math.sqrt(x);
    }

    private static float powf(float x, float y) {
        return (float) Math.pow(x, y);
    }

    private static float cosf(float x) {
        return (float) Math.cos(x);
    }

    private static float sinf(float x) {
        return (float) Math.sin(x);
    }

    // Rectangle CSG equation. Returns minimum signed distance from
    // space carved by
    // lowerLeft vertex and opposite rectangle vertex upperRight.
    static float BoxTest(Vec position, Vec lowerLeft, Vec upperRight) {
        lowerLeft = sub(position, lowerLeft);
        upperRight = sub(upperRight, position);
        return -min(
                min(
                        min(lowerLeft.x, upperRight.x),
                        min(lowerLeft.y, upperRight.y)),
                min(lowerLeft.z, upperRight.z)
        );
    }

    private static final int HIT_NONE = 0;
    private static final int HIT_LETTER = 1;
    private static final int HIT_WALL = 2;
    private static final int HIT_SUN = 3;

    private static final char[] letters =              // 15 two points lines
            ("5O5_" + "5W9W" + "5_9_" +        // P (without curve)
                    "AOEO" + "COC_" + "A_E_" +        // I
                    "IOQ_" + "I_QO" +               // X
                    "UOY_" + "Y_]O" + "WW[W" +        // A
                    "aOa_" + "aWeW" + "a_e_" + "cWiO") // R (without curve)
                    .toCharArray();

    // Two curves (for P and R in PixaR) with hard-coded locations.
    private static final Vec[] curves = new Vec[]{Vec(-11, 6), Vec(11, 6)};

    // Sample the world using Signed Distance Fields.
    private static float QueryDatabase(Vec position, int[] hitType) {
        float distance = Float.MAX_VALUE;
        Vec f = position.copy(); // Flattened position (z=0)
        f.z = 0;


        for (int i = 0; i < letters.length; i += 4) {
            Vec begin = mul(Vec(letters[i] - 79, letters[i + 1] - 79), .5f);
            Vec e = sub(mul(Vec(letters[i + 2] - 79, letters[i + 3] - 79), .5f), begin);
            Vec o = sub(f,
                    (add(
                            begin,
                            mul(e, min(
                                    -min(
                                            dot(sub(begin, f), e) / dot(e, e),
                                            0),
                                    1)))
                    ));
            distance = min(distance, dot(o, o)); // compare squared distance.
        }
        distance = sqrtf(distance); // Get real distance, not square distance.
        for (int i = 1; i >= 0; i--) {
            Vec o = sub(f, curves[i]);
            // I *think* this equivalent to the C++ 'conditional expression', see https://stackoverflow.com/a/16676940
            float temp = 0.0f;
            if (o.x > 0) {
                temp = fabsf(sqrtf(dot(o, o)) - 2);
            } else {
                o.y += o.y > 0 ? -2 : 2;
                temp = sqrtf(dot(o, o));
            }
            distance = min(distance, temp);
        }

        distance = powf(powf(distance, 8) + powf(position.z, 8), 0.125f) - 0.5f;
        hitType[0] = HIT_LETTER;

        float roomDist;
        roomDist = min(// min(A,B) = Union with Constructive solid geometry
                //-min carves an empty space
                -min(// Lower room
                        BoxTest(position, Vec(-30, -0.5f, -30), Vec(30, 18, 30)),
                        // Upper room
                        BoxTest(position, Vec(-25, 17, -25), Vec(25, 20, 25))
                ),
                BoxTest( // Ceiling "planks" spaced 8 units apart.
                        Vec(fmodf(fabsf(position.x), 8),
                                position.y,
                                position.z),
                        Vec(1.5f, 18.5f, -25),
                        Vec(6.5f, 20, 25)
                )
        );
        if (roomDist < distance) {
            distance = roomDist;
            hitType[0] = HIT_WALL;
        }

        float sun = 19.9f - position.y; // Everything above 19.9 is light source.
        if (sun < distance) {
            distance = sun;
            hitType[0] = HIT_SUN;
        }

        return distance;
    }

    // Perform signed sphere marching
    // Returns hitType 0, 1, 2, or 3 and update hit position/normal
    static int RayMarching(Vec origin, Vec direction, Vec[] hitPos, Vec[] hitNorm) {
        int[] hitType = {HIT_NONE};
        int noHitCount = 0;

        int[] no_use = {0};

        float d = 0f; // distance from closest object in world.


        // Signed distance marching
        for (float total_d = 0; total_d < 100; total_d += d) {
            hitPos[0] = add(origin, mul(direction, total_d));
            d = QueryDatabase(hitPos[0], hitType);

            if (d < .01 || ++noHitCount > 99) { // if we hit or don't hit for a while
                // update hitNorm
                Vec vec = Vec(QueryDatabase(add(hitPos[0], Vec(.01f, 0, 0)), no_use) - d,
                        QueryDatabase(add(hitPos[0], Vec(0, .01f, 0)), no_use) - d,
                        QueryDatabase(add(hitPos[0], Vec(0, 0, .01f)), no_use) - d);
                hitNorm[0] = invSqrt(vec);

                // return hitType, this should be
                return hitType[0];
            }
        }
        return HIT_NONE; // throw new RuntimeException("Should not be here: " + origin + ", " + direction + ", " + hitPos[0] + ", " + hitNorm[0] + ", d = " + d + ", noHitCount = " + noHitCount);
    }

    static Vec Trace(Vec origin, Vec direction) {
        Vec[] sampledPosition = {Vec(1)};
        Vec[] normal = {Vec(0)};
        Vec color = Vec(0);
        Vec attenuation = Vec(1);
        Vec lightDirection = invSqrt(Vec(.6f, .6f, 1f)); // Directional light

        for (int bounceCount = 2; bounceCount >= 0; bounceCount--) {
            int hitType = RayMarching(origin, direction, sampledPosition, normal);

            if (hitType == HIT_NONE) {
                break; // No hit. This is over, return color.
            }
            Vec norm = normal[0];
            if (hitType == HIT_LETTER) { // Specular bounce on a letter. No color acc.
                direction = add(direction, mul(norm, (dot(norm, direction) * -2)));
                origin = add(sampledPosition[0], mul(direction, 0.1f));
                attenuation = mul(attenuation, 0.2f); // Attenuation via distance traveled.
            }
            if (hitType == HIT_WALL) { // Wall hit uses color yellow?
                float incidence = dot(norm, lightDirection);
                float p = 6.283185f * randomVal();
                float c = randomVal();
                float s = sqrtf(1 - c);
                float g = norm.z < 0 ? -1 : 1;
                float u = -1 / (g + norm.z);
                float v = norm.x * norm.y * u;

                direction = add(add(mul(Vec(v, g + norm.y * norm.y * u, -norm.y), (cosf(p) * s)), mul(Vec(1 + g * norm.x * norm.x * u, g * v, -g * norm.x), (sinf(p) * s))), mul(norm, sqrtf(c)));

                origin = add(sampledPosition[0], mul(direction, .1f));

                attenuation = mul(attenuation, 0.2f);

                if (incidence > 0 &&
                        RayMarching(add(sampledPosition[0], mul(norm, .1f)),
                                lightDirection,
                                sampledPosition,
                                normal) == HIT_SUN) {
                    color = add(color, mul(mul(attenuation, Vec(500, 400, 100)), incidence));
                }
            }
            if (hitType == HIT_SUN) { //
                color = add(color, mul(attenuation, Vec(50, 80, 100)));
                break; // Sun Color
            }
        }

        return color;
    }

    public static void main(String[] args) throws Exception {
        long start = -System.currentTimeMillis();

        int w = 960, h = 540, samplesCount = 16;  //8;
        Vec position = Vec(-22f, 5f, 25f);
        Vec goal = invSqrt(sub(Vec(-3f, 4f, 0f), position));


        Vec left = mul(invSqrt(Vec(goal.z, 0, -goal.x)), (1.0f / w));

        // Cross-product to get the up vector
        Vec up = Vec(
                goal.y * left.z - goal.z * left.y,
                goal.z * left.x - goal.x * left.z,
                goal.x * left.y - goal.y * left.x);

        Path fileName = Paths.get(String.format("output-java-%d.ppm", samplesCount));
        System.out.println("File: " + fileName);
        if (Files.exists(fileName)) {
            Files.delete(fileName);
        }

        try (FileOutputStream fw = new FileOutputStream(fileName.toFile())) {
            fw.write(String.format("P6 %d %d 255 ", w, h).getBytes(StandardCharsets.US_ASCII));
            for (int y = h; y > 0; y--) {
                for (int x = w; x > 0; x--) {
                    Vec color = Vec(0);
                    for (int p = samplesCount; p > 0; p--) {
                        color = add(color, Trace(position, invSqrt(
                                add(add(goal, mul(left,
                                        ((x - (w / 2)) + randomVal()))), mul(up, ((y - (h / 2)) + randomVal()))))));
                    }

                    // Reinhard tone mapping
                    color = add(mul(color, (1.0f / samplesCount)), Vec(14.0f / 241));
                    Vec o = add(color, Vec(1));
                    color = mul(Vec(color.x / o.x, color.y / o.y, color.z / o.z), 255);

                    fw.write(new byte[]{(byte) (int) color.x, (byte) (int) color.y, (byte) (int) color.z});
                }
            }
        }

        start += System.currentTimeMillis();
        System.out.println(start / 1000 + " s");
    }
}