src/PlotterPaper.cc at main · doomsdayone/openmsx

◆ main openmsx / src / PlotterPaper.cc

16.1 KiB · C++ · last changed chore(): initial commit 1 month ago source blame history raw

#include "PlotterPaper.hh"

#include "ranges.hh"

#include <array>
#include <cmath>
#include <limits>
#include <map>
#include <span>

namespace openmsx {

// Table containing half circle-width for every y-position.
// Zoomed in factor 16x, so:
// * y-coordinates are 16x larger
// * reported width is also 16x larger
static std::span<int> get_dot_table(int radius16)
{
	static std::map<int, MemBuffer<int>> cache;
	auto [it, inserted] = cache.try_emplace(radius16);
	auto& table = it->second;
	if (inserted) { // was not yet present
		// Table covers full vertical extent: y from -(radius16+16) to +(radius16+16)
		table.resize(2 * (radius16 + 16) + 1);

		int offset = radius16 + 16;  // Center of table (y=0 relative to circle center)
		int r_sq = radius16 * radius16;

		// Compute half-widths for all y positions relative to circle center
		for (int y_rel = -(radius16 + 16); y_rel <= (radius16 + 16); ++y_rel) {
			int y_sq = y_rel * y_rel;

			if (y_sq > r_sq) {
				// Outside circle
				table[offset + y_rel] = 0;
			} else {
				// Inside circle: compute half-width
				// a = sqrt(r^2 - y^2)
				float a_float = std::sqrt(float(r_sq - y_sq));
				table[offset + y_rel] = int(a_float);
			}
		}
	}
	return table;
}

// Rasterize a stationary dot (pen dwells at a point)
// ink_color: color of ink (gl::vec3{r,g,b} for RGB, or float for grayscale)
//     For subtractive model: 0=no ink (white), 1=full ink (black)
//     To control ink amount, pre-multiply ink_color (e.g., ink_color * 0.5f)
// Uses 16×16 sub-pixel grid (256 shading levels)
template<typename InkType>
BoundingBox rasterize_dot(Canvas<InkType>& canvas, gl::vec2 center, float pen_radius, InkType ink_color)
{
	// Compute pixel bounding box
	BoundingBox bbox{
		.tl = max(gl::ivec2{0}, floor(center - gl::vec2(pen_radius))),
		.br = min(canvas.size(), ceil(center + gl::vec2(pen_radius)))
	};
	if (bbox.empty()) return bbox; // Early exit if no pixels to draw

	// Quantized center position (in 1/16 pixel units)
	int xPos16 = int(16.0f * center.x);
	int yPos16 = int(16.0f * center.y);
	int radius16 = int(16.0f * pen_radius);
	auto table = get_dot_table(radius16);

	// Starting y offset in lookup table
	int y = 16 * bbox.tl.y - yPos16 + 16 + radius16;

	for (int yy = bbox.tl.y; yy < bbox.br.y; ++yy, y += 16) {
		auto line = canvas.getLine(yy);
		auto table_slice = subspan<16>(table, y);

		// Calculate which pixel contains the Y-axis (x=0, circle center)
		int x_start = 16 * bbox.tl.x - xPos16;

		// Find the pixel that crosses the Y-axis
		// Y-axis at x=0, middle pixel has x in [-16, 0)
		int xx_mid = bbox.tl.x + (-x_start + 15) / 16;
		xx_mid = std::clamp(xx_mid, bbox.tl.x, bbox.br.x - 1);

		// Left half: pixels before middle (x+16 < 0)
		int x = x_start;
		for (int xx = bbox.tl.x; xx < xx_mid; ++xx) {
			int sum = 0;
			for (int i = 0; i < 16; ++i) {
				int a = table_slice[i];
				int coverage = std::clamp(16 + a + x, 0, 16);
				sum += coverage;
			}
			// Accumulate into float canvas (normalized to [0, 1])
			line[xx] += ink_color * (sum / 256.0f);
			x += 16;
		}

		// Middle pixel: exactly 1 (if Y-axis within bounding box)
		if (xx_mid < bbox.br.x) {
			int sum = 0;
			for (int i = 0; i < 16; ++i) {
				int a = table_slice[i];
				int left_coverage = std::clamp(16 + a + x, 0, 2 * a);
				int right_coverage = std::clamp(a - x, 0, 16);
				int coverage = (x < -a) ? left_coverage : right_coverage;
				sum += coverage;
			}
			line[xx_mid] += ink_color * (sum / 256.0f);
			x += 16;
		}

		// Right half: pixels after middle (x >= 0)
		for (int xx = xx_mid + 1; xx < bbox.br.x; ++xx) {
			int sum = 0;
			for (int i = 0; i < 16; ++i) {
				int a = table_slice[i];
				int coverage = std::clamp(a - x, 0, 16);
				sum += coverage;
			}
			line[xx] += ink_color * (sum / 256.0f);
			x += 16;
		}
	}
	return bbox;
}


// Analytical uniform disc: compute ink contribution from pen motion along a segment
// Does NOT include dwell time at endpoints - add those separately if needed
// Ink is automatically normalized: alpha=1.0 at center of long segment
// Caller should normalize endpoint dwells by dividing by (2*pen_radius) for consistency
// Uses 8x MSAA for anti-aliasing
template<typename InkType>
BoundingBox rasterize_segment_motion(Canvas<InkType>& canvas, gl::vec2 A, gl::vec2 B, float pen_radius, InkType ink_color)
{
	// MSAA 8x sample pattern (positions within pixel, in range [0, 1])
	// (0,0) = top-left corner, (1,1) = bottom-right corner, (0.5,0.5) = center
	// Rotated grid pattern for good quality anti-aliasing
	static constexpr std::array<gl::vec2, 8> msaa_8x = {{
		{0.125f, 0.375f}, {0.625f, 0.125f}, {0.375f, 0.875f}, {0.875f, 0.625f},
		{0.125f, 0.875f}, {0.875f, 0.125f}, {0.375f, 0.375f}, {0.625f, 0.625f}
	}};

	auto D = B - A;
	float L_sq = dot(D, D);
	if (L_sq < 1e-18f) return BoundingBox{}; // degenerate segment

	float L = std::sqrt(L_sq);
	gl::vec2 T = D * (1.0f / L); // unit tangent
	gl::vec2 N = {-T.y, T.x};    // unit normal

	// Pre-transform MSAA samples into (T, N) coordinate system (once per segment)
	// Exploits bilinearity: dot(a+b, c) = dot(a, c) + dot(b, c)
	std::array<gl::vec2, msaa_8x.size()> msaa_TN;
	for (size_t i = 0; i < msaa_8x.size(); ++i) {
		msaa_TN[i] = {dot(msaa_8x[i], T), dot(msaa_8x[i], N)};
	}

	// Vertical bounding box (using vector math, but only y component is used here)
	auto segment_min = min(A, B);
	auto segment_max = max(A, B);
	int y_min = std::max(0, int(std::floor(segment_min.y - pen_radius)));
	int y_max = std::min(canvas.size().y - 1, int(std::ceil(segment_max.y + pen_radius)));

	// Normalization: max ink is 2*pen_radius, scale by 1/(2*pen_radius*num_samples) for MSAA
	// Pre-multiply by 2 to optimize the common fast-path (avoids multiply in inner loop)
	constexpr float scale_factor = 1.0f / (2.0f * msaa_8x.size());
	const float pen_radius_sq = pen_radius * pen_radius;
	const InkType ink_contribution_double = ink_color * (2.0f * scale_factor / pen_radius);

	// Lambda to process a pixel at (x, dy) and write to line if non-zero
	// dy is the y-coordinate relative to A.y (computed once per scanline)
	auto process_pixel = [&](std::span<InkType> line, int x, float dy) {
		// Transform pixel position into (T, N) coordinate system (once per pixel)
		gl::vec2 pixel_minus_A = {float(x) - A.x, dy};
		gl::vec2 pixel_TN = {dot(pixel_minus_A, T), dot(pixel_minus_A, N)};

		// Accumulate half-ink (optimizes fast path: avoids *2 multiply for middle section)
		float half_ink_sum = 0.0f;
		for (const auto& sample_tn : msaa_TN) {
			// Exploit bilinearity: distances are just additions in (T,N) space
			float d_parallel = pixel_TN.x + sample_tn.x;  // Distance along line
			float d_perp = pixel_TN.y + sample_tn.y;      // Perpendicular distance (signed)
			float d_perp_sq = d_perp * d_perp;	      // Squared distance (positive)

			// Early rejection: too far perpendicular to line
			if (d_perp_sq > pen_radius_sq) continue;

			// Compute half-width of integration interval
			float h = std::sqrt(pen_radius_sq - d_perp_sq);

			// Fast path: middle section (most common case for long lines)
			// When h < d_parallel < L - h, the integration interval is 2*h
			// We accumulate h (not 2*h) and compensate with pre-doubled ink_contribution
			if (h < d_parallel && d_parallel < L - h) {
				half_ink_sum += h;
			} else {
				// Slow path: near start or end caps (needs clamping)
				float t_min = std::max(0.0f, d_parallel - h);
				float t_max = std::min(L, d_parallel + h);
				half_ink_sum += 0.5f * std::max(0.0f, t_max - t_min);
			}
		}
		if (half_ink_sum > 0.0f) {
			line[x] += ink_contribution_double * half_ink_sum;
		}
	};

	const float epsilon = 1e-6f;
	if (std::abs(T.y) > epsilon) {
		// Non-horizontal segment: use scanline intersection for tighter x bounds
		const float inv_Ty = 1.0f / T.y;
		const float dx_factor = (std::abs(N.x) > epsilon ? pen_radius / std::abs(N.x) : pen_radius);

		auto global_x_min = std::numeric_limits<int>::max();
		auto global_x_max = std::numeric_limits<int>::min();
		for (int y = y_min; y <= y_max; ++y) {
			auto line = canvas.getLine(y);
			float y_float = float(y);
			float dy = y_float - A.y;  // Relative to A (loop-invariant for x-loop)
			float Py = y_float + 0.5f;
			float t_center = (Py - A.y) * inv_Ty;
			float x_center = A.x + t_center * T.x;
			auto x_min = std::max(0, (int)std::floor(x_center - dx_factor));
			auto x_max = std::min(canvas.size().x - 1, (int)std::ceil(x_center + dx_factor));
			global_x_min = std::min(global_x_min, x_min);
			global_x_max = std::max(global_x_max, x_max);

			for (int x = x_min; x <= x_max; ++x) {
				process_pixel(line, x, dy);
			}
		}
		return BoundingBox{
			.tl = gl::ivec2{global_x_min,     y_min},
			.br = gl::ivec2{global_x_max + 1, y_max + 1}
		};
	} else {
		// Nearly horizontal segment: use axis-aligned bounding box
		int x_min = std::max(0, int(std::floor(segment_min.x - pen_radius)));
		int x_max = std::min(canvas.size().x - 1, int(std::ceil(segment_max.x + pen_radius)));

		for (int y = y_min; y <= y_max; ++y) {
			auto line = canvas.getLine(y);
			float dy = float(y) - A.y;  // Relative to A (loop-invariant for x-loop)
			for (int x = x_min; x <= x_max; ++x) {
				process_pixel(line, x, dy);
			}
		}
		return BoundingBox{
			.tl = gl::ivec2{x_min,     y_min},
			.br = gl::ivec2{x_max + 1, y_max + 1}
		};
	}
}

/*
// Rasterize a polyline with proper endpoint dwell times
template<typename InkType>
void rasterize_polyline(
	Canvas<InkType>& canvas, std::span<const gl::vec2> pts,
	float pen_radius, InkType ink_color, float dwell_time = 1.0f)
{
	if (pts.size() < 2) {
		if (pts.size() == 1) {
			// Single point: just a dot with specified dwell time
			rasterize_dot(canvas, pts[0], pen_radius, ink_color * dwell_time);
		}
		return;
	}

	// Normalize dwell to ink amount and pre-multiply with ink_color
	const float ink_amount = dwell_time / (2.0f * pen_radius);
	const InkType endpoint_ink = ink_color * ink_amount;

	// Dwell at start point (pen lands)
	rasterize_dot(canvas, pts[0], pen_radius, endpoint_ink);

	// All segments (just motion, no dwell at interior points)
	for (size_t i = 0; i + 1 < pts.size(); ++i) {
		rasterize_segment_motion(canvas, pts[i], pts[i + 1], pen_radius, ink_color);
	}

	// Dwell at end point (pen lifts)
	rasterize_dot(canvas, pts.back(), pen_radius, endpoint_ink);
}

// Convenience helper. Is this useful enough? Remove?
template <typename InkType>
void rasterize_segment(
	Canvas<InkType>& canvas, gl::vec2 A, gl::vec2 B,
	float pen_radius, InkType ink_color, float dwell_time = 1.0f)
{
	std::array pts = {A, B};
	rasterize_polyline(canvas, pts, pen_radius, ink_color, dwell_time);
}*/

////

// InkType: 'gl::vec3' for RGB, 'float' for grayscale
// Deduce pixel type from ink type
template<typename InkType> struct PixelTypeFor;
template<> struct PixelTypeFor<gl::vec3> { using type = RgbPixel; };
template<> struct PixelTypeFor<float> { using type = GrayPixel; };
template<typename InkType> using PixelTypeFor_t = typename PixelTypeFor<InkType>::type;

template<typename Sampler> using InkTypeForSampler = decltype(std::declval<Sampler>()(0, 0));

template<typename Sampler>
Canvas<PixelTypeFor_t<InkTypeForSampler<Sampler>>> ink_to_pixels(gl::ivec2 size, Sampler sampler)
{
	using InkType = InkTypeForSampler<Sampler>;
	using PixelType = PixelTypeFor_t<InkType>;

	// Saturation strength - models paper's limited ink absorption
	// Rational function: ink*(1+k) / (1+k*ink), scaled to 0-255 range
	// - Passes through (1,1): single stroke looks correct
	// - For overlaps (ink>1): output < ink (overlaps saturate smoothly)
	constexpr float k = 4.0f;
	constexpr float numerator = (1.0f + k) * 255.0f;

	Canvas<PixelType> output(size, PixelType(255));
	for (int y = 0; y < size.y; ++y) {
		auto out_line = output.getLine(y);
		for (int x = 0; x < size.x; ++x) {
			InkType accumulated_ink = sampler(x, y);
			if (accumulated_ink == InkType(0.0f)) {
				continue;  // Already initialized to white (255)
			}

			// Apply saturation function in 0-255 range (formula works for both float and gl::vec3)
			InkType saturated_ink_255 = accumulated_ink * numerator / (InkType(1.0f) + k * accumulated_ink);

			// Apply subtractive ink model: white paper (255) minus saturated ink
			InkType pixel_value = InkType(255.0f) - saturated_ink_255;
			using gl::clamp;  // Bring both clamp functions into scope ..
			using std::clamp; // .. overload resolution picks the right one
			out_line[x] = PixelType(clamp(pixel_value, InkType(0.0f), InkType(255.0f)));
		}
	}
	return output;
}

////

PlotterPaper::PlotterPaper(gl::ivec2 size)
	: paper(size)
	, damage{.tl = gl::ivec2{0}, .br = size}
{
}

void PlotterPaper::draw_dot(gl::vec2 pos, float radius, gl::vec3 color)
{
	anythingPlotted = true;
	auto box = rasterize_dot(paper, pos, radius, color);
	damage = BoundingBox::merge(damage, box);
}

void PlotterPaper::draw_motion(gl::vec2 from, gl::vec2 to, float radius, gl::vec3 color)
{
	anythingPlotted = true;
	auto box = rasterize_segment_motion(paper, from, to, radius, color);
	damage = BoundingBox::merge(damage, box);
}

Canvas<RgbPixel> PlotterPaper::getRGB() const
{
	auto out_size = size();
	auto sampler = [&](int x, int y) { return paper.getLine(y)[x]; };
	return ink_to_pixels(out_size, sampler);
}

gl::Texture& PlotterPaper::updateTexture(gl::ivec2 targetSize)
{
	assert(targetSize.x >= 1);
	assert(targetSize.y >= 1);

	auto fullSize = paper.size();
	auto clampedTargetSize = min(targetSize, gl::ivec2(1024)); // ensure below 2048 x 2028
	auto factor = max_component(trunc(fullSize / clampedTargetSize));
	auto newTexSize = max(fullSize / factor, gl::ivec2(1));
	bool sizeChanged = newTexSize != texSize;

	if (texture.get() && !sizeChanged && damage.empty()) {
		return texture; // already up to date
	}
	if (sizeChanged) {
		// full update if size changed
		damage = BoundingBox{.tl = gl::ivec2{0}, .br = paper.size()};
	}
	auto dstTl = damage.tl / factor;
	auto dstBr = damage.br / factor;
	auto dstSize = max(dstBr - dstTl, gl::ivec2(1));
	auto srcOffset = dstTl * factor;

	auto sample_direct = [&](int x, int y) { return paper.getLine(y + srcOffset.y)[x + srcOffset.x]; };
	auto down_sample_generic = [&](int x, int y) { // run-time 'factor'
		auto srcPos = gl::ivec2(x, y) * factor + srcOffset;
		auto sum = gl::vec3(0.0f);
		auto factorXY = min(paper.size(), gl::ivec2(factor));
		for (int sy = 0; sy < factorXY.y; ++sy) {
			auto line = subspan(paper.getLine(srcPos.y + sy), srcPos.x, factorXY.x);
			for (int sx = 0; sx < factorXY.x; ++sx) {
				sum += line[sx];
			}
		}
		return sum * (1.0f / float(factor * factor));
	};
	auto down_sample_N = [&]<int FACTOR>(int x, int y) { // compile-time 'FACTOR'
		auto srcPos = gl::ivec2(x, y) * factor + srcOffset;
		auto sum = gl::vec3(0.0f);
		for (int sy = 0; sy < FACTOR; ++sy) {
			auto line = subspan<FACTOR>(paper.getLine(srcPos.y + sy), srcPos.x);
			for (int sx = 0; sx < FACTOR; ++sx) {
				sum += line[sx];
			}
		}
		return sum * (1.0f / float(FACTOR * FACTOR));
	};


	auto rgb = [&]{
		// specializations for small factors, this gives a reasonably large speedup
		switch (factor) {
		case 1:  return ink_to_pixels(dstSize, sample_direct);
		case 2:  return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<2>(x, y); });
		case 3:  return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<3>(x, y); });
		case 4:  return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<4>(x, y); });
		case 5:  return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<5>(x, y); });
		case 6:  return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<6>(x, y); });
		case 7:  return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<7>(x, y); });
		case 8:  return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<8>(x, y); });
		default: return ink_to_pixels(dstSize, down_sample_generic);
		}
	}();

	if (!texture.get()) {
		texture.allocate();
		texture.setInterpolation(true);
	} else {
		texture.bind();
	}
	if (sizeChanged) {
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, newTexSize.x, newTexSize.y, 0, GL_RGB,
		             GL_UNSIGNED_BYTE, &rgb.getLine(0)[0].x);
	} else {
		glTexSubImage2D(GL_TEXTURE_2D, 0, dstTl.x, dstTl.y, dstSize.x, dstSize.y, GL_RGB,
		                 GL_UNSIGNED_BYTE, &rgb.getLine(0)[0].x);
	}

	damage.clear();
	texSize = newTexSize;
	return texture;
}

} // namespace openmsx

