Multithreaded version of @shelajev Java port of the postcard raytracer by Fabien Sanglard

PathtracerMT.java

package org.example.pathtracer.multithreaded;

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 java.util.concurrent.CountDownLatch;

import static org.example.pathtracer.multithreaded.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 PathtracerMT
{
	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;
	}

	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 = 560, 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);
		}

		byte[] imageData = new byte[w * h * 3];

		int processorCount = Runtime.getRuntime().availableProcessors();

		int linesPerProcessor = h / processorCount;

		final CountDownLatch latch = new CountDownLatch(processorCount);

		for (int i = 0; i < processorCount; i++)
		{
			final int startingLine = h - 1 - (i * linesPerProcessor);
			final int pixelBufferOffset = i * linesPerProcessor;
			
			final int threadID = i;

			Thread worker = new Thread(new Runnable()
			{
				@Override public void run()
				{
					int pixel = w * pixelBufferOffset * 3;

					// For each line
					for (int y = startingLine; y > startingLine - linesPerProcessor; 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);

							imageData[pixel++] = (byte) (int) color.x;
							imageData[pixel++] = (byte) (int) color.y;
							imageData[pixel++] = (byte) (int) color.z;
						}
					}

					latch.countDown();
				}
			});

			worker.start();
		}

		latch.await();

		try (FileOutputStream fw = new FileOutputStream(fileName.toFile()))
		{
			fw.write(String.format("P6 %d %d 255 ", w, h).getBytes(StandardCharsets.US_ASCII));
			fw.write(imageData);
		}

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