`/`	Search
`g` `c`	Code
`g` `h`	History
`g` `b`	Bookmarks
`g` `t`	Trust
`y`	Permalink to this state (files and trees)
`?`	This sheet

1	#include "PlotterPaper.hh"
2
3	#include "ranges.hh"
4
5	#include <array>
6	#include <cmath>
7	#include <limits>
8	#include <map>
9	#include <span>
10
11	namespace openmsx {
12
13	// Table containing half circle-width for every y-position.
14	// Zoomed in factor 16x, so:
15	// * y-coordinates are 16x larger
16	// * reported width is also 16x larger
17	static std::span<int> get_dot_table(int radius16)
18	{
19	static std::map<int, MemBuffer<int>> cache;
20	auto [it, inserted] = cache.try_emplace(radius16);
21	auto& table = it->second;
22	if (inserted) { // was not yet present
23	// Table covers full vertical extent: y from -(radius16+16) to +(radius16+16)
24	table.resize(2 * (radius16 + 16) + 1);
25
26	int offset = radius16 + 16; // Center of table (y=0 relative to circle center)
27	int r_sq = radius16 * radius16;
28
29	// Compute half-widths for all y positions relative to circle center
30	for (int y_rel = -(radius16 + 16); y_rel <= (radius16 + 16); ++y_rel) {
31	int y_sq = y_rel * y_rel;
32
33	if (y_sq > r_sq) {
34	// Outside circle
35	table[offset + y_rel] = 0;
36	} else {
37	// Inside circle: compute half-width
38	// a = sqrt(r^2 - y^2)
39	float a_float = std::sqrt(float(r_sq - y_sq));
40	table[offset + y_rel] = int(a_float);
41	}
42	}
43	}
44	return table;
45	}
46
47	// Rasterize a stationary dot (pen dwells at a point)
48	// ink_color: color of ink (gl::vec3{r,g,b} for RGB, or float for grayscale)
49	// For subtractive model: 0=no ink (white), 1=full ink (black)
50	// To control ink amount, pre-multiply ink_color (e.g., ink_color * 0.5f)
51	// Uses 16×16 sub-pixel grid (256 shading levels)
52	template<typename InkType>
53	BoundingBox rasterize_dot(Canvas<InkType>& canvas, gl::vec2 center, float pen_radius, InkType ink_color)
54	{
55	// Compute pixel bounding box
56	BoundingBox bbox{
57	.tl = max(gl::ivec2{0}, floor(center - gl::vec2(pen_radius))),
58	.br = min(canvas.size(), ceil(center + gl::vec2(pen_radius)))
59	};
60	if (bbox.empty()) return bbox; // Early exit if no pixels to draw
61
62	// Quantized center position (in 1/16 pixel units)
63	int xPos16 = int(16.0f * center.x);
64	int yPos16 = int(16.0f * center.y);
65	int radius16 = int(16.0f * pen_radius);
66	auto table = get_dot_table(radius16);
67
68	// Starting y offset in lookup table
69	int y = 16 * bbox.tl.y - yPos16 + 16 + radius16;
70
71	for (int yy = bbox.tl.y; yy < bbox.br.y; ++yy, y += 16) {
72	auto line = canvas.getLine(yy);
73	auto table_slice = subspan<16>(table, y);
74
75	// Calculate which pixel contains the Y-axis (x=0, circle center)
76	int x_start = 16 * bbox.tl.x - xPos16;
77
78	// Find the pixel that crosses the Y-axis
79	// Y-axis at x=0, middle pixel has x in [-16, 0)
80	int xx_mid = bbox.tl.x + (-x_start + 15) / 16;
81	xx_mid = std::clamp(xx_mid, bbox.tl.x, bbox.br.x - 1);
82
83	// Left half: pixels before middle (x+16 < 0)
84	int x = x_start;
85	for (int xx = bbox.tl.x; xx < xx_mid; ++xx) {
86	int sum = 0;
87	for (int i = 0; i < 16; ++i) {
88	int a = table_slice[i];
89	int coverage = std::clamp(16 + a + x, 0, 16);
90	sum += coverage;
91	}
92	// Accumulate into float canvas (normalized to [0, 1])
93	line[xx] += ink_color * (sum / 256.0f);
94	x += 16;
95	}
96
97	// Middle pixel: exactly 1 (if Y-axis within bounding box)
98	if (xx_mid < bbox.br.x) {
99	int sum = 0;
100	for (int i = 0; i < 16; ++i) {
101	int a = table_slice[i];
102	int left_coverage = std::clamp(16 + a + x, 0, 2 * a);
103	int right_coverage = std::clamp(a - x, 0, 16);
104	int coverage = (x < -a) ? left_coverage : right_coverage;
105	sum += coverage;
106	}
107	line[xx_mid] += ink_color * (sum / 256.0f);
108	x += 16;
109	}
110
111	// Right half: pixels after middle (x >= 0)
112	for (int xx = xx_mid + 1; xx < bbox.br.x; ++xx) {
113	int sum = 0;
114	for (int i = 0; i < 16; ++i) {
115	int a = table_slice[i];
116	int coverage = std::clamp(a - x, 0, 16);
117	sum += coverage;
118	}
119	line[xx] += ink_color * (sum / 256.0f);
120	x += 16;
121	}
122	}
123	return bbox;
124	}
125
126
127	// Analytical uniform disc: compute ink contribution from pen motion along a segment
128	// Does NOT include dwell time at endpoints - add those separately if needed
129	// Ink is automatically normalized: alpha=1.0 at center of long segment
130	// Caller should normalize endpoint dwells by dividing by (2*pen_radius) for consistency
131	// Uses 8x MSAA for anti-aliasing
132	template<typename InkType>
133	BoundingBox rasterize_segment_motion(Canvas<InkType>& canvas, gl::vec2 A, gl::vec2 B, float pen_radius, InkType ink_color)
134	{
135	// MSAA 8x sample pattern (positions within pixel, in range [0, 1])
136	// (0,0) = top-left corner, (1,1) = bottom-right corner, (0.5,0.5) = center
137	// Rotated grid pattern for good quality anti-aliasing
138	static constexpr std::array<gl::vec2, 8> msaa_8x = {{
139	{0.125f, 0.375f}, {0.625f, 0.125f}, {0.375f, 0.875f}, {0.875f, 0.625f},
140	{0.125f, 0.875f}, {0.875f, 0.125f}, {0.375f, 0.375f}, {0.625f, 0.625f}
141	}};
142
143	auto D = B - A;
144	float L_sq = dot(D, D);
145	if (L_sq < 1e-18f) return BoundingBox{}; // degenerate segment
146
147	float L = std::sqrt(L_sq);
148	gl::vec2 T = D * (1.0f / L); // unit tangent
149	gl::vec2 N = {-T.y, T.x}; // unit normal
150
151	// Pre-transform MSAA samples into (T, N) coordinate system (once per segment)
152	// Exploits bilinearity: dot(a+b, c) = dot(a, c) + dot(b, c)
153	std::array<gl::vec2, msaa_8x.size()> msaa_TN;
154	for (size_t i = 0; i < msaa_8x.size(); ++i) {
155	msaa_TN[i] = {dot(msaa_8x[i], T), dot(msaa_8x[i], N)};
156	}
157
158	// Vertical bounding box (using vector math, but only y component is used here)
159	auto segment_min = min(A, B);
160	auto segment_max = max(A, B);
161	int y_min = std::max(0, int(std::floor(segment_min.y - pen_radius)));
162	int y_max = std::min(canvas.size().y - 1, int(std::ceil(segment_max.y + pen_radius)));
163
164	// Normalization: max ink is 2pen_radius, scale by 1/(2pen_radius*num_samples) for MSAA
165	// Pre-multiply by 2 to optimize the common fast-path (avoids multiply in inner loop)
166	constexpr float scale_factor = 1.0f / (2.0f * msaa_8x.size());
167	const float pen_radius_sq = pen_radius * pen_radius;
168	const InkType ink_contribution_double = ink_color * (2.0f * scale_factor / pen_radius);
169
170	// Lambda to process a pixel at (x, dy) and write to line if non-zero
171	// dy is the y-coordinate relative to A.y (computed once per scanline)
172	auto process_pixel = [&](std::span<InkType> line, int x, float dy) {
173	// Transform pixel position into (T, N) coordinate system (once per pixel)
174	gl::vec2 pixel_minus_A = {float(x) - A.x, dy};
175	gl::vec2 pixel_TN = {dot(pixel_minus_A, T), dot(pixel_minus_A, N)};
176
177	// Accumulate half-ink (optimizes fast path: avoids *2 multiply for middle section)
178	float half_ink_sum = 0.0f;
179	for (const auto& sample_tn : msaa_TN) {
180	// Exploit bilinearity: distances are just additions in (T,N) space
181	float d_parallel = pixel_TN.x + sample_tn.x; // Distance along line
182	float d_perp = pixel_TN.y + sample_tn.y; // Perpendicular distance (signed)
183	float d_perp_sq = d_perp * d_perp; // Squared distance (positive)
184
185	// Early rejection: too far perpendicular to line
186	if (d_perp_sq > pen_radius_sq) continue;
187
188	// Compute half-width of integration interval
189	float h = std::sqrt(pen_radius_sq - d_perp_sq);
190
191	// Fast path: middle section (most common case for long lines)
192	// When h < d_parallel < L - h, the integration interval is 2*h
193	// We accumulate h (not 2*h) and compensate with pre-doubled ink_contribution
194	if (h < d_parallel && d_parallel < L - h) {
195	half_ink_sum += h;
196	} else {
197	// Slow path: near start or end caps (needs clamping)
198	float t_min = std::max(0.0f, d_parallel - h);
199	float t_max = std::min(L, d_parallel + h);
200	half_ink_sum += 0.5f * std::max(0.0f, t_max - t_min);
201	}
202	}
203	if (half_ink_sum > 0.0f) {
204	line[x] += ink_contribution_double * half_ink_sum;
205	}
206	};
207
208	const float epsilon = 1e-6f;
209	if (std::abs(T.y) > epsilon) {
210	// Non-horizontal segment: use scanline intersection for tighter x bounds
211	const float inv_Ty = 1.0f / T.y;
212	const float dx_factor = (std::abs(N.x) > epsilon ? pen_radius / std::abs(N.x) : pen_radius);
213
214	auto global_x_min = std::numeric_limits<int>::max();
215	auto global_x_max = std::numeric_limits<int>::min();
216	for (int y = y_min; y <= y_max; ++y) {
217	auto line = canvas.getLine(y);
218	float y_float = float(y);
219	float dy = y_float - A.y; // Relative to A (loop-invariant for x-loop)
220	float Py = y_float + 0.5f;
221	float t_center = (Py - A.y) * inv_Ty;
222	float x_center = A.x + t_center * T.x;
223	auto x_min = std::max(0, (int)std::floor(x_center - dx_factor));
224	auto x_max = std::min(canvas.size().x - 1, (int)std::ceil(x_center + dx_factor));
225	global_x_min = std::min(global_x_min, x_min);
226	global_x_max = std::max(global_x_max, x_max);
227
228	for (int x = x_min; x <= x_max; ++x) {
229	process_pixel(line, x, dy);
230	}
231	}
232	return BoundingBox{
233	.tl = gl::ivec2{global_x_min, y_min},
234	.br = gl::ivec2{global_x_max + 1, y_max + 1}
235	};
236	} else {
237	// Nearly horizontal segment: use axis-aligned bounding box
238	int x_min = std::max(0, int(std::floor(segment_min.x - pen_radius)));
239	int x_max = std::min(canvas.size().x - 1, int(std::ceil(segment_max.x + pen_radius)));
240
241	for (int y = y_min; y <= y_max; ++y) {
242	auto line = canvas.getLine(y);
243	float dy = float(y) - A.y; // Relative to A (loop-invariant for x-loop)
244	for (int x = x_min; x <= x_max; ++x) {
245	process_pixel(line, x, dy);
246	}
247	}
248	return BoundingBox{
249	.tl = gl::ivec2{x_min, y_min},
250	.br = gl::ivec2{x_max + 1, y_max + 1}
251	};
252	}
253	}
254
255	/*
256	// Rasterize a polyline with proper endpoint dwell times
257	template<typename InkType>
258	void rasterize_polyline(
259	Canvas<InkType>& canvas, std::span<const gl::vec2> pts,
260	float pen_radius, InkType ink_color, float dwell_time = 1.0f)
261	{
262	if (pts.size() < 2) {
263	if (pts.size() == 1) {
264	// Single point: just a dot with specified dwell time
265	rasterize_dot(canvas, pts[0], pen_radius, ink_color * dwell_time);
266	}
267	return;
268	}
269
270	// Normalize dwell to ink amount and pre-multiply with ink_color
271	const float ink_amount = dwell_time / (2.0f * pen_radius);
272	const InkType endpoint_ink = ink_color * ink_amount;
273
274	// Dwell at start point (pen lands)
275	rasterize_dot(canvas, pts[0], pen_radius, endpoint_ink);
276
277	// All segments (just motion, no dwell at interior points)
278	for (size_t i = 0; i + 1 < pts.size(); ++i) {
279	rasterize_segment_motion(canvas, pts[i], pts[i + 1], pen_radius, ink_color);
280	}
281
282	// Dwell at end point (pen lifts)
283	rasterize_dot(canvas, pts.back(), pen_radius, endpoint_ink);
284	}
285
286	// Convenience helper. Is this useful enough? Remove?
287	template <typename InkType>
288	void rasterize_segment(
289	Canvas<InkType>& canvas, gl::vec2 A, gl::vec2 B,
290	float pen_radius, InkType ink_color, float dwell_time = 1.0f)
291	{
292	std::array pts = {A, B};
293	rasterize_polyline(canvas, pts, pen_radius, ink_color, dwell_time);
294	}*/
295
296	////
297
298	// InkType: 'gl::vec3' for RGB, 'float' for grayscale
299	// Deduce pixel type from ink type
300	template<typename InkType> struct PixelTypeFor;
301	template<> struct PixelTypeFor<gl::vec3> { using type = RgbPixel; };
302	template<> struct PixelTypeFor<float> { using type = GrayPixel; };
303	template<typename InkType> using PixelTypeFor_t = typename PixelTypeFor<InkType>::type;
304
305	template<typename Sampler> using InkTypeForSampler = decltype(std::declval<Sampler>()(0, 0));
306
307	template<typename Sampler>
308	Canvas<PixelTypeFor_t<InkTypeForSampler<Sampler>>> ink_to_pixels(gl::ivec2 size, Sampler sampler)
309	{
310	using InkType = InkTypeForSampler<Sampler>;
311	using PixelType = PixelTypeFor_t<InkType>;
312
313	// Saturation strength - models paper's limited ink absorption
314	// Rational function: ink(1+k) / (1+kink), scaled to 0-255 range
315	// - Passes through (1,1): single stroke looks correct
316	// - For overlaps (ink>1): output < ink (overlaps saturate smoothly)
317	constexpr float k = 4.0f;
318	constexpr float numerator = (1.0f + k) * 255.0f;
319
320	Canvas<PixelType> output(size, PixelType(255));
321	for (int y = 0; y < size.y; ++y) {
322	auto out_line = output.getLine(y);
323	for (int x = 0; x < size.x; ++x) {
324	InkType accumulated_ink = sampler(x, y);
325	if (accumulated_ink == InkType(0.0f)) {
326	continue; // Already initialized to white (255)
327	}
328
329	// Apply saturation function in 0-255 range (formula works for both float and gl::vec3)
330	InkType saturated_ink_255 = accumulated_ink * numerator / (InkType(1.0f) + k * accumulated_ink);
331
332	// Apply subtractive ink model: white paper (255) minus saturated ink
333	InkType pixel_value = InkType(255.0f) - saturated_ink_255;
334	using gl::clamp; // Bring both clamp functions into scope ..
335	using std::clamp; // .. overload resolution picks the right one
336	out_line[x] = PixelType(clamp(pixel_value, InkType(0.0f), InkType(255.0f)));
337	}
338	}
339	return output;
340	}
341
342	////
343
344	PlotterPaper::PlotterPaper(gl::ivec2 size)
345	: paper(size)
346	, damage{.tl = gl::ivec2{0}, .br = size}
347	{
348	}
349
350	void PlotterPaper::draw_dot(gl::vec2 pos, float radius, gl::vec3 color)
351	{
352	anythingPlotted = true;
353	auto box = rasterize_dot(paper, pos, radius, color);
354	damage = BoundingBox::merge(damage, box);
355	}
356
357	void PlotterPaper::draw_motion(gl::vec2 from, gl::vec2 to, float radius, gl::vec3 color)
358	{
359	anythingPlotted = true;
360	auto box = rasterize_segment_motion(paper, from, to, radius, color);
361	damage = BoundingBox::merge(damage, box);
362	}
363
364	Canvas<RgbPixel> PlotterPaper::getRGB() const
365	{
366	auto out_size = size();
367	auto sampler = [&](int x, int y) { return paper.getLine(y)[x]; };
368	return ink_to_pixels(out_size, sampler);
369	}
370
371	gl::Texture& PlotterPaper::updateTexture(gl::ivec2 targetSize)
372	{
373	assert(targetSize.x >= 1);
374	assert(targetSize.y >= 1);
375
376	auto fullSize = paper.size();
377	auto clampedTargetSize = min(targetSize, gl::ivec2(1024)); // ensure below 2048 x 2028
378	auto factor = max_component(trunc(fullSize / clampedTargetSize));
379	auto newTexSize = max(fullSize / factor, gl::ivec2(1));
380	bool sizeChanged = newTexSize != texSize;
381
382	if (texture.get() && !sizeChanged && damage.empty()) {
383	return texture; // already up to date
384	}
385	if (sizeChanged) {
386	// full update if size changed
387	damage = BoundingBox{.tl = gl::ivec2{0}, .br = paper.size()};
388	}
389	auto dstTl = damage.tl / factor;
390	auto dstBr = damage.br / factor;
391	auto dstSize = max(dstBr - dstTl, gl::ivec2(1));
392	auto srcOffset = dstTl * factor;
393
394	auto sample_direct = [&](int x, int y) { return paper.getLine(y + srcOffset.y)[x + srcOffset.x]; };
395	auto down_sample_generic = [&](int x, int y) { // run-time 'factor'
396	auto srcPos = gl::ivec2(x, y) * factor + srcOffset;
397	auto sum = gl::vec3(0.0f);
398	auto factorXY = min(paper.size(), gl::ivec2(factor));
399	for (int sy = 0; sy < factorXY.y; ++sy) {
400	auto line = subspan(paper.getLine(srcPos.y + sy), srcPos.x, factorXY.x);
401	for (int sx = 0; sx < factorXY.x; ++sx) {
402	sum += line[sx];
403	}
404	}
405	return sum * (1.0f / float(factor * factor));
406	};
407	auto down_sample_N = [&]<int FACTOR>(int x, int y) { // compile-time 'FACTOR'
408	auto srcPos = gl::ivec2(x, y) * factor + srcOffset;
409	auto sum = gl::vec3(0.0f);
410	for (int sy = 0; sy < FACTOR; ++sy) {
411	auto line = subspan<FACTOR>(paper.getLine(srcPos.y + sy), srcPos.x);
412	for (int sx = 0; sx < FACTOR; ++sx) {
413	sum += line[sx];
414	}
415	}
416	return sum * (1.0f / float(FACTOR * FACTOR));
417	};
418
419
420	auto rgb = [&]{
421	// specializations for small factors, this gives a reasonably large speedup
422	switch (factor) {
423	case 1: return ink_to_pixels(dstSize, sample_direct);
424	case 2: return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<2>(x, y); });
425	case 3: return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<3>(x, y); });
426	case 4: return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<4>(x, y); });
427	case 5: return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<5>(x, y); });
428	case 6: return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<6>(x, y); });
429	case 7: return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<7>(x, y); });
430	case 8: return ink_to_pixels(dstSize, [&](int x, int y) { return down_sample_N.operator()<8>(x, y); });
431	default: return ink_to_pixels(dstSize, down_sample_generic);
432	}
433	}();
434
435	if (!texture.get()) {
436	texture.allocate();
437	texture.setInterpolation(true);
438	} else {
439	texture.bind();
440	}
441	if (sizeChanged) {
442	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, newTexSize.x, newTexSize.y, 0, GL_RGB,
443	GL_UNSIGNED_BYTE, &rgb.getLine(0)[0].x);
444	} else {
445	glTexSubImage2D(GL_TEXTURE_2D, 0, dstTl.x, dstTl.y, dstSize.x, dstSize.y, GL_RGB,
446	GL_UNSIGNED_BYTE, &rgb.getLine(0)[0].x);
447	}
448
449	damage.clear();
450	texSize = newTexSize;
451	return texture;
452	}
453
454	} // namespace openmsx

Keyboard