more work on crt-royale >_<

blurs/blur9fast-horizontal.slang

@@ -0,0 +1,98 @@
+#version 450
+
+layout(push_constant) uniform Push
+{
+	vec4 SourceSize;
+	vec4 OriginalSize;
+	vec4 OutputSize;
+	uint FrameCount;
+} params;
+
+layout(std140, set = 0, binding = 0) uniform UBO
+{
+	mat4 MVP;
+} global;
+
+/////////////////////////////////  MIT LICENSE  ////////////////////////////////
+
+//  Copyright (C) 2014 TroggleMonkey
+//
+//  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.
+
+
+/////////////////////////////  SETTINGS MANAGEMENT  ////////////////////////////
+
+//  PASS SETTINGS:
+//  gamma-management.h needs to know what kind of pipeline we're using and
+//  what pass this is in that pipeline.  This will become obsolete if/when we
+//  can #define things like this in the .cgp preset file.
+//#define GAMMA_ENCODE_EVERY_FBO
+//#define FIRST_PASS
+//#define LAST_PASS
+//#define SIMULATE_CRT_ON_LCD
+//#define SIMULATE_GBA_ON_LCD
+//#define SIMULATE_LCD_ON_CRT
+//#define SIMULATE_GBA_ON_CRT
+
+
+//////////////////////////////////  INCLUDES  //////////////////////////////////
+
+//  #included by vertex shader:
+//#include "../include/gamma-management.h"
+//#include "../include/blur-functions.h"
+#include "../crt/shaders/crt-royale/src/includes.h"
+
+#pragma stage vertex
+layout(location = 0) in vec4 Position;
+layout(location = 1) in vec2 TexCoord;
+layout(location = 0) out vec2 tex_uv;
+layout(location = 1) out vec2 blur_dxdy;
+
+void main()
+{
+   gl_Position = global.MVP * Position;
+   tex_uv = TexCoord;
+   
+   	//  Get the uv sample distance between output pixels.  Blurs are not generic
+    //  Gaussian resizers, and correct blurs require:
+    //  1.) params.OutputSize.xy == params.SourceSize.xy * 2^m, where m is an integer <= 0.
+    //  2.) mipmap_inputN = "true" for this pass in preset if m != 0
+    //  3.) filter_linearN = "true" except for 1x scale nearest neighbor blurs
+    //  Gaussian resizers would upsize using the distance between input texels
+    //  (not output pixels), but we avoid this and consistently blur at the
+    //  destination size.  Otherwise, combining statically calculated weights
+    //  with bilinear sample exploitation would result in terrible artifacts.
+    const vec2 dxdy_scale = params.SourceSize.xy * params.OutputSize.zw;
+	const vec2 dxdy = dxdy_scale * params.SourceSize.zw;
+    //  This blur is horizontal-only, so zero out the vertical offset:
+	blur_dxdy = vec2(dxdy.x, 0.0);
+}
+
+#pragma stage fragment
+layout(location = 0) in vec2 tex_uv;
+layout(location = 1) in vec2 blur_dxdy;
+layout(location = 0) out vec4 FragColor;
+layout(set = 0, binding = 2) uniform sampler2D Source;
+
+void main()
+{
+	vec3 color = tex2Dblur9fast(Source, tex_uv, blur_dxdy);
+    //  Encode and output the blurred image:
+   FragColor = vec4(1.0);//encode_output(vec4(color, 1.0));
+}

blurs/blur9fast-vertical.slang

@@ -0,0 +1,98 @@
+#version 450
+
+layout(push_constant) uniform Push
+{
+	vec4 SourceSize;
+	vec4 OriginalSize;
+	vec4 OutputSize;
+	uint FrameCount;
+} params;
+
+layout(std140, set = 0, binding = 0) uniform UBO
+{
+	mat4 MVP;
+} global;
+
+/////////////////////////////////  MIT LICENSE  ////////////////////////////////
+
+//  Copyright (C) 2014 TroggleMonkey
+//
+//  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.
+
+
+/////////////////////////////  SETTINGS MANAGEMENT  ////////////////////////////
+
+//  PASS SETTINGS:
+//  gamma-management.h needs to know what kind of pipeline we're using and
+//  what pass this is in that pipeline.  This will become obsolete if/when we
+//  can #define things like this in the .cgp preset file.
+//#define GAMMA_ENCODE_EVERY_FBO
+//#define FIRST_PASS
+//#define LAST_PASS
+//#define SIMULATE_CRT_ON_LCD
+//#define SIMULATE_GBA_ON_LCD
+//#define SIMULATE_LCD_ON_CRT
+//#define SIMULATE_GBA_ON_CRT
+
+
+//////////////////////////////////  INCLUDES  //////////////////////////////////
+
+//  #included by vertex shader:
+//#include "../include/gamma-management.h"
+//#include "../include/blur-functions.h"
+#include "../crt/shaders/crt-royale/src/includes.h"
+
+#pragma stage vertex
+layout(location = 0) in vec4 Position;
+layout(location = 1) in vec2 TexCoord;
+layout(location = 0) out vec2 tex_uv;
+layout(location = 1) out vec2 blur_dxdy;
+
+void main()
+{
+   gl_Position = global.MVP * Position;
+   tex_uv = TexCoord;
+   
+   	//  Get the uv sample distance between output pixels.  Blurs are not generic
+    //  Gaussian resizers, and correct blurs require:
+    //  1.) params.OutputSize.xy == params.SourceSize.xy * 2^m, where m is an integer <= 0.
+    //  2.) mipmap_inputN = "true" for this pass in preset if m != 0
+    //  3.) filter_linearN = "true" except for 1x scale nearest neighbor blurs
+    //  Gaussian resizers would upsize using the distance between input texels
+    //  (not output pixels), but we avoid this and consistently blur at the
+    //  destination size.  Otherwise, combining statically calculated weights
+    //  with bilinear sample exploitation would result in terrible artifacts.
+    const vec2 dxdy_scale = params.SourceSize.xy * params.OutputSize.zw;
+	const vec2 dxdy = dxdy_scale * params.SourceSize.zw;
+    //  This blur is vertical-only, so zero out the vertical offset:
+	blur_dxdy = vec2(0.0, dxdy.y);
+}
+
+#pragma stage fragment
+layout(location = 0) in vec2 tex_uv;
+layout(location = 1) in vec2 blur_dxdy;
+layout(location = 0) out vec4 FragColor;
+layout(set = 0, binding = 2) uniform sampler2D Source;
+
+void main()
+{
+	vec3 color = tex2Dblur9fast(Source, tex_uv, blur_dxdy);
+    //  Encode and output the blurred image:
+   FragColor = vec4(1.0);//encode_output(vec4(color, 1.0));
+}

crt/crt-royale-test.slangp

@@ -0,0 +1,92 @@
+# IMPORTANT:
+# Shader passes need to know details about the image in the mask_texture LUT
+# files, so set the following constants in user-preset-constants.h accordingly:
+# 1.) mask_triads_per_tile = (number of horizontal triads in mask texture LUT's)
+# 2.) mask_texture_small_size = (texture size of mask*texture_small LUT's)
+# 3.) mask_texture_large_size = (texture size of mask*texture_large LUT's)
+# 4.) mask_grille_avg_color = (avg. brightness of mask_grille_texture* LUT's, in [0, 1])
+# 5.) mask_slot_avg_color = (avg. brightness of mask_slot_texture* LUT's, in [0, 1])
+# 6.) mask_shadow_avg_color = (avg. brightness of mask_shadow_texture* LUT's, in [0, 1])
+# Shader passes also need to know certain scales set in this preset, but their
+# compilation model doesn't currently allow the preset file to tell them.  Make
+# sure to set the following constants in user-preset-constants.h accordingly too:
+# 1.) bloom_approx_scale_x = scale_x2
+# 2.) mask_resize_viewport_scale = vec2(scale_x6, scale_y5)
+# Finally, shader passes need to know the value of geom_max_aspect_ratio used to
+# calculate scale_y5 (among other values):
+# 1.) geom_max_aspect_ratio = (geom_max_aspect_ratio used to calculate scale_y5)
+
+shaders = "1"//"12"
+
+# Set an identifier, filename, and sampling traits for the phosphor mask texture.
+# Load an aperture grille, slot mask, and an EDP shadow mask, and load a small
+# non-mipmapped version and a large mipmapped version.
+# TODO: Test masks in other directories.
+textures = "mask_grille_texture_small;mask_grille_texture_large;mask_slot_texture_small;mask_slot_texture_large;mask_shadow_texture_small;mask_shadow_texture_large"
+mask_grille_texture_small = "shaders/crt-royale/TileableLinearApertureGrille15Wide8And5d5SpacingResizeTo64.png"
+mask_grille_texture_large = "shaders/crt-royale/TileableLinearApertureGrille15Wide8And5d5Spacing.png"
+mask_slot_texture_small = "shaders/crt-royale/TileableLinearSlotMaskTall15Wide9And4d5Horizontal9d14VerticalSpacingResizeTo64.png"
+mask_slot_texture_large = "shaders/crt-royale/TileableLinearSlotMaskTall15Wide9And4d5Horizontal9d14VerticalSpacing.png"
+mask_shadow_texture_small = "shaders/crt-royale/TileableLinearShadowMaskEDPResizeTo64.png"
+mask_shadow_texture_large = "shaders/crt-royale/TileableLinearShadowMaskEDP.png"
+mask_grille_texture_small_wrap_mode = "repeat"
+mask_grille_texture_large_wrap_mode = "repeat"
+mask_slot_texture_small_wrap_mode = "repeat"
+mask_slot_texture_large_wrap_mode = "repeat"
+mask_shadow_texture_small_wrap_mode = "repeat"
+mask_shadow_texture_large_wrap_mode = "repeat"
+mask_grille_texture_small_linear = "true"
+mask_grille_texture_large_linear = "true"
+mask_slot_texture_small_linear = "true"
+mask_slot_texture_large_linear = "true"
+mask_shadow_texture_small_linear = "true"
+mask_shadow_texture_large_linear = "true"
+mask_grille_texture_small_mipmap = "false"  # Mipmapping causes artifacts with manually resized masks without tex2Dlod
+mask_grille_texture_large_mipmap = "true"   # Essential for hardware-resized masks
+mask_slot_texture_small_mipmap = "false"    # Mipmapping causes artifacts with manually resized masks without tex2Dlod
+mask_slot_texture_large_mipmap = "true"     # Essential for hardware-resized masks
+mask_shadow_texture_small_mipmap = "false"  # Mipmapping causes artifacts with manually resized masks without tex2Dlod
+mask_shadow_texture_large_mipmap = "true"   # Essential for hardware-resized masks
+
+
+# Pass5: Lanczos-resize the phosphor mask vertically.  Set the absolute
+# scale_x5 == mask_texture_small_size.x (see IMPORTANT above).  Larger scales
+# will blur, and smaller scales could get nasty.  The vertical size must be
+# based on the viewport size and calculated carefully to avoid artifacts later.
+# First calculate the minimum number of mask tiles we need to draw.
+# Since curvature is computed after the scanline masking pass:
+#   num_resized_mask_tiles = 2.0;
+# If curvature were computed in the scanline masking pass (it's not):
+#   max_mask_texel_border = ~3.0 * (1/3.0 + 4.0*sqrt(2.0) + 0.5 + 1.0);
+#   max_mask_tile_border = max_mask_texel_border/
+#       (min_resized_phosphor_triad_size * mask_triads_per_tile);
+#   num_resized_mask_tiles = max(2.0, 1.0 + max_mask_tile_border * 2.0);
+#   At typical values (triad_size >= 2.0, mask_triads_per_tile == 8):
+#       num_resized_mask_tiles = ~3.8
+# Triad sizes are given in horizontal terms, so we need geom_max_aspect_ratio
+# to relate them to vertical resolution.  The widest we expect is:
+#   geom_max_aspect_ratio = 4.0/3.0  # Note: Shader passes need to know this!
+# The fewer triads we tile across the screen, the larger each triad will be as a
+# fraction of the viewport size, and the larger scale_y5 must be to draw a full
+# num_resized_mask_tiles.  Therefore, we must decide the smallest number of
+# triads we'll guarantee can be displayed on screen.  We'll set this according
+# to 3-pixel triads at 768p resolution (the lowest anyone's likely to use):
+#   min_allowed_viewport_triads = 768.0*geom_max_aspect_ratio / 3.0 = 341.333333
+# Now calculate the viewport scale that ensures we can draw resized_mask_tiles:
+#   min_scale_x = resized_mask_tiles * mask_triads_per_tile /
+#       min_allowed_viewport_triads
+#   scale_y5 = geom_max_aspect_ratio * min_scale_x
+#   # Some code might depend on equal scales:
+#   scale_x6 = scale_y5
+# Given our default geom_max_aspect_ratio and min_allowed_viewport_triads:
+#   scale_y5 = 4.0/3.0 * 2.0/(341.33333 / 8.0) = 0.0625
+# IMPORTANT: The scales MUST be calculated in this way.  If you wish to change
+# geom_max_aspect_ratio, update that constant in user-preset-constants.h!
+shader0 = "shaders/crt-royale/src/crt-royale-mask-resize-vertical.slang"
+filter_linear0 = "true"
+scale_type_x0 = "absolute"
+scale_x0 = "64"
+scale_type_y0 = "viewport"
+scale_y0 = "0.0625" # Safe for >= 341.333 horizontal triads at viewport size
+#srgb_framebuffer0 = "false" # mask_texture is already assumed linear
+

crt/crt-royale.slangp

@@ -0,0 +1,206 @@
+# IMPORTANT:
+# Shader passes need to know details about the image in the mask_texture LUT
+# files, so set the following constants in user-preset-constants.h accordingly:
+# 1.) mask_triads_per_tile = (number of horizontal triads in mask texture LUT's)
+# 2.) mask_texture_small_size = (texture size of mask*texture_small LUT's)
+# 3.) mask_texture_large_size = (texture size of mask*texture_large LUT's)
+# 4.) mask_grille_avg_color = (avg. brightness of mask_grille_texture* LUT's, in [0, 1])
+# 5.) mask_slot_avg_color = (avg. brightness of mask_slot_texture* LUT's, in [0, 1])
+# 6.) mask_shadow_avg_color = (avg. brightness of mask_shadow_texture* LUT's, in [0, 1])
+# Shader passes also need to know certain scales set in this preset, but their
+# compilation model doesn't currently allow the preset file to tell them.  Make
+# sure to set the following constants in user-preset-constants.h accordingly too:
+# 1.) bloom_approx_scale_x = scale_x2
+# 2.) mask_resize_viewport_scale = vec2(scale_x6, scale_y5)
+# Finally, shader passes need to know the value of geom_max_aspect_ratio used to
+# calculate scale_y5 (among other values):
+# 1.) geom_max_aspect_ratio = (geom_max_aspect_ratio used to calculate scale_y5)
+
+shaders = "5"//"12"
+
+# Set an identifier, filename, and sampling traits for the phosphor mask texture.
+# Load an aperture grille, slot mask, and an EDP shadow mask, and load a small
+# non-mipmapped version and a large mipmapped version.
+# TODO: Test masks in other directories.
+textures = "mask_grille_texture_small;mask_grille_texture_large;mask_slot_texture_small;mask_slot_texture_large;mask_shadow_texture_small;mask_shadow_texture_large"
+mask_grille_texture_small = "shaders/crt-royale/TileableLinearApertureGrille15Wide8And5d5SpacingResizeTo64.png"
+mask_grille_texture_large = "shaders/crt-royale/TileableLinearApertureGrille15Wide8And5d5Spacing.png"
+mask_slot_texture_small = "shaders/crt-royale/TileableLinearSlotMaskTall15Wide9And4d5Horizontal9d14VerticalSpacingResizeTo64.png"
+mask_slot_texture_large = "shaders/crt-royale/TileableLinearSlotMaskTall15Wide9And4d5Horizontal9d14VerticalSpacing.png"
+mask_shadow_texture_small = "shaders/crt-royale/TileableLinearShadowMaskEDPResizeTo64.png"
+mask_shadow_texture_large = "shaders/crt-royale/TileableLinearShadowMaskEDP.png"
+mask_grille_texture_small_wrap_mode = "repeat"
+mask_grille_texture_large_wrap_mode = "repeat"
+mask_slot_texture_small_wrap_mode = "repeat"
+mask_slot_texture_large_wrap_mode = "repeat"
+mask_shadow_texture_small_wrap_mode = "repeat"
+mask_shadow_texture_large_wrap_mode = "repeat"
+mask_grille_texture_small_linear = "true"
+mask_grille_texture_large_linear = "true"
+mask_slot_texture_small_linear = "true"
+mask_slot_texture_large_linear = "true"
+mask_shadow_texture_small_linear = "true"
+mask_shadow_texture_large_linear = "true"
+mask_grille_texture_small_mipmap = "false"  # Mipmapping causes artifacts with manually resized masks without tex2Dlod
+mask_grille_texture_large_mipmap = "true"   # Essential for hardware-resized masks
+mask_slot_texture_small_mipmap = "false"    # Mipmapping causes artifacts with manually resized masks without tex2Dlod
+mask_slot_texture_large_mipmap = "true"     # Essential for hardware-resized masks
+mask_shadow_texture_small_mipmap = "false"  # Mipmapping causes artifacts with manually resized masks without tex2Dlod
+mask_shadow_texture_large_mipmap = "true"   # Essential for hardware-resized masks
+
+
+# Pass0: Linearize the input based on CRT gamma and bob interlaced fields.
+# (Bobbing ensures we can immediately blur without getting artifacts.)
+shader0 = "shaders/crt-royale/src/crt-royale-first-pass-linearize-crt-gamma-bob-fields.slang"
+alias0 = "ORIG_LINEARIZED"
+filter_linear0 = "false"
+scale_type0 = "source"
+scale0 = "1.0"
+srgb_framebuffer0 = "true"
+
+# Pass1: Resample interlaced (and misconverged) scanlines vertically.
+# Separating vertical/horizontal scanline sampling is faster: It lets us
+# consider more scanlines while calculating weights for fewer pixels, and
+# it reduces our samples from vertical*horizontal to vertical+horizontal.
+# This has to come right after ORIG_LINEARIZED, because there's no
+# "original_source" scale_type we can use later.
+shader1 = "shaders/crt-royale/src/crt-royale-scanlines-vertical-interlacing.slang"
+alias1 = "VERTICAL_SCANLINES"
+filter_linear1 = "true"
+scale_type_x1 = "source"
+scale_x1 = "1.0"
+scale_type_y1 = "viewport"
+scale_y1 = "1.0"
+#float_framebuffer1 = "true"
+srgb_framebuffer1 = "true"
+
+# Pass2: Do a small resize blur of ORIG_LINEARIZED at an absolute size, and
+# account for convergence offsets.  We want to blur a predictable portion of the
+# screen to match the phosphor bloom, and absolute scale works best for
+# reliable results with a fixed-size bloom.  Picking a scale is tricky:
+# a.) 400x300 is a good compromise for the "fake-bloom" version: It's low enough
+#     to blur high-res/interlaced sources but high enough that resampling
+#     doesn't smear low-res sources too much.
+# b.) 320x240 works well for the "real bloom" version: It's 1-1.5% faster, and
+#     the only noticeable visual difference is a larger halation spread (which
+#     may be a good thing for people who like to crank it up).
+# Note the 4:3 aspect ratio assumes the input has cropped geom_overscan (so it's
+# *intended* for an ~4:3 aspect ratio).
+shader2 = "shaders/crt-royale/src/crt-royale-bloom-approx.slang"
+alias2 = "BLOOM_APPROX"
+filter_linear2 = "true"
+scale_type2 = "absolute"
+scale_x2 = "320"
+scale_y2 = "240"
+srgb_framebuffer2 = "true"
+
+# Pass3: Vertically blur the input for halation and refractive diffusion.
+# Base this on BLOOM_APPROX: This blur should be small and fast, and blurring
+# a constant portion of the screen is probably physically correct if the
+# viewport resolution is proportional to the simulated CRT size.
+shader3 = "../blurs/blur9fast-vertical.slang"
+filter_linear3 = "true"
+scale_type3 = "source"
+scale3 = "1.0"
+srgb_framebuffer3 = "true"
+
+# Pass4: Horizontally blur the input for halation and refractive diffusion.
+# Note: Using a one-pass 9x9 blur is about 1% slower.
+shader4 = "../blurs/blur9fast-horizontal.slang"
+alias4 = "HALATION_BLUR"
+filter_linear4 = "true"
+scale_type4 = "source"
+scale4 = "1.0"
+srgb_framebuffer4 = "true"
+
+# Pass5: Lanczos-resize the phosphor mask vertically.  Set the absolute
+# scale_x5 == mask_texture_small_size.x (see IMPORTANT above).  Larger scales
+# will blur, and smaller scales could get nasty.  The vertical size must be
+# based on the viewport size and calculated carefully to avoid artifacts later.
+# First calculate the minimum number of mask tiles we need to draw.
+# Since curvature is computed after the scanline masking pass:
+#   num_resized_mask_tiles = 2.0;
+# If curvature were computed in the scanline masking pass (it's not):
+#   max_mask_texel_border = ~3.0 * (1/3.0 + 4.0*sqrt(2.0) + 0.5 + 1.0);
+#   max_mask_tile_border = max_mask_texel_border/
+#       (min_resized_phosphor_triad_size * mask_triads_per_tile);
+#   num_resized_mask_tiles = max(2.0, 1.0 + max_mask_tile_border * 2.0);
+#   At typical values (triad_size >= 2.0, mask_triads_per_tile == 8):
+#       num_resized_mask_tiles = ~3.8
+# Triad sizes are given in horizontal terms, so we need geom_max_aspect_ratio
+# to relate them to vertical resolution.  The widest we expect is:
+#   geom_max_aspect_ratio = 4.0/3.0  # Note: Shader passes need to know this!
+# The fewer triads we tile across the screen, the larger each triad will be as a
+# fraction of the viewport size, and the larger scale_y5 must be to draw a full
+# num_resized_mask_tiles.  Therefore, we must decide the smallest number of
+# triads we'll guarantee can be displayed on screen.  We'll set this according
+# to 3-pixel triads at 768p resolution (the lowest anyone's likely to use):
+#   min_allowed_viewport_triads = 768.0*geom_max_aspect_ratio / 3.0 = 341.333333
+# Now calculate the viewport scale that ensures we can draw resized_mask_tiles:
+#   min_scale_x = resized_mask_tiles * mask_triads_per_tile /
+#       min_allowed_viewport_triads
+#   scale_y5 = geom_max_aspect_ratio * min_scale_x
+#   # Some code might depend on equal scales:
+#   scale_x6 = scale_y5
+# Given our default geom_max_aspect_ratio and min_allowed_viewport_triads:
+#   scale_y5 = 4.0/3.0 * 2.0/(341.33333 / 8.0) = 0.0625
+# IMPORTANT: The scales MUST be calculated in this way.  If you wish to change
+# geom_max_aspect_ratio, update that constant in user-preset-constants.h!
+shader5 = "shaders/crt-royale/src/crt-royale-mask-resize-vertical.slang"
+filter_linear5 = "true"
+scale_type_x5 = "absolute"
+scale_x5 = "64"
+scale_type_y5 = "viewport"
+scale_y5 = "0.0625" # Safe for >= 341.333 horizontal triads at viewport size
+#srgb_framebuffer5 = "false" # mask_texture is already assumed linear
+
+# Pass6: Lanczos-resize the phosphor mask horizontally.  scale_x6 = scale_y5.
+# TODO: Check again if the shaders actually require equal scales.
+shader6 = "shaders/crt-royale/src/crt-royale-mask-resize-horizontal.slang"
+alias6 = "MASK_RESIZE"
+filter_linear6 = "false"
+scale_type_x6 = "viewport"
+scale_x6 = "0.0625"
+scale_type_y6 = "source"
+scale_y6 = "1.0"
+#srgb_framebuffer6 = "false" # mask_texture is already assumed linear
+
+# Pass7: Resample (misconverged) scanlines horizontally, apply halation, and
+# apply the phosphor mask.
+shader7 = "shaders/crt-royale/src/crt-royale-scanlines-horizontal-apply-mask.slang"
+alias7 = "MASKED_SCANLINES"
+filter_linear7 = "true" # This could just as easily be nearest neighbor.
+scale_type7 = "viewport"
+scale7 = "1.0"
+#float_framebuffer7 = "true"
+srgb_framebuffer7 = "true"
+
+# Pass 8: Compute a brightpass.  This will require reading the final mask.
+shader8 = "shaders/crt-royale/src/crt-royale-brightpass.slang"
+alias8 = "BRIGHTPASS"
+filter_linear8 = "true" # This could just as easily be nearest neighbor.
+scale_type8 = "viewport"
+scale8 = "1.0"
+srgb_framebuffer8 = "true"
+
+# Pass 9: Blur the brightpass vertically
+shader9 = "shaders/crt-royale/src/crt-royale-bloom-vertical.slang"
+filter_linear9 = "true" # This could just as easily be nearest neighbor.
+scale_type9 = "source"
+scale9 = "1.0"
+srgb_framebuffer9 = "true"
+
+# Pass 10: Blur the brightpass horizontally and combine it with the dimpass:
+shader10 = "shaders/crt-royale/src/crt-royale-bloom-horizontal-reconstitute.slang"
+filter_linear10 = "true"
+scale_type10 = "source"
+scale10 = "1.0"
+srgb_framebuffer10 = "true"
+
+# Pass 11: Compute curvature/AA:
+shader11 = "shaders/crt-royale/src/crt-royale-geometry-aa-last-pass.slang"
+filter_linear11 = "true"
+scale_type11 = "viewport"
+mipmap_input11 = "true"
+texture_wrap_mode11 = "clamp_to_edge"
+

crt/shaders/crt-royale/src/blur-functions.h

@@ -1,4 +1,5 @@
-#define BLUR_FUNCTIONS
+#ifndef BLUR_FUNCTIONS_H
+#define BLUR_FUNCTIONS_H
 
 /////////////////////////////////  MIT LICENSE  ////////////////////////////////
 
@@ -278,4 +279,42 @@ vec3 tex2Dblur3x3resize(const sampler2D texture, const vec2 tex_uv,
     const vec2 dxdy)
 {
     return tex2Dblur3x3resize(texture, tex_uv, dxdy, blur3_std_dev);
-}
+}
+
+vec3 tex2Dblur9fast(const sampler2D tex, const vec2 tex_uv,
+    const vec2 dxdy, const float sigma)
+{
+    //  Requires:   Same as tex2Dblur11()
+    //  Returns:    A 1D 9x Gaussian blurred texture lookup using 1 nearest
+    //              neighbor and 4 linear taps.  It may be mipmapped depending
+    //              on settings and dxdy.
+    //  First get the texel weights and normalization factor as above.
+    const float denom_inv = 0.5/(sigma*sigma);
+    const float w0 = 1.0;
+    const float w1 = exp(-1.0 * denom_inv);
+    const float w2 = exp(-4.0 * denom_inv);
+    const float w3 = exp(-9.0 * denom_inv);
+    const float w4 = exp(-16.0 * denom_inv);
+    const float weight_sum_inv = 1.0 / (w0 + 2.0 * (w1 + w2 + w3 + w4));
+    //  Calculate combined weights and linear sample ratios between texel pairs.
+    const float w12 = w1 + w2;
+    const float w34 = w3 + w4;
+    const float w12_ratio = w2/w12;
+    const float w34_ratio = w4/w34;
+    //  Statically normalize weights, sum weighted samples, and return:
+    vec3 sum = vec3(0.0);
+    sum += w34 * tex2D_linearize(tex, tex_uv - (3.0 + w34_ratio) * dxdy).rgb;
+    sum += w12 * tex2D_linearize(tex, tex_uv - (1.0 + w12_ratio) * dxdy).rgb;
+    sum += w0 * tex2D_linearize(tex, tex_uv).rgb;
+    sum += w12 * tex2D_linearize(tex, tex_uv + (1.0 + w12_ratio) * dxdy).rgb;
+    sum += w34 * tex2D_linearize(tex, tex_uv + (3.0 + w34_ratio) * dxdy).rgb;
+    return sum * weight_sum_inv;
+}
+
+vec3 tex2Dblur9fast(const sampler2D tex, const vec2 tex_uv,
+    const vec2 dxdy)
+{
+    return tex2Dblur9fast(tex, tex_uv, dxdy, blur9_std_dev);
+}
+
+#endif  //  BLUR_FUNCTIONS_H

crt/shaders/crt-royale/src/crt-royale-bloom-approx.slang

@@ -345,5 +345,5 @@ void main()
         color = vec3(color_r.r, color_g.g, color_b.b);
     }
     //  Encode and output the blurred image:
-   FragColor = vec4(color, 1.0);
+   FragColor = vec4(1.0);//vec4(color, 1.0);
 }

crt/shaders/crt-royale/src/crt-royale-mask-resize-vertical.slang

@@ -0,0 +1,206 @@
+#version 450
+
+layout(push_constant) uniform Push
+{
+	vec4 SourceSize;
+	vec4 OriginalSize;
+	vec4 OutputSize;
+	uint FrameCount;
+} params;
+
+layout(std140, set = 0, binding = 0) uniform UBO
+{
+	mat4 MVP;
+} global;
+
+/////////////////////////////  GPL LICENSE NOTICE  /////////////////////////////
+
+//  crt-royale: A full-featured CRT shader, with cheese.
+//  Copyright (C) 2014 TroggleMonkey <trogglemonkey@gmx.com>
+//
+//  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 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., 59 Temple
+//  Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+/////////////////////////////  SETTINGS MANAGEMENT  ////////////////////////////
+
+//#include "../user-settings.h"
+//#include "derived-settings-and-constants.h"
+//#include "bind-shader-params.h"
+
+//////////////////////////////////  INCLUDES  //////////////////////////////////
+
+#include "includes.h"
+#include "phosphor-mask-resizing.h"
+
+#pragma stage vertex
+layout(location = 0) in vec4 Position;
+layout(location = 1) in vec2 TexCoord;
+layout(location = 0) out vec2 tex_uv;
+layout(location = 1) out vec2 src_tex_uv_wrap;
+layout(location = 2) out vec2 resize_magnification_scale;
+
+void main()
+{
+   gl_Position = global.MVP * Position;
+   tex_uv = TexCoord;
+   
+    //  First estimate the viewport size (the user will get the wrong number of
+    //  triads if it's wrong and mask_specify_num_triads is 1.0/true).
+    const float viewport_y = params.OutputSize.y / mask_resize_viewport_scale.y;
+    const float aspect_ratio = geom_aspect_ratio_x / geom_aspect_ratio_y;
+    const vec2 estimated_viewport_size =
+        vec2(viewport_y * aspect_ratio, viewport_y);
+    //  Estimate the output size of MASK_RESIZE (the next pass).  The estimated
+    //  x component shouldn't matter, because we're not using the x result, and
+    //  we're not swearing it's correct (if we did, the x result would influence
+    //  the y result to maintain the tile aspect ratio).
+    const vec2 estimated_mask_resize_output_size =
+        vec2(params.OutputSize.y * aspect_ratio, params.OutputSize.y);
+    //  Find the final intended [y] size of our resized phosphor mask tiles,
+    //  then the tile size for the current pass (resize y only):
+    const vec2 mask_resize_tile_size = get_resized_mask_tile_size(
+        estimated_viewport_size, estimated_mask_resize_output_size, false);
+    const vec2 pass_output_tile_size = vec2(min(
+        mask_resize_src_lut_size.x, params.OutputSize.x), mask_resize_tile_size.y);
+
+    //  We'll render resized tiles until filling the output FBO or meeting a
+    //  limit, so compute [wrapped] tile uv coords based on the output uv coords
+    //  and the number of tiles that will fit in the FBO.
+    const vec2 output_tiles_this_pass = params.OutputSize.xy / pass_output_tile_size;
+    const vec2 output_video_uv = tex_uv;
+    const vec2 tile_uv_wrap = output_video_uv * output_tiles_this_pass;
+
+    //  The input LUT is just a single mask tile, so texture uv coords are the
+    //  same as tile uv coords (save frac() for the fragment shader).  The
+    //  magnification scale is also straightforward:
+    src_tex_uv_wrap = tile_uv_wrap;
+    resize_magnification_scale =
+        pass_output_tile_size / mask_resize_src_lut_size;
+}
+
+#pragma stage fragment
+layout(location = 0) in vec2 tex_uv;
+layout(location = 1) in vec2 src_tex_uv_wrap;
+layout(location = 2) in vec2 resize_magnification_scale;
+layout(location = 0) out vec4 FragColor;
+layout(set = 0, binding = 2) uniform sampler2D Source;
+#ifdef PHOSPHOR_MASK_RESIZE_MIPMAPPED_LUT
+layout(set = 0, binding = 3) uniform sampler2D mask_grille_texture_large;
+layout(set = 0, binding = 4) uniform sampler2D mask_slot_texture_large;
+layout(set = 0, binding = 5) uniform sampler2D mask_shadow_texture_large;
+
+void main()
+{
+    //  Resize the input phosphor mask tile to the final vertical size it will
+    //  appear on screen.  Keep 1x horizontal size if possible (IN.output_size
+    //  >= mask_resize_src_lut_size), and otherwise linearly sample horizontally
+    //  to fit exactly one tile.  Lanczos-resizing the phosphor mask achieves
+    //  much sharper results than mipmapping, and vertically resizing first
+    //  minimizes the total number of taps required.  We output a number of
+    //  resized tiles >= mask_resize_num_tiles for easier tiled sampling later.
+    #ifdef PHOSPHOR_MASK_MANUALLY_RESIZE
+        //  Discard unneeded fragments in case our profile allows real branches.
+        const vec2 tile_uv_wrap = src_tex_uv_wrap;
+        if(get_mask_sample_mode() < 0.5 &&
+            tile_uv_wrap.y <= mask_resize_num_tiles)
+        {
+            const float src_dy = 1.0/mask_resize_src_lut_size.y;
+            const vec2 src_tex_uv = frac(src_tex_uv_wrap);
+            vec3 pixel_color;
+            //  If mask_type is static, this branch will be resolved statically.
+            if(mask_type < 0.5)
+            {
+                pixel_color = downsample_vertical_sinc_tiled(
+                    mask_grille_texture_large, src_tex_uv, mask_resize_src_lut_size,
+                    src_dy, resize_magnification_scale.y, 1.0);
+            }
+            else if(mask_type < 1.5)
+            {
+                pixel_color = downsample_vertical_sinc_tiled(
+                    mask_slot_texture_large, src_tex_uv, mask_resize_src_lut_size,
+                    src_dy, resize_magnification_scale.y, 1.0);
+            }
+            else
+            {
+                pixel_color = downsample_vertical_sinc_tiled(
+                    mask_shadow_texture_large, src_tex_uv, mask_resize_src_lut_size,
+                    src_dy, resize_magnification_scale.y, 1.0);
+            }
+            //  The input LUT was linear RGB, and so is our output:
+            FragColor = vec4(pixel_color, 1.0);
+        }
+        else
+        {
+            discard;
+        }
+    #else
+        discard;
+        FragColor = vec4(1.0);
+    #endif
+}
+#else
+layout(set = 0, binding = 3) uniform sampler2D mask_grille_texture_small;
+layout(set = 0, binding = 4) uniform sampler2D mask_slot_texture_small;
+layout(set = 0, binding = 5) uniform sampler2D mask_shadow_texture_small;
+
+void main()
+{
+    //  Resize the input phosphor mask tile to the final vertical size it will
+    //  appear on screen.  Keep 1x horizontal size if possible (IN.output_size
+    //  >= mask_resize_src_lut_size), and otherwise linearly sample horizontally
+    //  to fit exactly one tile.  Lanczos-resizing the phosphor mask achieves
+    //  much sharper results than mipmapping, and vertically resizing first
+    //  minimizes the total number of taps required.  We output a number of
+    //  resized tiles >= mask_resize_num_tiles for easier tiled sampling later.
+    #ifdef PHOSPHOR_MASK_MANUALLY_RESIZE
+        //  Discard unneeded fragments in case our profile allows real branches.
+        const vec2 tile_uv_wrap = src_tex_uv_wrap;
+        if(get_mask_sample_mode() < 0.5 &&
+            tile_uv_wrap.y <= mask_resize_num_tiles)
+        {
+            const float src_dy = 1.0/mask_resize_src_lut_size.y;
+            const vec2 src_tex_uv = frac(src_tex_uv_wrap);
+            vec3 pixel_color;
+            //  If mask_type is static, this branch will be resolved statically.
+            if(mask_type < 0.5)
+            {
+                pixel_color = downsample_vertical_sinc_tiled(
+                    mask_grille_texture_small, src_tex_uv, mask_resize_src_lut_size,
+                    src_dy, resize_magnification_scale.y, 1.0);
+            }
+            else if(mask_type < 1.5)
+            {
+                pixel_color = downsample_vertical_sinc_tiled(
+                    mask_slot_texture_small, src_tex_uv, mask_resize_src_lut_size,
+                    src_dy, resize_magnification_scale.y, 1.0);
+            }
+            else
+            {
+                pixel_color = downsample_vertical_sinc_tiled(
+                    mask_shadow_texture_small, src_tex_uv, mask_resize_src_lut_size,
+                    src_dy, resize_magnification_scale.y, 1.0);
+            }
+            //  The input LUT was linear RGB, and so is our output:
+            FragColor = vec4(pixel_color, 1.0);
+        }
+        else
+        {
+            discard;
+        }
+    #else
+        discard;
+        FragColor = vec4(1.0);
+    #endif
+}
+#endif

crt/shaders/crt-royale/src/gamma-management.h

@@ -1,3 +1,6 @@
+#ifndef GAMMA_MANAGEMENT_H
+#define GAMMA_MANAGEMENT_H
+
 ///////////////////////////////  BASE CONSTANTS  ///////////////////////////////
 
 //  Set standard gamma constants, but allow users to override them:
@@ -157,4 +160,6 @@ vec4 encode_output(const vec4 color)
 
 //#define tex2D_linearize(C, D, E) decode_input(vec4(texture(C, D, E)))
 //vec4 tex2D_linearize(const sampler2D tex, const vec2 tex_coords, const int texel_off)
-//{   return decode_input(vec4(texture(tex, tex_coords, texel_off)));    }
+//{   return decode_input(vec4(texture(tex, tex_coords, texel_off)));    }
+
+#endif  //  GAMMA_MANAGEMENT_H

crt/shaders/crt-royale/src/includes.h

@@ -7,4 +7,5 @@
 #include "gamma-management.h" //#include "../../../../include/gamma-management.h" <-move includes into crt-royale's src directory until it's actually working
 #include "blur-functions.h" //#include "../../../../include/blur-functions.h" <-move includes into crt-royale's src directory until it's actually working
 #include "scanline-functions.h"
-#include "bloom-functions.h"
+#include "bloom-functions.h"
+#include "phosphor-mask-resizing.h"

crt/shaders/crt-royale/src/phosphor-mask-resizing-old.h

@@ -0,0 +1,678 @@
+#ifndef PHOSPHOR_MASK_RESIZING_H
+#define PHOSPHOR_MASK_RESIZING_H
+
+/////////////////////////////  GPL LICENSE NOTICE  /////////////////////////////
+
+//  crt-royale: A full-featured CRT shader, with cheese.
+//  Copyright (C) 2014 TroggleMonkey <trogglemonkey@gmx.com>
+//
+//  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 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., 59 Temple
+//  Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+//////////////////////////////////  INCLUDES  //////////////////////////////////
+
+//#include "../user-settings.h"
+//#include "derived-settings-and-constants.h"
+#include "includes.h"
+
+/////////////////////////////  CODEPATH SELECTION  /////////////////////////////
+
+//  Choose a looping strategy based on what's allowed:
+//  Dynamic loops not allowed: Use a flat static loop.
+//  Dynamic loops accomodated: Coarsely branch around static loops.
+//  Dynamic loops assumed allowed: Use a flat dynamic loop.
+#ifndef DRIVERS_ALLOW_DYNAMIC_BRANCHES
+    #ifdef ACCOMODATE_POSSIBLE_DYNAMIC_LOOPS
+        #define BREAK_LOOPS_INTO_PIECES
+    #else
+        #define USE_SINGLE_STATIC_LOOP
+    #endif
+#endif  //  No else needed: Dynamic loops assumed.
+
+
+//////////////////////////////////  CONSTANTS  /////////////////////////////////
+
+//  The larger the resized tile, the fewer samples we'll need for downsizing.
+//  See if we can get a static min tile size > mask_min_allowed_tile_size:
+const float mask_min_allowed_tile_size = ceil(
+    mask_min_allowed_triad_size * mask_triads_per_tile);
+const float mask_min_expected_tile_size = 
+        mask_min_allowed_tile_size;
+//  Limit the number of sinc resize taps by the maximum minification factor:
+const float pi_over_lobes = pi/mask_sinc_lobes;
+const float max_sinc_resize_samples_float = 2.0 * mask_sinc_lobes *
+    mask_resize_src_lut_size.x/mask_min_expected_tile_size;
+//  Vectorized loops sample in multiples of 4.  Round up to be safe:
+const float max_sinc_resize_samples_m4 = ceil(
+    max_sinc_resize_samples_float * 0.25) * 4.0;
+
+
+/////////////////////////  RESAMPLING FUNCTION HELPERS  ////////////////////////
+
+float get_dynamic_loop_size(const float magnification_scale)
+{
+    //  Requires:   The following global constants must be defined:
+    //              1.) mask_sinc_lobes
+    //              2.) max_sinc_resize_samples_m4
+    //  Returns:    The minimum number of texture samples for a correct downsize
+    //              at magnification_scale.
+    //  We're downsizing, so the filter is sized across 2*lobes output pixels
+    //  (not 2*lobes input texels).  This impacts distance measurements and the
+    //  minimum number of input samples needed.
+    const float min_samples_float = 2.0 * mask_sinc_lobes / magnification_scale;
+    const float min_samples_m4 = ceil(min_samples_float * 0.25) * 4.0;
+    #ifdef DRIVERS_ALLOW_DYNAMIC_BRANCHES
+        const float max_samples_m4 = max_sinc_resize_samples_m4;
+    #else   // ifdef BREAK_LOOPS_INTO_PIECES
+        //  Simulating loops with branches imposes a 128-sample limit.
+        const float max_samples_m4 = min(128.0, max_sinc_resize_samples_m4);
+    #endif
+    return min(min_samples_m4, max_samples_m4);
+}
+
+vec2 get_first_texel_tile_uv_and_dist(const vec2 tex_uv, 
+    const vec2 texture_size, const float dr, 
+    const float input_tiles_per_texture_r, const float samples,
+    const bool vertical)
+{
+    //  Requires:   1.) dr == du == 1.0/texture_size.x or
+    //                  dr == dv == 1.0/texture_size.y
+    //                  (whichever direction we're resampling in).
+    //                  It's a scalar to save register space.
+    //              2.) input_tiles_per_texture_r is the number of input tiles
+    //                  that can fit in the input texture in the direction we're
+    //                  resampling this pass.
+    //              3.) vertical indicates whether we're resampling vertically
+    //                  this pass (or horizontally).
+    //  Returns:    Pack and return the first sample's tile_uv coord in [0, 1]
+    //              and its texel distance from the destination pixel, in the
+    //              resized dimension only.
+    //  We'll start with the topmost or leftmost sample and work down or right,
+    //  so get the first sample location and distance.  Modify both dimensions
+    //  as if we're doing a one-pass 2D resize; we'll throw away the unneeded
+    //  (and incorrect) dimension at the end.
+    const vec2 curr_texel = tex_uv * texture_size;
+    const vec2 prev_texel =
+        floor(curr_texel - vec2(under_half)) + vec2(0.5);
+    const vec2 first_texel = prev_texel - vec2(samples/2.0 - 1.0);
+    const vec2 first_texel_uv_wrap_2D = first_texel * dr;
+    const vec2 first_texel_dist_2D = curr_texel - first_texel;
+    //  Convert from tex_uv to tile_uv coords so we can sub fracts for fmods.
+    const vec2 first_texel_tile_uv_wrap_2D =
+        first_texel_uv_wrap_2D * input_tiles_per_texture_r;
+    //  Project wrapped coordinates to the [0, 1] range.  We'll do this with all
+    //  samples,but the first texel is special, since it might be negative.
+    const vec2 coord_negative =
+        vec2(first_texel_tile_uv_wrap_2D < vec2(0.0));
+    const vec2 first_texel_tile_uv_2D =
+        fract(first_texel_tile_uv_wrap_2D) + coord_negative;
+    //  Pack the first texel's tile_uv coord and texel distance in 1D:
+    const vec2 tile_u_and_dist =
+        vec2(first_texel_tile_uv_2D.x, first_texel_dist_2D.x);
+    const vec2 tile_v_and_dist =
+        vec2(first_texel_tile_uv_2D.y, first_texel_dist_2D.y);
+    return vertical ? tile_v_and_dist : tile_u_and_dist;
+    //return mix(tile_u_and_dist, tile_v_and_dist, float(vertical));
+}
+
+vec4 tex2Dlod0try(const sampler2D tex, const vec2 tex_uv)
+{
+    //  Mipmapping and anisotropic filtering get confused by sinc-resampling.
+    //  One [slow] workaround is to select the lowest mip level:
+    #ifdef ANISOTROPIC_RESAMPLING_COMPAT_TEX2DLOD
+        return tex2Dlod(tex, vec4(tex_uv, 0.0, 0.0));
+    #else
+        #ifdef ANISOTROPIC_RESAMPLING_COMPAT_TEX2DBIAS
+            return tex2Dbias(tex, vec4(tex_uv, 0.0, -16.0));
+        #else
+            return texture(tex, tex_uv);
+        #endif
+    #endif
+}
+
+
+//////////////////////////////  LOOP BODY MACROS  //////////////////////////////
+
+//  Using inline functions can exceed the temporary register limit, so we're
+//  stuck with #define macros (I'm TRULY sorry).  They're declared here instead
+//  of above to be closer to the actual invocation sites.  Steps:
+//  1.) Get the exact texel location.
+//  2.) Sample the phosphor mask (already assumed encoded in linear RGB).
+//  3.) Get the distance from the current pixel and sinc weight:
+//          sinc(dist) = sin(pi * dist)/(pi * dist)
+//      We can also use the slower/smoother Lanczos instead:
+//          L(x) = sinc(dist) * sinc(dist / lobes)
+//  4.) Accumulate the weight sum in weights, and accumulate the weighted texels
+//      in pixel_color (we'll normalize outside the loop at the end).
+//  We vectorize the loop to help reduce the Lanczos window's cost.
+
+    //  The r coord is the coord in the dimension we're resizing along (u or v),
+    //  and first_texel_tile_uv_rrrr is a vec4 of the first texel's u or v
+    //  tile_uv coord in [0, 1].  tex_uv_r will contain the tile_uv u or v coord
+    //  for four new texel samples.
+    #define CALCULATE_R_COORD_FOR_4_SAMPLES                                    \
+        const vec4 true_i = vec4(i_base + i) + vec4(0.0, 1.0, 2.0, 3.0); \
+        const vec4 tile_uv_r = fract(                                         \
+            first_texel_tile_uv_rrrr + true_i * tile_dr);                      \
+        const vec4 tex_uv_r = tile_uv_r * tile_size_uv_r;
+
+    #ifdef PHOSPHOR_MASK_RESIZE_LANCZOS_WINDOW
+        #define CALCULATE_SINC_RESAMPLE_WEIGHTS                                \
+            const vec4 pi_dist_over_lobes = pi_over_lobes * dist;            \
+            const vec4 weights = min(sin(pi_dist) * sin(pi_dist_over_lobes) /\
+                (pi_dist*pi_dist_over_lobes), vec4(1.0));
+    #else
+        #define CALCULATE_SINC_RESAMPLE_WEIGHTS                                \
+            const vec4 weights = min(sin(pi_dist)/pi_dist, vec4(1.0));
+    #endif
+
+    #define UPDATE_COLOR_AND_WEIGHT_SUMS                                       \
+        const vec4 dist = magnification_scale *                              \
+            abs(first_dist_unscaled - true_i);                                 \
+        const vec4 pi_dist = pi * dist;                                      \
+        CALCULATE_SINC_RESAMPLE_WEIGHTS;                                       \
+        pixel_color += new_sample0 * weights.xxx;                              \
+        pixel_color += new_sample1 * weights.yyy;                              \
+        pixel_color += new_sample2 * weights.zzz;                              \
+        pixel_color += new_sample3 * weights.www;                              \
+        weight_sum += weights;
+
+    #define VERTICAL_SINC_RESAMPLE_LOOP_BODY                                   \
+        CALCULATE_R_COORD_FOR_4_SAMPLES;                                       \
+        const vec3 new_sample0 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv.x, tex_uv_r.x)).rgb;                                 \
+        const vec3 new_sample1 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv.x, tex_uv_r.y)).rgb;                                 \
+        const vec3 new_sample2 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv.x, tex_uv_r.z)).rgb;                                 \
+        const vec3 new_sample3 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv.x, tex_uv_r.w)).rgb;                                 \
+        UPDATE_COLOR_AND_WEIGHT_SUMS;
+
+    #define HORIZONTAL_SINC_RESAMPLE_LOOP_BODY                                 \
+        CALCULATE_R_COORD_FOR_4_SAMPLES;                                       \
+        const vec3 new_sample0 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv_r.x, tex_uv.y)).rgb;                                 \
+        const vec3 new_sample1 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv_r.y, tex_uv.y)).rgb;                                 \
+        const vec3 new_sample2 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv_r.z, tex_uv.y)).rgb;                                 \
+        const vec3 new_sample3 = tex2Dlod0try(texture,                       \
+            vec2(tex_uv_r.w, tex_uv.y)).rgb;                                 \
+        UPDATE_COLOR_AND_WEIGHT_SUMS;
+
+
+////////////////////////////  RESAMPLING FUNCTIONS  ////////////////////////////
+
+vec3 downsample_vertical_sinc_tiled(const sampler2D texture,
+    const vec2 tex_uv, const vec2 texture_size, const float dr,
+    const float magnification_scale, const float tile_size_uv_r)
+{
+    //  Requires:   1.) dr == du == 1.0/texture_size.x or
+    //                  dr == dv == 1.0/texture_size.y
+    //                  (whichever direction we're resampling in).
+    //                  It's a scalar to save register space.
+    //              2.) tile_size_uv_r is the number of texels an input tile
+    //                  takes up in the input texture, in the direction we're
+    //                  resampling this pass.
+    //              3.) magnification_scale must be <= 1.0.
+    //  Returns:    Return a [Lanczos] sinc-resampled pixel of a vertically
+    //              downsized input tile embedded in an input texture.  (The
+    //              vertical version is special-cased though: It assumes the
+    //              tile size equals the [static] texture size, since it's used
+    //              on an LUT texture input containing one tile.  For more
+    //              generic use, eliminate the "static" in the parameters.)
+    //  The "r" in "dr," "tile_size_uv_r," etc. refers to the dimension
+    //  we're resizing along, e.g. "dy" in this case.
+    #ifdef USE_SINGLE_STATIC_LOOP
+        //  A static loop can be faster, but it might blur too much from using
+        //  more samples than it should.
+        const int samples = int(max_sinc_resize_samples_m4);
+    #else
+        const int samples = int(get_dynamic_loop_size(magnification_scale));
+    #endif
+
+    //  Get the first sample location (scalar tile uv coord along the resized
+    //  dimension) and distance from the output location (in texels):
+    const float input_tiles_per_texture_r = 1.0/tile_size_uv_r;
+    //  true = vertical resize:
+    const vec2 first_texel_tile_r_and_dist = get_first_texel_tile_uv_and_dist(
+        tex_uv, texture_size, dr, input_tiles_per_texture_r, samples, true);
+    const vec4 first_texel_tile_uv_rrrr = first_texel_tile_r_and_dist.xxxx;
+    const vec4 first_dist_unscaled = first_texel_tile_r_and_dist.yyyy;
+    //  Get the tile sample offset:
+    const float tile_dr = dr * input_tiles_per_texture_r;
+
+    //  Sum up each weight and weighted sample color, varying the looping
+    //  strategy based on our expected dynamic loop capabilities.  See the
+    //  loop body macros above.
+    int i_base = 0;
+    vec4 weight_sum = vec4(0.0);
+    vec3 pixel_color = vec3(0.0);
+    const int i_step = 4;
+    #ifdef BREAK_LOOPS_INTO_PIECES
+        if(samples - i_base >= 64)
+        {
+            for(int i = 0; i < 64; i += i_step)
+            {
+                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 64;
+        }
+        if(samples - i_base >= 32)
+        {
+            for(int i = 0; i < 32; i += i_step)
+            {
+                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 32;
+        }
+        if(samples - i_base >= 16)
+        {
+            for(int i = 0; i < 16; i += i_step)
+            {
+                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 16;
+        }
+        if(samples - i_base >= 8)
+        {
+            for(int i = 0; i < 8; i += i_step)
+            {
+                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 8;
+        }
+        if(samples - i_base >= 4)
+        {
+            for(int i = 0; i < 4; i += i_step)
+            {
+                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 4;
+        }
+        //  Do another 4-sample block for a total of 128 max samples.
+        if(samples - i_base > 0)
+        {
+            for(int i = 0; i < 4; i += i_step)
+            {
+                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+        }
+    #else
+        for(int i = 0; i < samples; i += i_step)
+        {
+            VERTICAL_SINC_RESAMPLE_LOOP_BODY;
+        }
+    #endif
+    //  Normalize so the weight_sum == 1.0, and return:
+    const vec2 weight_sum_reduce = weight_sum.xy + weight_sum.zw;
+    const vec3 scalar_weight_sum = vec3(weight_sum_reduce.x + 
+        weight_sum_reduce.y);
+    return (pixel_color/scalar_weight_sum);
+}
+
+vec3 downsample_horizontal_sinc_tiled(const sampler2D texture,
+    const vec2 tex_uv, const vec2 texture_size, const float dr,
+    const float magnification_scale, const float tile_size_uv_r)
+{
+    //  Differences from downsample_horizontal_sinc_tiled:
+    //  1.) The dr and tile_size_uv_r parameters are not static consts.
+    //  2.) The "vertical" parameter to get_first_texel_tile_uv_and_dist is
+    //      set to false instead of true.
+    //  3.) The horizontal version of the loop body is used.
+    //  TODO: If we can get guaranteed compile-time dead code elimination,
+    //  we can combine the vertical/horizontal downsampling functions by:
+    //  1.) Add an extra static const bool parameter called "vertical."
+    //  2.) Supply it with the result of get_first_texel_tile_uv_and_dist().
+    //  3.) Use a conditional assignment in the loop body macro.  This is the
+    //      tricky part: We DO NOT want to incur the extra conditional
+    //      assignment in the inner loop at runtime!
+    //  The "r" in "dr," "tile_size_uv_r," etc. refers to the dimension
+    //  we're resizing along, e.g. "dx" in this case.
+    #ifdef USE_SINGLE_STATIC_LOOP
+        //  If we have to load all samples, we might as well use them.
+        const int samples = int(max_sinc_resize_samples_m4);
+    #else
+        const int samples = int(get_dynamic_loop_size(magnification_scale));
+    #endif
+
+    //  Get the first sample location (scalar tile uv coord along resized
+    //  dimension) and distance from the output location (in texels):
+    const float input_tiles_per_texture_r = 1.0/tile_size_uv_r;
+    //  false = horizontal resize:
+    const vec2 first_texel_tile_r_and_dist = get_first_texel_tile_uv_and_dist(
+        tex_uv, texture_size, dr, input_tiles_per_texture_r, samples, false);
+    const vec4 first_texel_tile_uv_rrrr = first_texel_tile_r_and_dist.xxxx;
+    const vec4 first_dist_unscaled = first_texel_tile_r_and_dist.yyyy;
+    //  Get the tile sample offset:
+    const float tile_dr = dr * input_tiles_per_texture_r;
+
+    //  Sum up each weight and weighted sample color, varying the looping
+    //  strategy based on our expected dynamic loop capabilities.  See the
+    //  loop body macros above.
+    int i_base = 0;
+    vec4 weight_sum = vec4(0.0);
+    vec3 pixel_color = vec3(0.0);
+    const int i_step = 4;
+    #ifdef BREAK_LOOPS_INTO_PIECES
+        if(samples - i_base >= 64)
+        {
+            for(int i = 0; i < 64; i += i_step)
+            {
+                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 64;
+        }
+        if(samples - i_base >= 32)
+        {
+            for(int i = 0; i < 32; i += i_step)
+            {
+                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 32;
+        }
+        if(samples - i_base >= 16)
+        {
+            for(int i = 0; i < 16; i += i_step)
+            {
+                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 16;
+        }
+        if(samples - i_base >= 8)
+        {
+            for(int i = 0; i < 8; i += i_step)
+            {
+                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 8;
+        }
+        if(samples - i_base >= 4)
+        {
+            for(int i = 0; i < 4; i += i_step)
+            {
+                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+            i_base += 4;
+        }
+        //  Do another 4-sample block for a total of 128 max samples.
+        if(samples - i_base > 0)
+        {
+            for(int i = 0; i < 4; i += i_step)
+            {
+                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
+            }
+        }
+    #else
+        for(int i = 0; i < samples; i += i_step)
+        {
+            HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
+        }
+    #endif
+    //  Normalize so the weight_sum == 1.0, and return:
+    const vec2 weight_sum_reduce = weight_sum.xy + weight_sum.zw;
+    const vec3 scalar_weight_sum = vec3(weight_sum_reduce.x +
+        weight_sum_reduce.y);
+    return (pixel_color/scalar_weight_sum);
+}
+
+
+////////////////////////////  TILE SIZE CALCULATION  ///////////////////////////
+
+vec2 get_resized_mask_tile_size(const vec2 estimated_viewport_size,
+    const vec2 estimated_mask_resize_output_size,
+    const bool solemnly_swear_same_inputs_for_every_pass)
+{
+    //  Requires:   The following global constants must be defined according to
+    //              certain constraints:
+    //              1.) mask_resize_num_triads: Must be high enough that our
+    //                  mask sampling method won't have artifacts later
+    //                  (long story; see derived-settings-and-constants.h)
+    //              2.) mask_resize_src_lut_size: Texel size of our mask LUT
+    //              3.) mask_triads_per_tile: Num horizontal triads in our LUT
+    //              4.) mask_min_allowed_triad_size: User setting (the more
+    //                  restrictive it is, the faster the resize will go)
+    //              5.) mask_min_allowed_tile_size_x < mask_resize_src_lut_size.x
+    //              6.) mask_triad_size_desired_{runtime, static}
+    //              7.) mask_num_triads_desired_{runtime, static}
+    //              8.) mask_specify_num_triads must be 0.0/1.0 (false/true)
+    //              The function parameters must be defined as follows:
+    //              1.) estimated_viewport_size == (final viewport size);
+    //                  If mask_specify_num_triads is 1.0/true and the viewport
+    //                  estimate is wrong, the number of triads will differ from
+    //                  the user's preference by about the same factor.
+    //              2.) estimated_mask_resize_output_size: Must equal the
+    //                  output size of the MASK_RESIZE pass.
+    //                  Exception: The x component may be estimated garbage if
+    //                  and only if the caller throws away the x result.
+    //              3.) solemnly_swear_same_inputs_for_every_pass: Set to false,
+    //                  unless you can guarantee that every call across every
+    //                  pass will use the same sizes for the other parameters.
+    //              When calling this across multiple passes, always use the
+    //              same y viewport size/scale, and always use the same x
+    //              viewport size/scale when using the x result.
+    //  Returns:    Return the final size of a manually resized mask tile, after
+    //              constraining the desired size to avoid artifacts.  Under
+    //              unusual circumstances, tiles may become stretched vertically
+    //              (see wall of text below).
+    //  Stated tile properties must be correct:
+    const float tile_aspect_ratio_inv =
+        mask_resize_src_lut_size.y/mask_resize_src_lut_size.x;
+    const float tile_aspect_ratio = 1.0/tile_aspect_ratio_inv;
+    const vec2 tile_aspect = vec2(1.0, tile_aspect_ratio_inv);
+    //  If mask_specify_num_triads is 1.0/true and estimated_viewport_size.x is
+    //  wrong, the user preference will be misinterpreted:
+    const float desired_tile_size_x = mask_triads_per_tile * mix(
+        mask_triad_size_desired,
+        estimated_viewport_size.x / mask_num_triads_desired,
+        mask_specify_num_triads);
+    if(get_mask_sample_mode() > 0.5)
+    {
+        //  We don't need constraints unless we're sampling MASK_RESIZE.
+        return desired_tile_size_x * tile_aspect;
+    }
+    //  Make sure we're not upsizing:
+    const float temp_tile_size_x =
+        min(desired_tile_size_x, mask_resize_src_lut_size.x);
+    //  Enforce min_tile_size and max_tile_size in both dimensions:
+    const vec2 temp_tile_size = temp_tile_size_x * tile_aspect;
+    const vec2 min_tile_size =
+        mask_min_allowed_tile_size * tile_aspect;
+    const vec2 max_tile_size =
+        estimated_mask_resize_output_size / mask_resize_num_tiles;
+    const vec2 clamped_tile_size =
+        clamp(temp_tile_size, min_tile_size, max_tile_size);
+    //  Try to maintain tile_aspect_ratio.  This is the tricky part:
+    //  If we're currently resizing in the y dimension, the x components
+    //  could be MEANINGLESS.  (If estimated_mask_resize_output_size.x is
+    //  bogus, then so is max_tile_size.x and clamped_tile_size.x.)
+    //  We can't adjust the y size based on clamped_tile_size.x.  If it
+    //  clamps when it shouldn't, it won't clamp again when later passes
+    //  call this function with the correct sizes, and the discrepancy will
+    //  break the sampling coords in MASKED_SCANLINES.  Instead, we'll limit
+    //  the x size based on the y size, but not vice versa, unless the
+    //  caller swears the parameters were the same (correct) in every pass.
+    //  As a result, triads could appear vertically stretched if:
+    //  a.) mask_resize_src_lut_size.x > mask_resize_src_lut_size.y: Wide
+    //      LUT's might clamp x more than y (all provided LUT's are square)
+    //  b.) true_viewport_size.x < true_viewport_size.y: The user is playing
+    //      with a vertically oriented screen (not accounted for anyway)
+    //  c.) mask_resize_viewport_scale.x < masked_resize_viewport_scale.y:
+    //      Viewport scales are equal by default.
+    //  If any of these are the case, you can fix the stretching by setting:
+    //      mask_resize_viewport_scale.x = mask_resize_viewport_scale.y *
+    //          (1.0 / min_expected_aspect_ratio) *
+    //          (mask_resize_src_lut_size.x / mask_resize_src_lut_size.y)
+    const float x_tile_size_from_y =
+        clamped_tile_size.y * tile_aspect_ratio;
+    const float y_tile_size_from_x = mix(clamped_tile_size.y,
+        clamped_tile_size.x * tile_aspect_ratio_inv,
+        float(solemnly_swear_same_inputs_for_every_pass));
+    const vec2 reclamped_tile_size = vec2(
+        min(clamped_tile_size.x, x_tile_size_from_y),
+        min(clamped_tile_size.y, y_tile_size_from_x));
+    //  We need integer tile sizes in both directions for tiled sampling to
+    //  work correctly.  Use floor (to make sure we don't round up), but be
+    //  careful to avoid a rounding bug where floor decreases whole numbers:
+    const vec2 final_resized_tile_size =
+        floor(reclamped_tile_size + vec2(FIX_ZERO(0.0)));
+    return final_resized_tile_size;
+}
+
+
+/////////////////////////  FINAL MASK SAMPLING HELPERS  ////////////////////////
+
+vec4 get_mask_sampling_parameters(const vec2 mask_resize_texture_size,
+    const vec2 mask_resize_video_size, const vec2 true_viewport_size,
+    out vec2 mask_tiles_per_screen)
+{
+    //  Requires:   1.) Requirements of get_resized_mask_tile_size() must be
+    //                  met, particularly regarding global constants.
+    //              The function parameters must be defined as follows:
+    //              1.) mask_resize_texture_size == MASK_RESIZE.texture_size
+    //                  if get_mask_sample_mode() is 0 (otherwise anything)
+    //              2.) mask_resize_video_size == MASK_RESIZE.video_size
+    //                  if get_mask_sample_mode() is 0 (otherwise anything)
+    //              3.) true_viewport_size == IN.output_size for a pass set to
+    //                  1.0 viewport scale (i.e. it must be correct)
+    //  Returns:    Return a vec4 containing:
+    //                  xy: tex_uv coords for the start of the mask tile
+    //                  zw: tex_uv size of the mask tile from start to end
+    //              mask_tiles_per_screen is an out parameter containing the
+    //              number of mask tiles that will fit on the screen.
+    //  First get the final resized tile size.  The viewport size and mask
+    //  resize viewport scale must be correct, but don't solemnly swear they
+    //  were correct in both mask resize passes unless you know it's true.
+    //  (We can better ensure a correct tile aspect ratio if the parameters are
+    //  guaranteed correct in all passes...but if we lie, we'll get inconsistent
+    //  sizes across passes, resulting in broken texture coordinates.)
+    const float mask_sample_mode = get_mask_sample_mode();
+    const vec2 mask_resize_tile_size = get_resized_mask_tile_size(
+        true_viewport_size, mask_resize_video_size, false);
+    if(mask_sample_mode < 0.5)
+    {
+        //  Sample MASK_RESIZE: The resized tile is a fracttion of the texture
+        //  size and starts at a nonzero offset to allow for border texels:
+        const vec2 mask_tile_uv_size = mask_resize_tile_size /
+            mask_resize_texture_size;
+        const vec2 skipped_tiles = mask_start_texels/mask_resize_tile_size;
+        const vec2 mask_tile_start_uv = skipped_tiles * mask_tile_uv_size;
+        //  mask_tiles_per_screen must be based on the *true* viewport size:
+        mask_tiles_per_screen = true_viewport_size / mask_resize_tile_size;
+        return vec4(mask_tile_start_uv, mask_tile_uv_size);
+    }
+    else
+    {
+        //  If we're tiling at the original size (1:1 pixel:texel), redefine a
+        //  "tile" to be the full texture containing many triads.  Otherwise,
+        //  we're hardware-resampling an LUT, and the texture truly contains a
+        //  single unresized phosphor mask tile anyway.
+        const vec2 mask_tile_uv_size = vec2(1.0);
+        const vec2 mask_tile_start_uv = vec2(0.0);
+        if(mask_sample_mode > 1.5)
+        {
+            //  Repeat the full LUT at a 1:1 pixel:texel ratio without resizing:
+            mask_tiles_per_screen = true_viewport_size/mask_texture_large_size;
+        }
+        else
+        {
+            //  Hardware-resize the original LUT:
+            mask_tiles_per_screen = true_viewport_size / mask_resize_tile_size;
+        }
+        return vec4(mask_tile_start_uv, mask_tile_uv_size);
+    }
+}
+
+vec2 fix_tiling_discontinuities_normalized(const vec2 tile_uv,
+    vec2 duv_dx, vec2 duv_dy)
+{
+    //  Requires:   1.) duv_dx == ddx(tile_uv)
+    //              2.) duv_dy == ddy(tile_uv)
+    //              3.) tile_uv contains tile-relative uv coords in [0, 1],
+    //                  such that (0.5, 0.5) is the center of a tile, etc.
+    //                  ("Tile" can mean texture, the video embedded in the
+    //                  texture, or some other "tile" embedded in a texture.)
+    //  Returns:    Return new tile_uv coords that contain no discontinuities
+    //              across a 2x2 pixel quad.
+    //  Description:
+    //  When uv coords wrap from 1.0 to 0.0, they create a discontinuity in the
+    //  derivatives, which we assume happened if the absolute difference between
+    //  any fragment in a 2x2 block is > ~half a tile.  If the current block has
+    //  a u or v discontinuity and the current fragment is in the first half of
+    //  the tile along that axis (i.e. it wrapped from 1.0 to 0.0), add a tile
+    //  to that coord to make the 2x2 block continuous.  (It will now have a
+    //  coord > 1.0 in the padding area beyond the tile.)  This function takes
+    //  derivatives as parameters so the caller can reuse them.
+    //  In case we're using high-quality (nVidia-style) derivatives, ensure
+    //  diagonically opposite fragments see each other for correctness:
+    duv_dx = abs(duv_dx) + abs(ddy(duv_dx));
+    duv_dy = abs(duv_dy) + abs(ddx(duv_dy));
+    const vec2 pixel_in_first_half_tile = vec2(tile_uv < vec2(0.5));
+    const vec2 jump_exists = vec2(duv_dx + duv_dy > vec2(0.5));
+    return tile_uv + jump_exists * pixel_in_first_half_tile;
+}
+
+vec2 convert_phosphor_tile_uv_wrap_to_tex_uv(const vec2 tile_uv_wrap,
+    const vec4 mask_tile_start_uv_and_size)
+{
+    //  Requires:   1.) tile_uv_wrap contains tile-relative uv coords, where the
+    //                  tile spans from [0, 1], such that (0.5, 0.5) is at the
+    //                  tile center.  The input coords can range from [0, inf],
+    //                  and their fracttional parts map to a repeated tile.
+    //                  ("Tile" can mean texture, the video embedded in the
+    //                  texture, or some other "tile" embedded in a texture.)
+    //              2.) mask_tile_start_uv_and_size.xy contains tex_uv coords
+    //                  for the start of the embedded tile in the full texture.
+    //              3.) mask_tile_start_uv_and_size.zw contains the [fracttional]
+    //                  tex_uv size of the embedded tile in the full texture.
+    //  Returns:    Return tex_uv coords (used for texture sampling)
+    //              corresponding to tile_uv_wrap.
+    if(get_mask_sample_mode() < 0.5)
+    {
+        //  Manually repeat the resized mask tile to fill the screen:
+        //  First get fracttional tile_uv coords.  Using fract/fmod on coords
+        //  confuses anisotropic filtering; fix it as user options dictate.
+        //  derived-settings-and-constants.h disables incompatible options.
+        #ifdef ANISOTROPIC_TILING_COMPAT_TILE_FLAT_TWICE
+            vec2 tile_uv = fract(tile_uv_wrap * 0.5) * 2.0;
+        #else
+            vec2 tile_uv = fract(tile_uv_wrap);
+        #endif
+        #ifdef ANISOTROPIC_TILING_COMPAT_FIX_DISCONTINUITIES
+            const vec2 tile_uv_dx = ddx(tile_uv);
+            const vec2 tile_uv_dy = ddy(tile_uv);
+            tile_uv = fix_tiling_discontinuities_normalized(tile_uv,
+                tile_uv_dx, tile_uv_dy);
+        #endif
+        //  The tile is embedded in a padded FBO, and it may start at a
+        //  nonzero offset if border texels are used to avoid artifacts:
+        const vec2 mask_tex_uv = mask_tile_start_uv_and_size.xy +
+            tile_uv * mask_tile_start_uv_and_size.zw;
+        return mask_tex_uv;
+    }
+    else
+    {
+        //  Sample from the input phosphor mask texture with hardware tiling.
+        //  If we're tiling at the original size (mode 2), the "tile" is the
+        //  whole texture, and it contains a large number of triads mapped with
+        //  a 1:1 pixel:texel ratio.  OTHERWISE, the texture contains a single
+        //  unresized tile.  tile_uv_wrap already has correct coords for both!
+        return tile_uv_wrap;
+    }
+}
+
+
+#endif  //  PHOSPHOR_MASK_RESIZING_H
+

crt/shaders/crt-royale/src/phosphor-mask-resizing.h

@@ -22,8 +22,9 @@
 
 //////////////////////////////////  INCLUDES  //////////////////////////////////
 
-#include "../user-settings.h"
-#include "derived-settings-and-constants.h"
+//#include "../user-settings.h"
+//#include "derived-settings-and-constants.h"
+#include "includes.h"
 
 /////////////////////////////  CODEPATH SELECTION  /////////////////////////////
 
@@ -59,7 +60,7 @@ const float max_sinc_resize_samples_m4 = ceil(
 
 /////////////////////////  RESAMPLING FUNCTION HELPERS  ////////////////////////
 
-inline float get_dynamic_loop_size(const float magnification_scale)
+float get_dynamic_loop_size(const float magnification_scale)
 {
     //  Requires:   The following global constants must be defined:
     //              1.) mask_sinc_lobes
@@ -107,7 +108,7 @@ vec2 get_first_texel_tile_uv_and_dist(const vec2 tex_uv,
     const vec2 first_texel = prev_texel - vec2(samples/2.0 - 1.0);
     const vec2 first_texel_uv_wrap_2D = first_texel * dr;
     const vec2 first_texel_dist_2D = curr_texel - first_texel;
-    //  Convert from tex_uv to tile_uv coords so we can sub fracs for fmods.
+    //  Convert from tex_uv to tile_uv coords so we can sub fracts for fmods.
     const vec2 first_texel_tile_uv_wrap_2D =
         first_texel_uv_wrap_2D * input_tiles_per_texture_r;
     //  Project wrapped coordinates to the [0, 1] range.  We'll do this with all
@@ -115,17 +116,17 @@ vec2 get_first_texel_tile_uv_and_dist(const vec2 tex_uv,
     const vec2 coord_negative =
         vec2(first_texel_tile_uv_wrap_2D < vec2(0.0));
     const vec2 first_texel_tile_uv_2D =
-        frac(first_texel_tile_uv_wrap_2D) + coord_negative;
+        fract(first_texel_tile_uv_wrap_2D) + coord_negative;
     //  Pack the first texel's tile_uv coord and texel distance in 1D:
     const vec2 tile_u_and_dist =
         vec2(first_texel_tile_uv_2D.x, first_texel_dist_2D.x);
     const vec2 tile_v_and_dist =
         vec2(first_texel_tile_uv_2D.y, first_texel_dist_2D.y);
     return vertical ? tile_v_and_dist : tile_u_and_dist;
-    //return lerp(tile_u_and_dist, tile_v_and_dist, float(vertical));
+    //return mix(tile_u_and_dist, tile_v_and_dist, float(vertical));
 }
 
-inline vec4 tex2Dlod0try(const sampler2D tex, const vec2 tex_uv)
+vec4 tex2Dlod0try(const sampler2D tex, const vec2 tex_uv)
 {
     //  Mipmapping and anisotropic filtering get confused by sinc-resampling.
     //  One [slow] workaround is to select the lowest mip level:
@@ -162,7 +163,7 @@ inline vec4 tex2Dlod0try(const sampler2D tex, const vec2 tex_uv)
     //  for four new texel samples.
     #define CALCULATE_R_COORD_FOR_4_SAMPLES                                    \
         const vec4 true_i = vec4(i_base + i) + vec4(0.0, 1.0, 2.0, 3.0); \
-        const vec4 tile_uv_r = frac(                                         \
+        const vec4 tile_uv_r = fract(                                         \
             first_texel_tile_uv_rrrr + true_i * tile_dr);                      \
         const vec4 tex_uv_r = tile_uv_r * tile_size_uv_r;
 
@@ -214,220 +215,6 @@ inline vec4 tex2Dlod0try(const sampler2D tex, const vec2 tex_uv)
 
 ////////////////////////////  RESAMPLING FUNCTIONS  ////////////////////////////
 
-vec3 downsample_vertical_sinc_tiled(const sampler2D texture,
-    const vec2 tex_uv, const vec2 texture_size, const float dr,
-    const float magnification_scale, const float tile_size_uv_r)
-{
-    //  Requires:   1.) dr == du == 1.0/texture_size.x or
-    //                  dr == dv == 1.0/texture_size.y
-    //                  (whichever direction we're resampling in).
-    //                  It's a scalar to save register space.
-    //              2.) tile_size_uv_r is the number of texels an input tile
-    //                  takes up in the input texture, in the direction we're
-    //                  resampling this pass.
-    //              3.) magnification_scale must be <= 1.0.
-    //  Returns:    Return a [Lanczos] sinc-resampled pixel of a vertically
-    //              downsized input tile embedded in an input texture.  (The
-    //              vertical version is special-cased though: It assumes the
-    //              tile size equals the [static] texture size, since it's used
-    //              on an LUT texture input containing one tile.  For more
-    //              generic use, eliminate the "static" in the parameters.)
-    //  The "r" in "dr," "tile_size_uv_r," etc. refers to the dimension
-    //  we're resizing along, e.g. "dy" in this case.
-    #ifdef USE_SINGLE_STATIC_LOOP
-        //  A static loop can be faster, but it might blur too much from using
-        //  more samples than it should.
-        const int samples = int(max_sinc_resize_samples_m4);
-    #else
-        const int samples = int(get_dynamic_loop_size(magnification_scale));
-    #endif
-
-    //  Get the first sample location (scalar tile uv coord along the resized
-    //  dimension) and distance from the output location (in texels):
-    const float input_tiles_per_texture_r = 1.0/tile_size_uv_r;
-    //  true = vertical resize:
-    const vec2 first_texel_tile_r_and_dist = get_first_texel_tile_uv_and_dist(
-        tex_uv, texture_size, dr, input_tiles_per_texture_r, samples, true);
-    const vec4 first_texel_tile_uv_rrrr = first_texel_tile_r_and_dist.xxxx;
-    const vec4 first_dist_unscaled = first_texel_tile_r_and_dist.yyyy;
-    //  Get the tile sample offset:
-    const float tile_dr = dr * input_tiles_per_texture_r;
-
-    //  Sum up each weight and weighted sample color, varying the looping
-    //  strategy based on our expected dynamic loop capabilities.  See the
-    //  loop body macros above.
-    int i_base = 0;
-    vec4 weight_sum = vec4(0.0);
-    vec3 pixel_color = vec3(0.0);
-    const int i_step = 4;
-    #ifdef BREAK_LOOPS_INTO_PIECES
-        if(samples - i_base >= 64)
-        {
-            for(int i = 0; i < 64; i += i_step)
-            {
-                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 64;
-        }
-        if(samples - i_base >= 32)
-        {
-            for(int i = 0; i < 32; i += i_step)
-            {
-                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 32;
-        }
-        if(samples - i_base >= 16)
-        {
-            for(int i = 0; i < 16; i += i_step)
-            {
-                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 16;
-        }
-        if(samples - i_base >= 8)
-        {
-            for(int i = 0; i < 8; i += i_step)
-            {
-                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 8;
-        }
-        if(samples - i_base >= 4)
-        {
-            for(int i = 0; i < 4; i += i_step)
-            {
-                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 4;
-        }
-        //  Do another 4-sample block for a total of 128 max samples.
-        if(samples - i_base > 0)
-        {
-            for(int i = 0; i < 4; i += i_step)
-            {
-                VERTICAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-        }
-    #else
-        for(int i = 0; i < samples; i += i_step)
-        {
-            VERTICAL_SINC_RESAMPLE_LOOP_BODY;
-        }
-    #endif
-    //  Normalize so the weight_sum == 1.0, and return:
-    const vec2 weight_sum_reduce = weight_sum.xy + weight_sum.zw;
-    const vec3 scalar_weight_sum = vec3(weight_sum_reduce.x + 
-        weight_sum_reduce.y);
-    return (pixel_color/scalar_weight_sum);
-}
-
-vec3 downsample_horizontal_sinc_tiled(const sampler2D texture,
-    const vec2 tex_uv, const vec2 texture_size, const float dr,
-    const float magnification_scale, const float tile_size_uv_r)
-{
-    //  Differences from downsample_horizontal_sinc_tiled:
-    //  1.) The dr and tile_size_uv_r parameters are not static consts.
-    //  2.) The "vertical" parameter to get_first_texel_tile_uv_and_dist is
-    //      set to false instead of true.
-    //  3.) The horizontal version of the loop body is used.
-    //  TODO: If we can get guaranteed compile-time dead code elimination,
-    //  we can combine the vertical/horizontal downsampling functions by:
-    //  1.) Add an extra static const bool parameter called "vertical."
-    //  2.) Supply it with the result of get_first_texel_tile_uv_and_dist().
-    //  3.) Use a conditional assignment in the loop body macro.  This is the
-    //      tricky part: We DO NOT want to incur the extra conditional
-    //      assignment in the inner loop at runtime!
-    //  The "r" in "dr," "tile_size_uv_r," etc. refers to the dimension
-    //  we're resizing along, e.g. "dx" in this case.
-    #ifdef USE_SINGLE_STATIC_LOOP
-        //  If we have to load all samples, we might as well use them.
-        const int samples = int(max_sinc_resize_samples_m4);
-    #else
-        const int samples = int(get_dynamic_loop_size(magnification_scale));
-    #endif
-
-    //  Get the first sample location (scalar tile uv coord along resized
-    //  dimension) and distance from the output location (in texels):
-    const float input_tiles_per_texture_r = 1.0/tile_size_uv_r;
-    //  false = horizontal resize:
-    const vec2 first_texel_tile_r_and_dist = get_first_texel_tile_uv_and_dist(
-        tex_uv, texture_size, dr, input_tiles_per_texture_r, samples, false);
-    const vec4 first_texel_tile_uv_rrrr = first_texel_tile_r_and_dist.xxxx;
-    const vec4 first_dist_unscaled = first_texel_tile_r_and_dist.yyyy;
-    //  Get the tile sample offset:
-    const float tile_dr = dr * input_tiles_per_texture_r;
-
-    //  Sum up each weight and weighted sample color, varying the looping
-    //  strategy based on our expected dynamic loop capabilities.  See the
-    //  loop body macros above.
-    int i_base = 0;
-    vec4 weight_sum = vec4(0.0);
-    vec3 pixel_color = vec3(0.0);
-    const int i_step = 4;
-    #ifdef BREAK_LOOPS_INTO_PIECES
-        if(samples - i_base >= 64)
-        {
-            for(int i = 0; i < 64; i += i_step)
-            {
-                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 64;
-        }
-        if(samples - i_base >= 32)
-        {
-            for(int i = 0; i < 32; i += i_step)
-            {
-                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 32;
-        }
-        if(samples - i_base >= 16)
-        {
-            for(int i = 0; i < 16; i += i_step)
-            {
-                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 16;
-        }
-        if(samples - i_base >= 8)
-        {
-            for(int i = 0; i < 8; i += i_step)
-            {
-                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 8;
-        }
-        if(samples - i_base >= 4)
-        {
-            for(int i = 0; i < 4; i += i_step)
-            {
-                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-            i_base += 4;
-        }
-        //  Do another 4-sample block for a total of 128 max samples.
-        if(samples - i_base > 0)
-        {
-            for(int i = 0; i < 4; i += i_step)
-            {
-                HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
-            }
-        }
-    #else
-        for(int i = 0; i < samples; i += i_step)
-        {
-            HORIZONTAL_SINC_RESAMPLE_LOOP_BODY;
-        }
-    #endif
-    //  Normalize so the weight_sum == 1.0, and return:
-    const vec2 weight_sum_reduce = weight_sum.xy + weight_sum.zw;
-    const vec3 scalar_weight_sum = vec3(weight_sum_reduce.x +
-        weight_sum_reduce.y);
-    return (pixel_color/scalar_weight_sum);
-}
-
-
 ////////////////////////////  TILE SIZE CALCULATION  ///////////////////////////
 
 vec2 get_resized_mask_tile_size(const vec2 estimated_viewport_size,
@@ -473,7 +260,7 @@ vec2 get_resized_mask_tile_size(const vec2 estimated_viewport_size,
     const vec2 tile_aspect = vec2(1.0, tile_aspect_ratio_inv);
     //  If mask_specify_num_triads is 1.0/true and estimated_viewport_size.x is
     //  wrong, the user preference will be misinterpreted:
-    const float desired_tile_size_x = mask_triads_per_tile * lerp(
+    const float desired_tile_size_x = mask_triads_per_tile * mix(
         mask_triad_size_desired,
         estimated_viewport_size.x / mask_num_triads_desired,
         mask_specify_num_triads);
@@ -516,7 +303,7 @@ vec2 get_resized_mask_tile_size(const vec2 estimated_viewport_size,
     //          (mask_resize_src_lut_size.x / mask_resize_src_lut_size.y)
     const float x_tile_size_from_y =
         clamped_tile_size.y * tile_aspect_ratio;
-    const float y_tile_size_from_x = lerp(clamped_tile_size.y,
+    const float y_tile_size_from_x = mix(clamped_tile_size.y,
         clamped_tile_size.x * tile_aspect_ratio_inv,
         float(solemnly_swear_same_inputs_for_every_pass));
     const vec2 reclamped_tile_size = vec2(
@@ -531,147 +318,6 @@ vec2 get_resized_mask_tile_size(const vec2 estimated_viewport_size,
 }
 
 
-/////////////////////////  FINAL MASK SAMPLING HELPERS  ////////////////////////
-
-vec4 get_mask_sampling_parameters(const vec2 mask_resize_texture_size,
-    const vec2 mask_resize_video_size, const vec2 true_viewport_size,
-    out vec2 mask_tiles_per_screen)
-{
-    //  Requires:   1.) Requirements of get_resized_mask_tile_size() must be
-    //                  met, particularly regarding global constants.
-    //              The function parameters must be defined as follows:
-    //              1.) mask_resize_texture_size == MASK_RESIZE.texture_size
-    //                  if get_mask_sample_mode() is 0 (otherwise anything)
-    //              2.) mask_resize_video_size == MASK_RESIZE.video_size
-    //                  if get_mask_sample_mode() is 0 (otherwise anything)
-    //              3.) true_viewport_size == IN.output_size for a pass set to
-    //                  1.0 viewport scale (i.e. it must be correct)
-    //  Returns:    Return a vec4 containing:
-    //                  xy: tex_uv coords for the start of the mask tile
-    //                  zw: tex_uv size of the mask tile from start to end
-    //              mask_tiles_per_screen is an out parameter containing the
-    //              number of mask tiles that will fit on the screen.
-    //  First get the final resized tile size.  The viewport size and mask
-    //  resize viewport scale must be correct, but don't solemnly swear they
-    //  were correct in both mask resize passes unless you know it's true.
-    //  (We can better ensure a correct tile aspect ratio if the parameters are
-    //  guaranteed correct in all passes...but if we lie, we'll get inconsistent
-    //  sizes across passes, resulting in broken texture coordinates.)
-    const float mask_sample_mode = get_mask_sample_mode();
-    const vec2 mask_resize_tile_size = get_resized_mask_tile_size(
-        true_viewport_size, mask_resize_video_size, false);
-    if(mask_sample_mode < 0.5)
-    {
-        //  Sample MASK_RESIZE: The resized tile is a fraction of the texture
-        //  size and starts at a nonzero offset to allow for border texels:
-        const vec2 mask_tile_uv_size = mask_resize_tile_size /
-            mask_resize_texture_size;
-        const vec2 skipped_tiles = mask_start_texels/mask_resize_tile_size;
-        const vec2 mask_tile_start_uv = skipped_tiles * mask_tile_uv_size;
-        //  mask_tiles_per_screen must be based on the *true* viewport size:
-        mask_tiles_per_screen = true_viewport_size / mask_resize_tile_size;
-        return vec4(mask_tile_start_uv, mask_tile_uv_size);
-    }
-    else
-    {
-        //  If we're tiling at the original size (1:1 pixel:texel), redefine a
-        //  "tile" to be the full texture containing many triads.  Otherwise,
-        //  we're hardware-resampling an LUT, and the texture truly contains a
-        //  single unresized phosphor mask tile anyway.
-        const vec2 mask_tile_uv_size = vec2(1.0);
-        const vec2 mask_tile_start_uv = vec2(0.0);
-        if(mask_sample_mode > 1.5)
-        {
-            //  Repeat the full LUT at a 1:1 pixel:texel ratio without resizing:
-            mask_tiles_per_screen = true_viewport_size/mask_texture_large_size;
-        }
-        else
-        {
-            //  Hardware-resize the original LUT:
-            mask_tiles_per_screen = true_viewport_size / mask_resize_tile_size;
-        }
-        return vec4(mask_tile_start_uv, mask_tile_uv_size);
-    }
-}
-
-vec2 fix_tiling_discontinuities_normalized(const vec2 tile_uv,
-    vec2 duv_dx, vec2 duv_dy)
-{
-    //  Requires:   1.) duv_dx == ddx(tile_uv)
-    //              2.) duv_dy == ddy(tile_uv)
-    //              3.) tile_uv contains tile-relative uv coords in [0, 1],
-    //                  such that (0.5, 0.5) is the center of a tile, etc.
-    //                  ("Tile" can mean texture, the video embedded in the
-    //                  texture, or some other "tile" embedded in a texture.)
-    //  Returns:    Return new tile_uv coords that contain no discontinuities
-    //              across a 2x2 pixel quad.
-    //  Description:
-    //  When uv coords wrap from 1.0 to 0.0, they create a discontinuity in the
-    //  derivatives, which we assume happened if the absolute difference between
-    //  any fragment in a 2x2 block is > ~half a tile.  If the current block has
-    //  a u or v discontinuity and the current fragment is in the first half of
-    //  the tile along that axis (i.e. it wrapped from 1.0 to 0.0), add a tile
-    //  to that coord to make the 2x2 block continuous.  (It will now have a
-    //  coord > 1.0 in the padding area beyond the tile.)  This function takes
-    //  derivatives as parameters so the caller can reuse them.
-    //  In case we're using high-quality (nVidia-style) derivatives, ensure
-    //  diagonically opposite fragments see each other for correctness:
-    duv_dx = abs(duv_dx) + abs(ddy(duv_dx));
-    duv_dy = abs(duv_dy) + abs(ddx(duv_dy));
-    const vec2 pixel_in_first_half_tile = vec2(tile_uv < vec2(0.5));
-    const vec2 jump_exists = vec2(duv_dx + duv_dy > vec2(0.5));
-    return tile_uv + jump_exists * pixel_in_first_half_tile;
-}
-
-vec2 convert_phosphor_tile_uv_wrap_to_tex_uv(const vec2 tile_uv_wrap,
-    const vec4 mask_tile_start_uv_and_size)
-{
-    //  Requires:   1.) tile_uv_wrap contains tile-relative uv coords, where the
-    //                  tile spans from [0, 1], such that (0.5, 0.5) is at the
-    //                  tile center.  The input coords can range from [0, inf],
-    //                  and their fractional parts map to a repeated tile.
-    //                  ("Tile" can mean texture, the video embedded in the
-    //                  texture, or some other "tile" embedded in a texture.)
-    //              2.) mask_tile_start_uv_and_size.xy contains tex_uv coords
-    //                  for the start of the embedded tile in the full texture.
-    //              3.) mask_tile_start_uv_and_size.zw contains the [fractional]
-    //                  tex_uv size of the embedded tile in the full texture.
-    //  Returns:    Return tex_uv coords (used for texture sampling)
-    //              corresponding to tile_uv_wrap.
-    if(get_mask_sample_mode() < 0.5)
-    {
-        //  Manually repeat the resized mask tile to fill the screen:
-        //  First get fractional tile_uv coords.  Using frac/fmod on coords
-        //  confuses anisotropic filtering; fix it as user options dictate.
-        //  derived-settings-and-constants.h disables incompatible options.
-        #ifdef ANISOTROPIC_TILING_COMPAT_TILE_FLAT_TWICE
-            vec2 tile_uv = frac(tile_uv_wrap * 0.5) * 2.0;
-        #else
-            vec2 tile_uv = frac(tile_uv_wrap);
-        #endif
-        #ifdef ANISOTROPIC_TILING_COMPAT_FIX_DISCONTINUITIES
-            const vec2 tile_uv_dx = ddx(tile_uv);
-            const vec2 tile_uv_dy = ddy(tile_uv);
-            tile_uv = fix_tiling_discontinuities_normalized(tile_uv,
-                tile_uv_dx, tile_uv_dy);
-        #endif
-        //  The tile is embedded in a padded FBO, and it may start at a
-        //  nonzero offset if border texels are used to avoid artifacts:
-        const vec2 mask_tex_uv = mask_tile_start_uv_and_size.xy +
-            tile_uv * mask_tile_start_uv_and_size.zw;
-        return mask_tex_uv;
-    }
-    else
-    {
-        //  Sample from the input phosphor mask texture with hardware tiling.
-        //  If we're tiling at the original size (mode 2), the "tile" is the
-        //  whole texture, and it contains a large number of triads mapped with
-        //  a 1:1 pixel:texel ratio.  OTHERWISE, the texture contains a single
-        //  unresized tile.  tile_uv_wrap already has correct coords for both!
-        return tile_uv_wrap;
-    }
-}
-
 
 #endif  //  PHOSPHOR_MASK_RESIZING_H
 

crt/shaders/crt-royale/src/scanline-functions.h

@@ -1,3 +1,6 @@
+#ifndef SCANLINE_FUNCTIONS_H
+#define SCANLINE_FUNCTIONS_H
+
 /////////////////////////////  GPL LICENSE NOTICE  /////////////////////////////
 
 //  crt-royale: A full-featured CRT shader, with cheese.
@@ -350,4 +353,6 @@ vec3 scanline_contrib(vec3 dist, vec3 color,
                 dist, color, pixel_height, sigma_range);
         }
     }
-}
+}
+
+#endif  //  SCANLINE_FUNCTIONS_H

crt/shaders/crt-royale/src/special-functions.h

@@ -1,3 +1,7 @@
+#ifndef SPECIAL_FUNCTIONS_H
+#define SPECIAL_FUNCTIONS_H
+
+
 /////////////////////////////////  MIT LICENSE  ////////////////////////////////
 
 //  Copyright (C) 2014 TroggleMonkey
@@ -489,4 +493,6 @@ float normalized_ligamma(const float s, const float z)
 	const float s_inv = 1.0/s;
 	const float gamma_s_inv = 1.0/gamma_impl(s, s_inv);
 	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
-}
+}
+
+#endif  //  SPECIAL_FUNCTIONS_H

include/blur-functions.h

@@ -0,0 +1,281 @@
+#define BLUR_FUNCTIONS
+
+/////////////////////////////////  MIT LICENSE  ////////////////////////////////
+
+//  Copyright (C) 2014 TroggleMonkey
+//
+//  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.
+
+/////////////////////////////////  DESCRIPTION  ////////////////////////////////
+
+//  This file provides reusable one-pass and separable (two-pass) blurs.
+//  Requires:   All blurs share these requirements (dxdy requirement is split):
+//              1.) All requirements of gamma-management.h must be satisfied!
+//              2.) filter_linearN must == "true" in your .cgp preset unless
+//                  you're using tex2DblurNresize at 1x scale.
+//              3.) mipmap_inputN must == "true" in your .cgp preset if
+//                  IN.output_size < IN.video_size.
+//              4.) IN.output_size == IN.video_size / pow(2, M), where M is some
+//                  positive integer.  tex2Dblur*resize can resize arbitrarily
+//                  (and the blur will be done after resizing), but arbitrary
+//                  resizes "fail" with other blurs due to the way they mix
+//                  static weights with bilinear sample exploitation.
+//              5.) In general, dxdy should contain the uv pixel spacing:
+//                      dxdy = (IN.video_size/IN.output_size)/IN.texture_size
+//              6.) For separable blurs (tex2DblurNresize and tex2DblurNfast),
+//                  zero out the dxdy component in the unblurred dimension:
+//                      dxdy = vec2(dxdy.x, 0.0) or vec2(0.0, dxdy.y)
+//              Many blurs share these requirements:
+//              1.) One-pass blurs require scale_xN == scale_yN or scales > 1.0,
+//                  or they will blur more in the lower-scaled dimension.
+//              2.) One-pass shared sample blurs require ddx(), ddy(), and
+//                  tex2Dlod() to be supported by the current Cg profile, and
+//                  the drivers must support high-quality derivatives.
+//              3.) One-pass shared sample blurs require:
+//                      tex_uv.w == log2(IN.video_size/IN.output_size).y;
+//              Non-wrapper blurs share this requirement:
+//              1.) sigma is the intended standard deviation of the blur
+//              Wrapper blurs share this requirement, which is automatically
+//              met (unless OVERRIDE_BLUR_STD_DEVS is #defined; see below):
+//              1.) blurN_std_dev must be global static const float values
+//                  specifying standard deviations for Nx blurs in units
+//                  of destination pixels
+//  Optional:   1.) The including file (or an earlier included file) may
+//                  optionally #define USE_BINOMIAL_BLUR_STD_DEVS to replace
+//                  default standard deviations with those matching a binomial
+//                  distribution.  (See below for details/properties.)
+//              2.) The including file (or an earlier included file) may
+//                  optionally #define OVERRIDE_BLUR_STD_DEVS and override:
+//                      static const float blur3_std_dev
+//                      static const float blur4_std_dev
+//                      static const float blur5_std_dev
+//                      static const float blur6_std_dev
+//                      static const float blur7_std_dev
+//                      static const float blur8_std_dev
+//                      static const float blur9_std_dev
+//                      static const float blur10_std_dev
+//                      static const float blur11_std_dev
+//                      static const float blur12_std_dev
+//                      static const float blur17_std_dev
+//                      static const float blur25_std_dev
+//                      static const float blur31_std_dev
+//                      static const float blur43_std_dev
+//              3.) The including file (or an earlier included file) may
+//                  optionally #define OVERRIDE_ERROR_BLURRING and override:
+//                      static const float error_blurring
+//                  This tuning value helps mitigate weighting errors from one-
+//                  pass shared-sample blurs sharing bilinear samples between
+//                  fragments.  Values closer to 0.0 have "correct" blurriness
+//                  but allow more artifacts, and values closer to 1.0 blur away
+//                  artifacts by sampling closer to halfway between texels.
+//              UPDATE 6/21/14: The above static constants may now be overridden
+//              by non-static uniform constants.  This permits exposing blur
+//              standard deviations as runtime GUI shader parameters.  However,
+//              using them keeps weights from being statically computed, and the
+//              speed hit depends on the blur: On my machine, uniforms kill over
+//              53% of the framerate with tex2Dblur12x12shared, but they only
+//              drop the framerate by about 18% with tex2Dblur11fast.
+//  Quality and Performance Comparisons:
+//  For the purposes of the following discussion, "no sRGB" means
+//  GAMMA_ENCODE_EVERY_FBO is #defined, and "sRGB" means it isn't.
+//  1.) tex2DblurNfast is always faster than tex2DblurNresize.
+//  2.) tex2DblurNresize functions are the only ones that can arbitrarily resize
+//      well, because they're the only ones that don't exploit bilinear samples.
+//      This also means they're the only functions which can be truly gamma-
+//      correct without linear (or sRGB FBO) input, but only at 1x scale.
+//  3.) One-pass shared sample blurs only have a speed advantage without sRGB.
+//      They also have some inaccuracies due to their shared-[bilinear-]sample
+//      design, which grow increasingly bothersome for smaller blurs and higher-
+//      frequency source images (relative to their resolution).  I had high
+//      hopes for them, but their most realistic use case is limited to quickly
+//      reblurring an already blurred input at full resolution.  Otherwise:
+//      a.) If you're blurring a low-resolution source, you want a better blur.
+//      b.) If you're blurring a lower mipmap, you want a better blur.
+//      c.) If you're blurring a high-resolution, high-frequency source, you
+//          want a better blur.
+//  4.) The one-pass blurs without shared samples grow slower for larger blurs,
+//      but they're competitive with separable blurs at 5x5 and smaller, and
+//      even tex2Dblur7x7 isn't bad if you're wanting to conserve passes.
+//  Here are some framerates from a GeForce 8800GTS.  The first pass resizes to
+//  viewport size (4x in this test) and linearizes for sRGB codepaths, and the
+//  remaining passes perform 6 full blurs.  Mipmapped tests are performed at the
+//  same scale, so they just measure the cost of mipmapping each FBO (only every
+//  other FBO is mipmapped for separable blurs, to mimic realistic usage).
+//  Mipmap      Neither     sRGB+Mipmap sRGB        Function
+//  76.0        92.3        131.3       193.7       tex2Dblur3fast
+//  63.2        74.4        122.4       175.5       tex2Dblur3resize
+//  93.7        121.2       159.3       263.2       tex2Dblur3x3
+//  59.7        68.7        115.4       162.1       tex2Dblur3x3resize
+//  63.2        74.4        122.4       175.5       tex2Dblur5fast
+//  49.3        54.8        100.0       132.7       tex2Dblur5resize
+//  59.7        68.7        115.4       162.1       tex2Dblur5x5
+//  64.9        77.2        99.1        137.2       tex2Dblur6x6shared
+//  55.8        63.7        110.4       151.8       tex2Dblur7fast
+//  39.8        43.9        83.9        105.8       tex2Dblur7resize
+//  40.0        44.2        83.2        104.9       tex2Dblur7x7
+//  56.4        65.5        71.9        87.9        tex2Dblur8x8shared
+//  49.3        55.1        99.9        132.5       tex2Dblur9fast
+//  33.3        36.2        72.4        88.0        tex2Dblur9resize
+//  27.8        29.7        61.3        72.2        tex2Dblur9x9
+//  37.2        41.1        52.6        60.2        tex2Dblur10x10shared
+//  44.4        49.5        91.3        117.8       tex2Dblur11fast
+//  28.8        30.8        63.6        75.4        tex2Dblur11resize
+//  33.6        36.5        40.9        45.5        tex2Dblur12x12shared
+//  TODO: Fill in benchmarks for new untested blurs.
+//                                                  tex2Dblur17fast
+//                                                  tex2Dblur25fast
+//                                                  tex2Dblur31fast
+//                                                  tex2Dblur43fast
+//                                                  tex2Dblur3x3resize
+
+/////////////////////////////  SETTINGS MANAGEMENT  ////////////////////////////
+
+//  Set static standard deviations, but allow users to override them with their
+//  own constants (even non-static uniforms if they're okay with the speed hit):
+#ifndef OVERRIDE_BLUR_STD_DEVS
+    //  blurN_std_dev values are specified in terms of dxdy strides.
+    #ifdef USE_BINOMIAL_BLUR_STD_DEVS
+        //  By request, we can define standard deviations corresponding to a
+        //  binomial distribution with p = 0.5 (related to Pascal's triangle).
+        //  This distribution works such that blurring multiple times should
+        //  have the same result as a single larger blur.  These values are
+        //  larger than default for blurs up to 6x and smaller thereafter.
+        const float blur3_std_dev = 0.84931640625;
+        const float blur4_std_dev = 0.84931640625;
+        const float blur5_std_dev = 1.0595703125;
+        const float blur6_std_dev = 1.06591796875;
+        const float blur7_std_dev = 1.17041015625;
+        const float blur8_std_dev = 1.1720703125;
+        const float blur9_std_dev = 1.2259765625;
+        const float blur10_std_dev = 1.21982421875;
+        const float blur11_std_dev = 1.25361328125;
+        const float blur12_std_dev = 1.2423828125;
+        const float blur17_std_dev = 1.27783203125;
+        const float blur25_std_dev = 1.2810546875;
+        const float blur31_std_dev = 1.28125;
+        const float blur43_std_dev = 1.28125;
+    #else
+        //  The defaults are the largest values that keep the largest unused
+        //  blur term on each side <= 1.0/256.0.  (We could get away with more
+        //  or be more conservative, but this compromise is pretty reasonable.)
+        const float blur3_std_dev = 0.62666015625;
+        const float blur4_std_dev = 0.66171875;
+        const float blur5_std_dev = 0.9845703125;
+        const float blur6_std_dev = 1.02626953125;
+        const float blur7_std_dev = 1.36103515625;
+        const float blur8_std_dev = 1.4080078125;
+        const float blur9_std_dev = 1.7533203125;
+        const float blur10_std_dev = 1.80478515625;
+        const float blur11_std_dev = 2.15986328125;
+        const float blur12_std_dev = 2.215234375;
+        const float blur17_std_dev = 3.45535583496;
+        const float blur25_std_dev = 5.3409576416;
+        const float blur31_std_dev = 6.86488037109;
+        const float blur43_std_dev = 10.1852050781;
+    #endif  //  USE_BINOMIAL_BLUR_STD_DEVS
+#endif  //  OVERRIDE_BLUR_STD_DEVS
+
+#ifndef OVERRIDE_ERROR_BLURRING
+    //  error_blurring should be in [0.0, 1.0].  Higher values reduce ringing
+    //  in shared-sample blurs but increase blurring and feature shifting.
+    const float error_blurring = 0.5;
+#endif
+
+//  Make a length squared helper macro (for usage with static constants):
+#define LENGTH_SQ(vec) (dot(vec, vec))
+
+///////////////////////////////////  HELPERS  //////////////////////////////////
+
+vec4 uv2_to_uv4(vec2 tex_uv)
+{
+    //  Make a vec2 uv offset safe for adding to vec4 tex2Dlod coords:
+    return vec4(tex_uv, 0.0, 0.0);
+}
+
+//  Make a length squared helper macro (for usage with static constants):
+#define LENGTH_SQ(vec) (dot(vec, vec))
+
+float get_fast_gaussian_weight_sum_inv(const float sigma)
+{
+    //  We can use the Gaussian integral to calculate the asymptotic weight for
+    //  the center pixel.  Since the unnormalized center pixel weight is 1.0,
+    //  the normalized weight is the same as the weight sum inverse.  Given a
+    //  large enough blur (9+), the asymptotic weight sum is close and faster:
+    //      center_weight = 0.5 *
+    //          (erf(0.5/(sigma*sqrt(2.0))) - erf(-0.5/(sigma*sqrt(2.0))))
+    //      erf(-x) == -erf(x), so we get 0.5 * (2.0 * erf(blah blah)):
+    //  However, we can get even faster results with curve-fitting.  These are
+    //  also closer than the asymptotic results, because they were constructed
+    //  from 64 blurs sizes from [3, 131) and 255 equally-spaced sigmas from
+    //  (0, blurN_std_dev), so the results for smaller sigmas are biased toward
+    //  smaller blurs.  The max error is 0.0031793913.
+    //  Relative FPS: 134.3 with erf, 135.8 with curve-fitting.
+    //static const float temp = 0.5/sqrt(2.0);
+    //return erf(temp/sigma);
+    return min(exp(exp(0.348348412457428/
+        (sigma - 0.0860587260734721))), 0.399334576340352/sigma);
+}
+
+////////////////////  ARBITRARILY RESIZABLE ONE-PASS BLURS  ////////////////////
+
+vec3 tex2Dblur3x3resize(const sampler2D tex, const vec2 tex_uv,
+    const vec2 dxdy, const float sigma)
+{
+    //  Requires:   Global requirements must be met (see file description).
+    //  Returns:    A 3x3 Gaussian blurred mipmapped texture lookup of the
+    //              resized input.
+    //  Description:
+    //  This is the only arbitrarily resizable one-pass blur; tex2Dblur5x5resize
+    //  would perform like tex2Dblur9x9, MUCH slower than tex2Dblur5resize.
+    const float denom_inv = 0.5/(sigma*sigma);
+    //  Load each sample.  We need all 3x3 samples.  Quad-pixel communication
+    //  won't help either: This should perform like tex2Dblur5x5, but sharing a
+    //  4x4 sample field would perform more like tex2Dblur8x8shared (worse).
+    const vec2 sample4_uv = tex_uv;
+    const vec2 dx = vec2(dxdy.x, 0.0);
+    const vec2 dy = vec2(0.0, dxdy.y);
+    const vec2 sample1_uv = sample4_uv - dy;
+    const vec2 sample7_uv = sample4_uv + dy;
+    const vec3 sample0 = tex2D_linearize(tex, sample1_uv - dx).rgb;
+    const vec3 sample1 = tex2D_linearize(tex, sample1_uv).rgb;
+    const vec3 sample2 = tex2D_linearize(tex, sample1_uv + dx).rgb;
+    const vec3 sample3 = tex2D_linearize(tex, sample4_uv - dx).rgb;
+    const vec3 sample4 = tex2D_linearize(tex, sample4_uv).rgb;
+    const vec3 sample5 = tex2D_linearize(tex, sample4_uv + dx).rgb;
+    const vec3 sample6 = tex2D_linearize(tex, sample7_uv - dx).rgb;
+    const vec3 sample7 = tex2D_linearize(tex, sample7_uv).rgb;
+    const vec3 sample8 = tex2D_linearize(tex, sample7_uv + dx).rgb;
+    //  Statically compute Gaussian sample weights:
+    const float w4 = 1.0;
+    const float w1_3_5_7 = exp(-LENGTH_SQ(vec2(1.0, 0.0)) * denom_inv);
+    const float w0_2_6_8 = exp(-LENGTH_SQ(vec2(1.0, 1.0)) * denom_inv);
+    const float weight_sum_inv = 1.0/(w4 + 4.0 * (w1_3_5_7 + w0_2_6_8));
+    //  Weight and sum the samples:
+    const vec3 sum = w4 * sample4 +
+        w1_3_5_7 * (sample1 + sample3 + sample5 + sample7) +
+        w0_2_6_8 * (sample0 + sample2 + sample6 + sample8);
+    return sum * weight_sum_inv;
+}
+
+//  Resizable one-pass blurs:
+vec3 tex2Dblur3x3resize(const sampler2D texture, const vec2 tex_uv,
+    const vec2 dxdy)
+{
+    return tex2Dblur3x3resize(texture, tex_uv, dxdy, blur3_std_dev);
+}

include/gamma-management-old.h

@@ -0,0 +1,547 @@
+#ifndef GAMMA_MANAGEMENT_H
+#define GAMMA_MANAGEMENT_H
+
+/////////////////////////////////  MIT LICENSE  ////////////////////////////////
+
+//  Copyright (C) 2014 TroggleMonkey
+//
+//  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.
+
+/////////////////////////////////  DESCRIPTION  ////////////////////////////////
+
+//  This file provides gamma-aware tex*D*() and encode_output() functions.
+//  Requires:   Before #include-ing this file, the including file must #define
+//              the following macros when applicable and follow their rules:
+//              1.) #define FIRST_PASS if this is the first pass.
+//              2.) #define LAST_PASS if this is the last pass.
+//              3.) If sRGB is available, set srgb_framebufferN = "true" for
+//                  every pass except the last in your .cgp preset.
+//              4.) If sRGB isn't available but you want gamma-correctness with
+//                  no banding, #define GAMMA_ENCODE_EVERY_FBO each pass.
+//              5.) #define SIMULATE_CRT_ON_LCD if desired (precedence over 5-7)
+//              6.) #define SIMULATE_GBA_ON_LCD if desired (precedence over 6-7)
+//              7.) #define SIMULATE_LCD_ON_CRT if desired (precedence over 7)
+//              8.) #define SIMULATE_GBA_ON_CRT if desired (precedence over -)
+//              If an option in [5, 8] is #defined in the first or last pass, it
+//              should be #defined for both.  It shouldn't make a difference
+//              whether it's #defined for intermediate passes or not.
+//  Optional:   The including file (or an earlier included file) may optionally
+//              #define a number of macros indicating it will override certain
+//              macros and associated constants are as follows:
+//              static constants with either static or uniform constants.  The
+//              1.) OVERRIDE_STANDARD_GAMMA: The user must first define:
+//                  static const float ntsc_gamma
+//                  static const float pal_gamma
+//                  static const float crt_reference_gamma_high
+//                  static const float crt_reference_gamma_low
+//                  static const float lcd_reference_gamma
+//                  static const float crt_office_gamma
+//                  static const float lcd_office_gamma
+//              2.) OVERRIDE_DEVICE_GAMMA: The user must first define:
+//                  static const float crt_gamma
+//                  static const float gba_gamma
+//                  static const float lcd_gamma
+//              3.) OVERRIDE_FINAL_GAMMA: The user must first define:
+//                  static const float input_gamma
+//                  static const float intermediate_gamma
+//                  static const float output_gamma
+//                  (intermediate_gamma is for GAMMA_ENCODE_EVERY_FBO.)
+//              4.) OVERRIDE_ALPHA_ASSUMPTIONS: The user must first define:
+//                  static const bool assume_opaque_alpha
+//              The gamma constant overrides must be used in every pass or none,
+//              and OVERRIDE_FINAL_GAMMA bypasses all of the SIMULATE* macros.
+//              OVERRIDE_ALPHA_ASSUMPTIONS may be set on a per-pass basis.
+//  Usage:      After setting macros appropriately, ignore gamma correction and
+//              replace all tex*D*() calls with equivalent gamma-aware
+//              tex*D*_linearize calls, except:
+//              1.) When you read an LUT, use regular tex*D or a gamma-specified
+//                  function, depending on its gamma encoding:
+//                      tex*D*_linearize_gamma (takes a runtime gamma parameter)
+//              2.) If you must read pass0's original input in a later pass, use
+//                  tex2D_linearize_ntsc_gamma.  If you want to read pass0's
+//                  input with gamma-corrected bilinear filtering, consider
+//                  creating a first linearizing pass and reading from the input
+//                  of pass1 later.
+//              Then, return encode_output(color) from every fragment shader.
+//              Finally, use the global gamma_aware_bilinear boolean if you want
+//              to statically branch based on whether bilinear filtering is
+//              gamma-correct or not (e.g. for placing Gaussian blur samples).
+//
+//  Detailed Policy:
+//  tex*D*_linearize() functions enforce a consistent gamma-management policy
+//  based on the FIRST_PASS and GAMMA_ENCODE_EVERY_FBO settings.  They assume
+//  their input texture has the same encoding characteristics as the input for
+//  the current pass (which doesn't apply to the exceptions listed above).
+//  Similarly, encode_output() enforces a policy based on the LAST_PASS and
+//  GAMMA_ENCODE_EVERY_FBO settings.  Together, they result in one of the
+//  following two pipelines.
+//  Typical pipeline with intermediate sRGB framebuffers:
+//      linear_color = pow(pass0_encoded_color, input_gamma);
+//      intermediate_output = linear_color;     //  Automatic sRGB encoding
+//      linear_color = intermediate_output;     //  Automatic sRGB decoding
+//      final_output = pow(intermediate_output, 1.0/output_gamma);
+//  Typical pipeline without intermediate sRGB framebuffers:
+//      linear_color = pow(pass0_encoded_color, input_gamma);
+//      intermediate_output = pow(linear_color, 1.0/intermediate_gamma);
+//      linear_color = pow(intermediate_output, intermediate_gamma);
+//      final_output = pow(intermediate_output, 1.0/output_gamma);
+//  Using GAMMA_ENCODE_EVERY_FBO is much slower, but it's provided as a way to
+//  easily get gamma-correctness without banding on devices where sRGB isn't
+//  supported.
+//
+//  Use This Header to Maximize Code Reuse:
+//  The purpose of this header is to provide a consistent interface for texture
+//  reads and output gamma-encoding that localizes and abstracts away all the
+//  annoying details.  This greatly reduces the amount of code in each shader
+//  pass that depends on the pass number in the .cgp preset or whether sRGB
+//  FBO's are being used: You can trivially change the gamma behavior of your
+//  whole pass by commenting or uncommenting 1-3 #defines.  To reuse the same
+//  code in your first, Nth, and last passes, you can even put it all in another
+//  header file and #include it from skeleton .cg files that #define the
+//  appropriate pass-specific settings.
+//
+//  Rationale for Using Three Macros:
+//  This file uses GAMMA_ENCODE_EVERY_FBO instead of an opposite macro like
+//  SRGB_PIPELINE to ensure sRGB is assumed by default, which hopefully imposes
+//  a lower maintenance burden on each pass.  At first glance it seems we could
+//  accomplish everything with two macros: GAMMA_CORRECT_IN / GAMMA_CORRECT_OUT.
+//  This works for simple use cases where input_gamma == output_gamma, but it
+//  breaks down for more complex scenarios like CRT simulation, where the pass
+//  number determines the gamma encoding of the input and output.
+
+
+///////////////////////////////  BASE CONSTANTS  ///////////////////////////////
+
+//  Set standard gamma constants, but allow users to override them:
+#ifndef OVERRIDE_STANDARD_GAMMA
+    //  Standard encoding gammas:
+    const float ntsc_gamma = 2.2;    //  Best to use NTSC for PAL too?
+    const float pal_gamma = 2.8;     //  Never actually 2.8 in practice
+    //  Typical device decoding gammas (only use for emulating devices):
+    //  CRT/LCD reference gammas are higher than NTSC and Rec.709 video standard
+    //  gammas: The standards purposely undercorrected for an analog CRT's
+    //  assumed 2.5 reference display gamma to maintain contrast in assumed
+    //  [dark] viewing conditions: http://www.poynton.com/PDFs/GammaFAQ.pdf
+    //  These unstated assumptions about display gamma and perceptual rendering
+    //  intent caused a lot of confusion, and more modern CRT's seemed to target
+    //  NTSC 2.2 gamma with circuitry.  LCD displays seem to have followed suit
+    //  (they struggle near black with 2.5 gamma anyway), especially PC/laptop
+    //  displays designed to view sRGB in bright environments.  (Standards are
+    //  also in flux again with BT.1886, but it's underspecified for displays.)
+    const float crt_reference_gamma_high = 2.5;  //  In (2.35, 2.55)
+    const float crt_reference_gamma_low = 2.35;  //  In (2.35, 2.55)
+    const float lcd_reference_gamma = 2.5;       //  To match CRT
+    const float crt_office_gamma = 2.2;  //  Circuitry-adjusted for NTSC
+    const float lcd_office_gamma = 2.2;  //  Approximates sRGB
+#endif  //  OVERRIDE_STANDARD_GAMMA
+
+//  Assuming alpha == 1.0 might make it easier for users to avoid some bugs,
+//  but only if they're aware of it.
+#ifndef OVERRIDE_ALPHA_ASSUMPTIONS
+    const bool assume_opaque_alpha = false;
+#endif
+
+
+///////////////////////  DERIVED CONSTANTS AS FUNCTIONS  ///////////////////////
+
+//  gamma-management.h should be compatible with overriding gamma values with
+//  runtime user parameters, but we can only define other global constants in
+//  terms of static constants, not uniform user parameters.  To get around this
+//  limitation, we need to define derived constants using functions.
+
+//  Set device gamma constants, but allow users to override them:
+#ifdef OVERRIDE_DEVICE_GAMMA
+    //  The user promises to globally define the appropriate constants:
+    float get_crt_gamma()    {   return crt_gamma;   }
+    float get_gba_gamma()    {   return gba_gamma;   }
+    float get_lcd_gamma()    {   return lcd_gamma;   }
+#else
+    float get_crt_gamma()    {   return crt_reference_gamma_high;    }
+    float get_gba_gamma()    {   return 3.5; }   //  Game Boy Advance; in (3.0, 4.0)
+    float get_lcd_gamma()    {   return lcd_office_gamma;            }
+#endif  //  OVERRIDE_DEVICE_GAMMA
+
+//  Set decoding/encoding gammas for the first/lass passes, but allow overrides:
+#ifdef OVERRIDE_FINAL_GAMMA
+    //  The user promises to globally define the appropriate constants:
+    float get_intermediate_gamma()   {   return intermediate_gamma;  }
+    float get_input_gamma()          {   return input_gamma;         }
+    float get_output_gamma()         {   return output_gamma;        }
+#else
+    //  If we gamma-correct every pass, always use ntsc_gamma between passes to
+    //  ensure middle passes don't need to care if anything is being simulated:
+    float get_intermediate_gamma()   {   return ntsc_gamma;          }
+    #ifdef SIMULATE_CRT_ON_LCD
+        float get_input_gamma()      {   return get_crt_gamma();     }
+        float get_output_gamma()     {   return get_lcd_gamma();     }
+    #else
+    #ifdef SIMULATE_GBA_ON_LCD
+        float get_input_gamma()      {   return get_gba_gamma();     }
+        float get_output_gamma()     {   return get_lcd_gamma();     }
+    #else
+    #ifdef SIMULATE_LCD_ON_CRT
+        float get_input_gamma()      {   return get_lcd_gamma();     }
+        float get_output_gamma()     {   return get_crt_gamma();     }
+    #else
+    #ifdef SIMULATE_GBA_ON_CRT
+        float get_input_gamma()      {   return get_gba_gamma();     }
+        float get_output_gamma()     {   return get_crt_gamma();     }
+    #else   //  Don't simulate anything:
+        float get_input_gamma()      {   return ntsc_gamma;          }
+        float get_output_gamma()     {   return ntsc_gamma;          }
+    #endif  //  SIMULATE_GBA_ON_CRT
+    #endif  //  SIMULATE_LCD_ON_CRT
+    #endif  //  SIMULATE_GBA_ON_LCD
+    #endif  //  SIMULATE_CRT_ON_LCD
+#endif  //  OVERRIDE_FINAL_GAMMA
+
+//  Set decoding/encoding gammas for the current pass.  Use static constants for
+//  linearize_input and gamma_encode_output, because they aren't derived, and
+//  they let the compiler do dead-code elimination.
+#ifndef GAMMA_ENCODE_EVERY_FBO
+    #ifdef FIRST_PASS
+        const bool linearize_input = true;
+        float get_pass_input_gamma()     {   return get_input_gamma();   }
+    #else
+        const bool linearize_input = false;
+        float get_pass_input_gamma()     {   return 1.0;                 }
+    #endif
+    #ifdef LAST_PASS
+        const bool gamma_encode_output = true;
+        float get_pass_output_gamma()    {   return get_output_gamma();  }
+    #else
+        const bool gamma_encode_output = false;
+        float get_pass_output_gamma()    {   return 1.0;                 }
+    #endif
+#else
+    const bool linearize_input = true;
+    const bool gamma_encode_output = true;
+    #ifdef FIRST_PASS
+        float get_pass_input_gamma()     {   return get_input_gamma();   }
+    #else
+        float get_pass_input_gamma()     {   return get_intermediate_gamma();    }
+    #endif
+    #ifdef LAST_PASS
+        float get_pass_output_gamma()    {   return get_output_gamma();  }
+    #else
+        float get_pass_output_gamma()    {   return get_intermediate_gamma();    }
+    #endif
+#endif
+
+//  Users might want to know if bilinear filtering will be gamma-correct:
+const bool gamma_aware_bilinear = !linearize_input;
+
+
+//////////////////////  COLOR ENCODING/DECODING FUNCTIONS  /////////////////////
+
+vec4 encode_output(const vec4 color)
+{
+    if(gamma_encode_output)
+    {
+        if(assume_opaque_alpha)
+        {
+            return vec4(pow(color.rgb, vec3(1.0/get_pass_output_gamma())), 1.0);
+        }
+        else
+        {
+            return vec4(pow(color.rgb, vec3(1.0/get_pass_output_gamma())), color.a);
+        }
+    }
+    else
+    {
+        return color;
+    }
+}
+
+vec4 decode_input(const vec4 color)
+{
+    if(linearize_input)
+    {
+        if(assume_opaque_alpha)
+        {
+            return vec4(pow(color.rgb, vec3(get_pass_input_gamma())), 1.0);
+        }
+        else
+        {
+            return vec4(pow(color.rgb, vec3(get_pass_input_gamma())), color.a);
+        }
+    }
+    else
+    {
+        return color;
+    }
+}
+
+vec4 decode_gamma_input(const vec4 color, const vec3 gamma)
+{
+    if(assume_opaque_alpha)
+    {
+        return vec4(pow(color.rgb, vec3(gamma)), 1.0);
+    }
+    else
+    {
+        return vec4(pow(color.rgb, vec3(gamma)), color.a);
+    }
+}
+
+
+///////////////////////////  TEXTURE LOOKUP WRAPPERS  //////////////////////////
+
+//  "SMART" LINEARIZING TEXTURE LOOKUP FUNCTIONS:
+//  Provide a wide array of linearizing texture lookup wrapper functions.  The
+//  Cg shader spec Retroarch uses only allows for 2D textures, but 1D and 3D
+//  lookups are provided for completeness in case that changes someday.  Nobody
+//  is likely to use the *fetch and *proj functions, but they're included just
+//  in case.  The only tex*D texture sampling functions omitted are:
+//      - tex*Dcmpbias
+//      - tex*Dcmplod
+//      - tex*DARRAY*
+//      - tex*DMS*
+//      - Variants returning integers
+//  Standard line length restrictions are ignored below for vertical brevity.
+
+/*
+//  tex1D:
+vec4 tex1D_linearize(const sampler1D texture, const float tex_coords)
+{   return decode_input(tex1D(texture, tex_coords));   }
+
+vec4 tex1D_linearize(const sampler1D texture, const vec2 tex_coords)
+{   return decode_input(tex1D(texture, tex_coords));   }
+
+vec4 tex1D_linearize(const sampler1D texture, const float tex_coords, const int texel_off)
+{   return decode_input(tex1D(texture, tex_coords, texel_off));    }
+
+vec4 tex1D_linearize(const sampler1D texture, const vec2 tex_coords, const int texel_off)
+{   return decode_input(tex1D(texture, tex_coords, texel_off));    }
+
+vec4 tex1D_linearize(const sampler1D texture, const float tex_coords, const float dx, const float dy)
+{   return decode_input(tex1D(texture, tex_coords, dx, dy));   }
+
+vec4 tex1D_linearize(const sampler1D texture, const vec2 tex_coords, const float dx, const float dy)
+{   return decode_input(tex1D(texture, tex_coords, dx, dy));   }
+
+vec4 tex1D_linearize(const sampler1D texture, const float tex_coords, const float dx, const float dy, const int texel_off)
+{   return decode_input(tex1D(texture, tex_coords, dx, dy, texel_off));    }
+
+vec4 tex1D_linearize(const sampler1D texture, const vec2 tex_coords, const float dx, const float dy, const int texel_off)
+{   return decode_input(tex1D(texture, tex_coords, dx, dy, texel_off));    }
+
+//  tex1Dbias:
+vec4 tex1Dbias_linearize(const sampler1D texture, const vec4 tex_coords)
+{   return decode_input(tex1Dbias(texture, tex_coords));   }
+
+vec4 tex1Dbias_linearize(const sampler1D texture, const vec4 tex_coords, const int texel_off)
+{   return decode_input(tex1Dbias(texture, tex_coords, texel_off));    }
+
+//  tex1Dfetch:
+vec4 tex1Dfetch_linearize(const sampler1D texture, const int4 tex_coords)
+{   return decode_input(tex1Dfetch(texture, tex_coords));  }
+
+vec4 tex1Dfetch_linearize(const sampler1D texture, const int4 tex_coords, const int texel_off)
+{   return decode_input(tex1Dfetch(texture, tex_coords, texel_off));   }
+
+//  tex1Dlod:
+vec4 tex1Dlod_linearize(const sampler1D texture, const vec4 tex_coords)
+{   return decode_input(tex1Dlod(texture, tex_coords));    }
+
+vec4 tex1Dlod_linearize(const sampler1D texture, const vec4 tex_coords, const int texel_off)
+{   return decode_input(tex1Dlod(texture, tex_coords, texel_off));     }
+
+//  tex1Dproj:
+vec4 tex1Dproj_linearize(const sampler1D texture, const vec2 tex_coords)
+{   return decode_input(tex1Dproj(texture, tex_coords));   }
+
+vec4 tex1Dproj_linearize(const sampler1D texture, const vec3 tex_coords)
+{   return decode_input(tex1Dproj(texture, tex_coords));   }
+
+vec4 tex1Dproj_linearize(const sampler1D texture, const vec2 tex_coords, const int texel_off)
+{   return decode_input(tex1Dproj(texture, tex_coords, texel_off));    }
+
+vec4 tex1Dproj_linearize(const sampler1D texture, const vec3 tex_coords, const int texel_off)
+{   return decode_input(tex1Dproj(texture, tex_coords, texel_off));    }
+*/
+
+//  tex2D:
+vec4 tex2D_linearize(const sampler2D tex, const vec2 tex_coords)
+{   return decode_input(vec4(texture(tex, tex_coords)));   }
+
+vec4 tex2D_linearize(const sampler2D tex, const vec3 tex_coords)
+{   return decode_input(vec4(texture(tex, tex_coords)));   }
+
+vec4 tex2D_linearize(const sampler2D tex, const vec2 tex_coords, const int texel_off)
+{   return decode_input(vec4(texture(tex, tex_coords, texel_off)));    }
+
+vec4 tex2D_linearize(const sampler2D tex, const vec3 tex_coords, const int texel_off)
+{   return decode_input(vec4(texture(tex, tex_coords, texel_off)));    }
+
+vec4 tex2D_linearize(const sampler2D tex, const vec2 tex_coords, const vec2 dx, const vec2 dy)
+{   return decode_input(vec4(texture(tex, tex_coords, dx, dy)));   }
+
+vec4 tex2D_linearize(const sampler2D tex, const vec3 tex_coords, const vec2 dx, const vec2 dy)
+{   return decode_input(vec4(texture(tex, tex_coords, dx, dy)));   }
+
+vec4 tex2D_linearize(const sampler2D tex, const vec2 tex_coords, const vec2 dx, const vec2 dy, const int texel_off)
+{   return decode_input(vec4(texture(tex, tex_coords, dx, dy, texel_off)));    }
+
+vec4 tex2D_linearize(const sampler2D tex, const vec3 tex_coords, const vec2 dx, const vec2 dy, const int texel_off)
+{   return decode_input(vec4(texture(tex, tex_coords, dx, dy, texel_off)));    }
+
+//  tex2Dbias:
+vec4 tex2Dbias_linearize(const sampler2D tex, const vec4 tex_coords)
+{   return decode_input(vec4(tex2Dbias(tex, tex_coords)));   }
+
+vec4 tex2Dbias_linearize(const sampler2D tex, const vec4 tex_coords, const int texel_off)
+{   return decode_input(vec4(tex2Dbias(tex, tex_coords, texel_off)));    }
+
+//  tex2Dfetch:
+vec4 tex2Dfetch_linearize(const sampler2D tex, const ivec4 tex_coords)
+{   return decode_input(vec4(texture2Dfetch(tex, tex_coords)));  }
+
+vec4 tex2Dfetch_linearize(const sampler2D tex, const ivec4 tex_coords, const int texel_off)
+{   return decode_input(vec4(texture2Dfetch(tex, tex_coords, texel_off)));   }
+
+//  tex2Dlod:
+vec4 tex2Dlod_linearize(const sampler2D tex, const vec4 tex_coords)
+{   return decode_input(vec4(texture2Dlod(tex, tex_coords)));    }
+
+vec4 tex2Dlod_linearize(const sampler2D tex, const vec4 tex_coords, const int texel_off)
+{   return decode_input(vec4(texture2Dlod(tex, tex_coords, texel_off)));     }
+
+//  tex2Dproj:
+vec4 tex2Dproj_linearize(const sampler2D tex, const vec3 tex_coords)
+{   return decode_input(vec4(tex2Dproj(tex, tex_coords)));   }
+
+vec4 tex2Dproj_linearize(const sampler2D tex, const vec4 tex_coords)
+{   return decode_input(vec4(tex2Dproj(tex, tex_coords)));   }
+
+vec4 tex2Dproj_linearize(const sampler2D tex, const vec3 tex_coords, const int texel_off)
+{   return decode_input(vec4(tex2Dproj(tex, tex_coords, texel_off)));    }
+
+vec4 tex2Dproj_linearize(const sampler2D tex, const vec4 tex_coords, const int texel_off)
+{   return decode_input(vec4(tex2Dproj(tex, tex_coords, texel_off)));    }
+
+/*
+//  tex3D:
+vec4 tex3D_linearize(const sampler3D texture, const vec3 tex_coords)
+{   return decode_input(tex3D(texture, tex_coords));   }
+
+vec4 tex3D_linearize(const sampler3D texture, const vec3 tex_coords, const int texel_off)
+{   return decode_input(tex3D(texture, tex_coords, texel_off));    }
+
+vec4 tex3D_linearize(const sampler3D texture, const vec3 tex_coords, const vec3 dx, const vec3 dy)
+{   return decode_input(tex3D(texture, tex_coords, dx, dy));   }
+
+vec4 tex3D_linearize(const sampler3D texture, const vec3 tex_coords, const vec3 dx, const vec3 dy, const int texel_off)
+{   return decode_input(tex3D(texture, tex_coords, dx, dy, texel_off));    }
+
+//  tex3Dbias:
+vec4 tex3Dbias_linearize(const sampler3D texture, const vec4 tex_coords)
+{   return decode_input(tex3Dbias(texture, tex_coords));   }
+
+vec4 tex3Dbias_linearize(const sampler3D texture, const vec4 tex_coords, const int texel_off)
+{   return decode_input(tex3Dbias(texture, tex_coords, texel_off));    }
+
+//  tex3Dfetch:
+vec4 tex3Dfetch_linearize(const sampler3D texture, const int4 tex_coords)
+{   return decode_input(tex3Dfetch(texture, tex_coords));  }
+
+vec4 tex3Dfetch_linearize(const sampler3D texture, const int4 tex_coords, const int texel_off)
+{   return decode_input(tex3Dfetch(texture, tex_coords, texel_off));   }
+
+//  tex3Dlod:
+vec4 tex3Dlod_linearize(const sampler3D texture, const vec4 tex_coords)
+{   return decode_input(tex3Dlod(texture, tex_coords));    }
+
+vec4 tex3Dlod_linearize(const sampler3D texture, const vec4 tex_coords, const int texel_off)
+{   return decode_input(tex3Dlod(texture, tex_coords, texel_off));     }
+
+//  tex3Dproj:
+vec4 tex3Dproj_linearize(const sampler3D texture, const vec4 tex_coords)
+{   return decode_input(tex3Dproj(texture, tex_coords));   }
+
+vec4 tex3Dproj_linearize(const sampler3D texture, const vec4 tex_coords, const int texel_off)
+{   return decode_input(tex3Dproj(texture, tex_coords, texel_off));    }
+*/
+
+
+//  NONSTANDARD "SMART" LINEARIZING TEXTURE LOOKUP FUNCTIONS:
+//  This narrow selection of nonstandard tex2D* functions can be useful:
+
+//  tex2Dlod0: Automatically fill in the tex2D LOD parameter for mip level 0.
+vec4 tex2Dlod0_linearize(const sampler2D texture, const vec2 tex_coords)
+{   return decode_input(vec4(texture2Dlod(texture, vec4(tex_coords, 0.0, 0.0))));   }
+
+vec4 tex2Dlod0_linearize(const sampler2D texture, const vec2 tex_coords, const int texel_off)
+{   return decode_input(vec4(texture2Dlod(texture, vec4(tex_coords, 0.0, 0.0), texel_off)));    }
+
+
+//  MANUALLY LINEARIZING TEXTURE LOOKUP FUNCTIONS:
+//  Provide a narrower selection of tex2D* wrapper functions that decode an
+//  input sample with a specified gamma value.  These are useful for reading
+//  LUT's and for reading the input of pass0 in a later pass.
+
+//  tex2D:
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec2 tex_coords, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords), vec3(gamma)));   }
+
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec3 tex_coords, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords), vec3(gamma)));   }
+
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec2 tex_coords, const int texel_off, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords, texel_off), vec3(gamma)));    }
+
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec3 tex_coords, const int texel_off, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords, texel_off), vec3(gamma)));    }
+
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec2 tex_coords, const vec2 dx, const vec2 dy, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords, dx, dy), vec3(gamma)));   }
+
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec3 tex_coords, const vec2 dx, const vec2 dy, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords, dx, dy), vec3(gamma)));   }
+
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec2 tex_coords, const vec2 dx, const vec2 dy, const int texel_off, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords, dx, dy, texel_off), vec3(gamma)));    }
+
+vec4 tex2D_linearize_gamma(const sampler2D tex, const vec3 tex_coords, const vec2 dx, const vec2 dy, const int texel_off, const vec3 gamma)
+{   return decode_gamma_input(vec4(texture(tex, tex_coords, dx, dy, texel_off), vec3(gamma)));    }
+
+//  tex2Dbias:
+vec4 tex2Dbias_linearize_gamma(const sampler2D tex, const vec4 tex_coords, const vec3 gamma)
+{   return decode_gamma_input(vec4(tex2Dbias(tex, tex_coords), vec3(gamma)));   }
+
+vec4 tex2Dbias_linearize_gamma(const sampler2D tex, const vec4 tex_coords, const int texel_off, const vec3 gamma)
+{   return decode_gamma_input(vec4(tex2Dbias(tex, tex_coords, texel_off), vec3(gamma)));    }
+
+//  tex2Dfetch:
+vec4 tex2Dfetch_linearize_gamma(const sampler2D tex, const int4 tex_coords, const vec3 gamma)
+{   return decode_gamma_input(vec4(tex2Dfetch(tex, tex_coords), vec3(gamma)));  }
+
+vec4 tex2Dfetch_linearize_gamma(const sampler2D tex, const int4 tex_coords, const int texel_off, const vec3 gamma)
+{   return decode_gamma_input(vec4(tex2Dfetch(tex, tex_coords, texel_off), vec3(gamma)));   }
+
+//  tex2Dlod:
+vec4 tex2Dlod_linearize_gamma(const sampler2D tex, const vec4 tex_coords, const vec3 gamma)
+{   return decode_gamma_input(vec4(tex2Dlod(tex, tex_coords), vec3(gamma)));    }
+
+vec4 tex2Dlod_linearize_gamma(const sampler2D tex, const vec4 tex_coords, const int texel_off, const vec3 gamma)
+{   return decode_gamma_input(vec4(tex2Dlod(tex, tex_coords, texel_off), vec3(gamma)));     }
+
+
+#endif  //  GAMMA_MANAGEMENT_H
+

include/gamma-management.h

@@ -0,0 +1,160 @@
+///////////////////////////////  BASE CONSTANTS  ///////////////////////////////
+
+//  Set standard gamma constants, but allow users to override them:
+#ifndef OVERRIDE_STANDARD_GAMMA
+    //  Standard encoding gammas:
+    const float ntsc_gamma = 2.2;    //  Best to use NTSC for PAL too?
+    const float pal_gamma = 2.8;     //  Never actually 2.8 in practice
+    //  Typical device decoding gammas (only use for emulating devices):
+    //  CRT/LCD reference gammas are higher than NTSC and Rec.709 video standard
+    //  gammas: The standards purposely undercorrected for an analog CRT's
+    //  assumed 2.5 reference display gamma to maintain contrast in assumed
+    //  [dark] viewing conditions: http://www.poynton.com/PDFs/GammaFAQ.pdf
+    //  These unstated assumptions about display gamma and perceptual rendering
+    //  intent caused a lot of confusion, and more modern CRT's seemed to target
+    //  NTSC 2.2 gamma with circuitry.  LCD displays seem to have followed suit
+    //  (they struggle near black with 2.5 gamma anyway), especially PC/laptop
+    //  displays designed to view sRGB in bright environments.  (Standards are
+    //  also in flux again with BT.1886, but it's underspecified for displays.)
+    const float crt_reference_gamma_high = 2.5;  //  In (2.35, 2.55)
+    const float crt_reference_gamma_low = 2.35;  //  In (2.35, 2.55)
+    const float lcd_reference_gamma = 2.5;       //  To match CRT
+    const float crt_office_gamma = 2.2;  //  Circuitry-adjusted for NTSC
+    const float lcd_office_gamma = 2.2;  //  Approximates sRGB
+#endif  //  OVERRIDE_STANDARD_GAMMA
+
+//  Assuming alpha == 1.0 might make it easier for users to avoid some bugs,
+//  but only if they're aware of it.
+#ifndef OVERRIDE_ALPHA_ASSUMPTIONS
+    const bool assume_opaque_alpha = false;
+#endif
+
+
+///////////////////////  DERIVED CONSTANTS AS FUNCTIONS  ///////////////////////
+
+//  gamma-management.h should be compatible with overriding gamma values with
+//  runtime user parameters, but we can only define other global constants in
+//  terms of static constants, not uniform user parameters.  To get around this
+//  limitation, we need to define derived constants using functions.
+
+//  Set device gamma constants, but allow users to override them:
+#ifdef OVERRIDE_DEVICE_GAMMA
+    //  The user promises to globally define the appropriate constants:
+    float get_crt_gamma()    {   return crt_gamma;   }
+    float get_gba_gamma()    {   return gba_gamma;   }
+    float get_lcd_gamma()    {   return lcd_gamma;   }
+#else
+    float get_crt_gamma()    {   return crt_reference_gamma_high;    }
+    float get_gba_gamma()    {   return 3.5; }   //  Game Boy Advance; in (3.0, 4.0)
+    float get_lcd_gamma()    {   return lcd_office_gamma;            }
+#endif  //  OVERRIDE_DEVICE_GAMMA
+
+//  Set decoding/encoding gammas for the first/lass passes, but allow overrides:
+#ifdef OVERRIDE_FINAL_GAMMA
+    //  The user promises to globally define the appropriate constants:
+    float get_intermediate_gamma()   {   return intermediate_gamma;  }
+    float get_input_gamma()          {   return input_gamma;         }
+    float get_output_gamma()         {   return output_gamma;        }
+#else
+    //  If we gamma-correct every pass, always use ntsc_gamma between passes to
+    //  ensure middle passes don't need to care if anything is being simulated:
+    float get_intermediate_gamma()   {   return ntsc_gamma;          }
+    #ifdef SIMULATE_CRT_ON_LCD
+        float get_input_gamma()      {   return get_crt_gamma();     }
+        float get_output_gamma()     {   return get_lcd_gamma();     }
+    #else
+    #ifdef SIMULATE_GBA_ON_LCD
+        float get_input_gamma()      {   return get_gba_gamma();     }
+        float get_output_gamma()     {   return get_lcd_gamma();     }
+    #else
+    #ifdef SIMULATE_LCD_ON_CRT
+        float get_input_gamma()      {   return get_lcd_gamma();     }
+        float get_output_gamma()     {   return get_crt_gamma();     }
+    #else
+    #ifdef SIMULATE_GBA_ON_CRT
+        float get_input_gamma()      {   return get_gba_gamma();     }
+        float get_output_gamma()     {   return get_crt_gamma();     }
+    #else   //  Don't simulate anything:
+        float get_input_gamma()      {   return ntsc_gamma;          }
+        float get_output_gamma()     {   return ntsc_gamma;          }
+    #endif  //  SIMULATE_GBA_ON_CRT
+    #endif  //  SIMULATE_LCD_ON_CRT
+    #endif  //  SIMULATE_GBA_ON_LCD
+    #endif  //  SIMULATE_CRT_ON_LCD
+#endif  //  OVERRIDE_FINAL_GAMMA
+
+#ifndef GAMMA_ENCODE_EVERY_FBO
+    #ifdef FIRST_PASS
+        const bool linearize_input = true;
+        float get_pass_input_gamma()     {   return get_input_gamma();   }
+    #else
+        const bool linearize_input = false;
+        float get_pass_input_gamma()     {   return 1.0;                 }
+    #endif
+    #ifdef LAST_PASS
+        const bool gamma_encode_output = true;
+        float get_pass_output_gamma()    {   return get_output_gamma();  }
+    #else
+        const bool gamma_encode_output = false;
+        float get_pass_output_gamma()    {   return 1.0;                 }
+    #endif
+#else
+    const bool linearize_input = true;
+    const bool gamma_encode_output = true;
+    #ifdef FIRST_PASS
+        float get_pass_input_gamma()     {   return get_input_gamma();   }
+    #else
+        float get_pass_input_gamma()     {   return get_intermediate_gamma();    }
+    #endif
+    #ifdef LAST_PASS
+        float get_pass_output_gamma()    {   return get_output_gamma();  }
+    #else
+        float get_pass_output_gamma()    {   return get_intermediate_gamma();    }
+    #endif
+#endif
+
+vec4 decode_input(const vec4 color)
+{
+    if(linearize_input)
+    {
+        if(assume_opaque_alpha)
+        {
+            return vec4(pow(color.rgb, vec3(get_pass_input_gamma())), 1.0);
+        }
+        else
+        {
+            return vec4(pow(color.rgb, vec3(get_pass_input_gamma())), color.a);
+        }
+    }
+    else
+    {
+        return color;
+    }
+}
+
+vec4 encode_output(const vec4 color)
+{
+    if(gamma_encode_output)
+    {
+        if(assume_opaque_alpha)
+        {
+            return vec4(pow(color.rgb, vec3(1.0/get_pass_output_gamma())), 1.0);
+        }
+        else
+        {
+            return vec4(pow(color.rgb, vec3(1.0/get_pass_output_gamma())), color.a);
+        }
+    }
+    else
+    {
+        return color;
+    }
+}
+
+#define tex2D_linearize(C, D) decode_input(vec4(texture(C, D)))
+//vec4 tex2D_linearize(const sampler2D tex, const vec2 tex_coords)
+//{   return decode_input(vec4(texture(tex, tex_coords)));   }
+
+//#define tex2D_linearize(C, D, E) decode_input(vec4(texture(C, D, E)))
+//vec4 tex2D_linearize(const sampler2D tex, const vec2 tex_coords, const int texel_off)
+//{   return decode_input(vec4(texture(tex, tex_coords, texel_off)));    }

include/quad-pixel-communication.h

@@ -0,0 +1,243 @@
+#ifndef QUAD_PIXEL_COMMUNICATION_H
+#define QUAD_PIXEL_COMMUNICATION_H
+
+/////////////////////////////////  MIT LICENSE  ////////////////////////////////
+
+//  Copyright (C) 2014 TroggleMonkey*
+//
+//  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.
+
+/////////////////////////////////  DISCLAIMER  /////////////////////////////////
+
+//  *This code was inspired by "Shader Amortization using Pixel Quad Message
+//  Passing" by Eric Penner, published in GPU Pro 2, Chapter VI.2.  My intent
+//  is not to plagiarize his fundamentally similar code and assert my own
+//  copyright, but the algorithmic helper functions require so little code that
+//  implementations can't vary by much except bugfixes and conventions.  I just
+//  wanted to license my own particular code here to avoid ambiguity and make it
+//  clear that as far as I'm concerned, people can do as they please with it.
+
+/////////////////////////////////  DESCRIPTION  ////////////////////////////////
+
+//  Given screen pixel numbers, derive a "quad vector" describing a fragment's
+//  position in its 2x2 pixel quad.  Given that vector, obtain the values of any
+//  variable at neighboring fragments.
+//  Requires:   Using this file in general requires:
+//              1.) ddx() and ddy() are present in the current Cg profile.
+//              2.) The GPU driver is using fine/high-quality derivatives.
+//                  Functions will give incorrect results if this is not true,
+//                  so a test function is included.
+
+
+/////////////////////  QUAD-PIXEL COMMUNICATION PRIMITIVES  ////////////////////
+
+vec4 get_quad_vector_naive(const vec4 output_pixel_num_wrt_uvxy)
+{
+    //  Requires:   Two measures of the current fragment's output pixel number
+    //              in the range ([0, IN.output_size.x), [0, IN.output_size.y)):
+    //              1.) output_pixel_num_wrt_uvxy.xy increase with uv coords.
+    //              2.) output_pixel_num_wrt_uvxy.zw increase with screen xy.
+    //  Returns:    Two measures of the fragment's position in its 2x2 quad:
+    //              1.) The .xy components are its 2x2 placement with respect to
+    //                  uv direction (the origin (0, 0) is at the top-left):
+    //                  top-left     = (-1.0, -1.0) top-right    = ( 1.0, -1.0)
+    //                  bottom-left  = (-1.0,  1.0) bottom-right = ( 1.0,  1.0)
+    //                  You need this to arrange/weight shared texture samples.
+    //              2.) The .zw components are its 2x2 placement with respect to
+    //                  screen xy direction (IN.position); the origin varies.
+    //                  quad_gather needs this measure to work correctly.
+    //              Note: quad_vector.zw = quad_vector.xy * vec2(
+    //                      ddx(output_pixel_num_wrt_uvxy.x),
+    //                      ddy(output_pixel_num_wrt_uvxy.y));
+    //  Caveats:    This function assumes the GPU driver always starts 2x2 pixel
+    //              quads at even pixel numbers.  This assumption can be wrong
+    //              for odd output resolutions (nondeterministically so).
+    const vec4 pixel_odd = frac(output_pixel_num_wrt_uvxy * 0.5) * 2.0;
+    const vec4 quad_vector = pixel_odd * 2.0 - vec4(1.0);
+    return quad_vector;
+}
+
+vec4 get_quad_vector(const vec4 output_pixel_num_wrt_uvxy)
+{
+    //  Requires:   Same as get_quad_vector_naive() (see that first).
+    //  Returns:    Same as get_quad_vector_naive() (see that first), but it's
+    //              correct even if the 2x2 pixel quad starts at an odd pixel,
+    //              which can occur at odd resolutions.
+    const vec4 quad_vector_guess =
+        get_quad_vector_naive(output_pixel_num_wrt_uvxy);
+    //  If quad_vector_guess.zw doesn't increase with screen xy, we know
+    //  the 2x2 pixel quad starts at an odd pixel:
+    const vec2 odd_start_mirror = 0.5 * vec2(ddx(quad_vector_guess.z),
+                                                ddy(quad_vector_guess.w));
+    return quad_vector_guess * odd_start_mirror.xyxy;
+}
+
+vec4 get_quad_vector(const vec2 output_pixel_num_wrt_uv)
+{
+    //  Requires:   1.) ddx() and ddy() are present in the current Cg profile.
+    //              2.) output_pixel_num_wrt_uv must increase with uv coords and
+    //                  measure the current fragment's output pixel number in:
+    //                      ([0, IN.output_size.x), [0, IN.output_size.y))
+    //  Returns:    Same as get_quad_vector_naive() (see that first), but it's
+    //              correct even if the 2x2 pixel quad starts at an odd pixel,
+    //              which can occur at odd resolutions.
+    //  Caveats:    This function requires less information than the version
+    //              taking a vec4, but it's potentially slower.
+    //  Do screen coords increase with or against uv?  Get the direction
+    //  with respect to (uv.x, uv.y) for (screen.x, screen.y) in {-1, 1}.
+    const vec2 screen_uv_mirror = vec2(ddx(output_pixel_num_wrt_uv.x),
+                                        ddy(output_pixel_num_wrt_uv.y));
+    const vec2 pixel_odd_wrt_uv = frac(output_pixel_num_wrt_uv * 0.5) * 2.0;
+    const vec2 quad_vector_uv_guess = (pixel_odd_wrt_uv - vec2(0.5)) * 2.0;
+    const vec2 quad_vector_screen_guess = quad_vector_uv_guess * screen_uv_mirror;
+    //  If quad_vector_screen_guess doesn't increase with screen xy, we know
+    //  the 2x2 pixel quad starts at an odd pixel:
+    const vec2 odd_start_mirror = 0.5 * vec2(ddx(quad_vector_screen_guess.x),
+                                                ddy(quad_vector_screen_guess.y));
+    const vec4 quad_vector_guess = vec4(
+        quad_vector_uv_guess, quad_vector_screen_guess);
+    return quad_vector_guess * odd_start_mirror.xyxy;
+}
+
+void quad_gather(const vec4 quad_vector, const vec4 curr,
+    out vec4 adjx, out vec4 adjy, out vec4 diag)
+{
+    //  Requires:   1.) ddx() and ddy() are present in the current Cg profile.
+    //              2.) The GPU driver is using fine/high-quality derivatives.
+    //              3.) quad_vector describes the current fragment's location in
+    //                  its 2x2 pixel quad using get_quad_vector()'s conventions.
+    //              4.) curr is any vector you wish to get neighboring values of.
+    //  Returns:    Values of an input vector (curr) at neighboring fragments
+    //              adjacent x, adjacent y, and diagonal (via out parameters).
+    adjx = curr - ddx(curr) * quad_vector.z;
+    adjy = curr - ddy(curr) * quad_vector.w;
+    diag = adjx - ddy(adjx) * quad_vector.w;
+}
+
+void quad_gather(const vec4 quad_vector, const vec3 curr,
+    out vec3 adjx, out vec3 adjy, out vec3 diag)
+{
+    //  vec3 version
+    adjx = curr - ddx(curr) * quad_vector.z;
+    adjy = curr - ddy(curr) * quad_vector.w;
+    diag = adjx - ddy(adjx) * quad_vector.w;
+}
+
+void quad_gather(const vec4 quad_vector, const vec2 curr,
+    out vec2 adjx, out vec2 adjy, out vec2 diag)
+{
+    //  vec2 version
+    adjx = curr - ddx(curr) * quad_vector.z;
+    adjy = curr - ddy(curr) * quad_vector.w;
+    diag = adjx - ddy(adjx) * quad_vector.w;
+}
+
+vec4 quad_gather(const vec4 quad_vector, const float curr)
+{
+    //  Float version:
+    //  Returns:    return.x == current
+    //              return.y == adjacent x
+    //              return.z == adjacent y
+    //              return.w == diagonal
+    vec4 all = vec4(curr);
+    all.y = all.x - ddx(all.x) * quad_vector.z;
+    all.zw = all.xy - ddy(all.xy) * quad_vector.w;
+    return all;
+}
+
+vec4 quad_gather_sum(const vec4 quad_vector, const vec4 curr)
+{
+    //  Requires:   Same as quad_gather()
+    //  Returns:    Sum of an input vector (curr) at all fragments in a quad.
+    vec4 adjx, adjy, diag;
+    quad_gather(quad_vector, curr, adjx, adjy, diag);
+    return (curr + adjx + adjy + diag);
+}
+
+vec3 quad_gather_sum(const vec4 quad_vector, const vec3 curr)
+{
+    //  vec3 version:
+    vec3 adjx, adjy, diag;
+    quad_gather(quad_vector, curr, adjx, adjy, diag);
+    return (curr + adjx + adjy + diag);
+}
+
+vec2 quad_gather_sum(const vec4 quad_vector, const vec2 curr)
+{
+    //  vec2 version:
+    vec2 adjx, adjy, diag;
+    quad_gather(quad_vector, curr, adjx, adjy, diag);
+    return (curr + adjx + adjy + diag);
+}
+
+float quad_gather_sum(const vec4 quad_vector, const float curr)
+{
+    //  Float version:
+    const vec4 all_values = quad_gather(quad_vector, curr);
+    return (all_values.x + all_values.y + all_values.z + all_values.w);
+}
+
+bool fine_derivatives_working(const vec4 quad_vector, vec4 curr)
+{
+    //  Requires:   1.) ddx() and ddy() are present in the current Cg profile.
+    //              2.) quad_vector describes the current fragment's location in
+    //                  its 2x2 pixel quad using get_quad_vector()'s conventions.
+    //              3.) curr must be a test vector with non-constant derivatives
+    //                  (its value should change nonlinearly across fragments).
+    //  Returns:    true if fine/hybrid/high-quality derivatives are used, or
+    //              false if coarse derivatives are used or inconclusive
+    //  Usage:      Test whether quad-pixel communication is working!
+    //  Method:     We can confirm fine derivatives are used if the following
+    //              holds (ever, for any value at any fragment):
+    //                  (ddy(curr) != ddy(adjx)) or (ddx(curr) != ddx(adjy))
+    //              The more values we test (e.g. test a vec4 two ways), the
+    //              easier it is to demonstrate fine derivatives are working.
+    //  TODO: Check for floating point exact comparison issues!
+    vec4 ddx_curr = ddx(curr);
+    vec4 ddy_curr = ddy(curr);
+    vec4 adjx = curr - ddx_curr * quad_vector.z;
+    vec4 adjy = curr - ddy_curr * quad_vector.w;
+    bool ddy_different = any(ddy_curr != ddy(adjx));
+    bool ddx_different = any(ddx_curr != ddx(adjy));
+    return any(bool2(ddy_different, ddx_different));
+}
+
+bool fine_derivatives_working_fast(const vec4 quad_vector, float curr)
+{
+    //  Requires:   Same as fine_derivatives_working()
+    //  Returns:    Same as fine_derivatives_working()
+    //  Usage:      This is faster than fine_derivatives_working() but more
+    //              likely to return false negatives, so it's less useful for
+    //              offline testing/debugging.  It's also useless as the basis
+    //              for dynamic runtime branching as of May 2014: Derivatives
+    //              (and quad-pixel communication) are currently disallowed in
+    //              branches.  However, future GPU's may allow you to use them
+    //              in dynamic branches if you promise the branch condition
+    //              evaluates the same for every fragment in the quad (and/or if
+    //              the driver enforces that promise by making a single fragment
+    //              control branch decisions).  If that ever happens, this
+    //              version may become a more economical choice.
+    float ddx_curr = ddx(curr);
+    float ddy_curr = ddy(curr);
+    float adjx = curr - ddx_curr * quad_vector.z;
+    return (ddy_curr != ddy(adjx));
+}
+
+#endif  //  QUAD_PIXEL_COMMUNICATION_H
+

include/special-functions-old.h

@@ -0,0 +1,498 @@
+#ifndef SPECIAL_FUNCTIONS_H
+#define SPECIAL_FUNCTIONS_H
+
+/////////////////////////////////  MIT LICENSE  ////////////////////////////////
+
+//  Copyright (C) 2014 TroggleMonkey
+//
+//  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.
+
+
+/////////////////////////////////  DESCRIPTION  ////////////////////////////////
+
+//  This file implements the following mathematical special functions:
+//  1.) erf() = 2/sqrt(pi) * indefinite_integral(e**(-x**2))
+//  2.) gamma(s), a real-numbered extension of the integer factorial function
+//  It also implements normalized_ligamma(s, z), a normalized lower incomplete
+//  gamma function for s < 0.5 only.  Both gamma() and normalized_ligamma() can
+//  be called with an _impl suffix to use an implementation version with a few
+//  extra precomputed parameters (which may be useful for the caller to reuse).
+//  See below for details.
+//
+//  Design Rationale:
+//  Pretty much every line of code in this file is duplicated four times for
+//  different input types (vec4/vec3/vec2/float).  This is unfortunate,
+//  but Cg doesn't allow function templates.  Macros would be far less verbose,
+//  but they would make the code harder to document and read.  I don't expect
+//  these functions will require a whole lot of maintenance changes unless
+//  someone ever has need for more robust incomplete gamma functions, so code
+//  duplication seems to be the lesser evil in this case.
+
+
+///////////////////////////  GAUSSIAN ERROR FUNCTION  //////////////////////////
+
+vec4 erf6(vec4 x)
+{
+    //  Requires:   x is the standard parameter to erf().
+    //  Returns:    Return an Abramowitz/Stegun approximation of erf(), where:
+    //                  erf(x) = 2/sqrt(pi) * integral(e**(-x**2))
+    //              This approximation has a max absolute error of 2.5*10**-5
+    //              with solid numerical robustness and efficiency.  See:
+	//                  https://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions
+	const vec4 one = vec4(1.0);
+	const vec4 sign_x = sign(x);
+	const vec4 t = one/(one + 0.47047*abs(x));
+	const vec4 result = one - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+vec3 erf6(const vec3 x)
+{
+    //  vec3 version:
+	const vec3 one = vec3(1.0);
+	const vec3 sign_x = sign(x);
+	const vec3 t = one/(one + 0.47047*abs(x));
+	const vec3 result = one - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+vec2 erf6(const vec2 x)
+{
+    //  vec2 version:
+	const vec2 one = vec2(1.0);
+	const vec2 sign_x = sign(x);
+	const vec2 t = one/(one + 0.47047*abs(x));
+	const vec2 result = one - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+float erf6(const float x)
+{
+    //  Float version:
+	const float sign_x = sign(x);
+	const float t = 1.0/(1.0 + 0.47047*abs(x));
+	const float result = 1.0 - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+vec4 erft(const vec4 x)
+{
+    //  Requires:   x is the standard parameter to erf().
+    //  Returns:    Approximate erf() with the hyperbolic tangent.  The error is
+    //              visually noticeable, but it's blazing fast and perceptually
+    //              close...at least on ATI hardware.  See:
+    //                  http://www.maplesoft.com/applications/view.aspx?SID=5525&view=html
+    //  Warning:    Only use this if your hardware drivers correctly implement
+    //              tanh(): My nVidia 8800GTS returns garbage output.
+	return tanh(1.202760580 * x);
+}
+
+vec3 erft(const vec3 x)
+{
+    //  vec3 version:
+	return tanh(1.202760580 * x);
+}
+
+vec2 erft(const vec2 x)
+{
+    //  vec2 version:
+	return tanh(1.202760580 * x);
+}
+
+float erft(const float x)
+{
+    //  Float version:
+	return tanh(1.202760580 * x);
+}
+
+vec4 erf(const vec4 x)
+{
+    //  Requires:   x is the standard parameter to erf().
+    //  Returns:    Some approximation of erf(x), depending on user settings.
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+vec3 erf(const vec3 x)
+{
+    //  vec3 version:
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+vec2 erf(const vec2 x)
+{
+    //  vec2 version:
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+float erf(const float x)
+{
+    //  Float version:
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+
+///////////////////////////  COMPLETE GAMMA FUNCTION  //////////////////////////
+
+vec4 gamma_impl(const vec4 s, const vec4 s_inv)
+{
+    //  Requires:   1.) s is the standard parameter to the gamma function, and
+    //                  it should lie in the [0, 36] range.
+    //              2.) s_inv = 1.0/s.  This implementation function requires
+    //                  the caller to precompute this value, giving users the
+    //                  opportunity to reuse it.
+    //  Returns:    Return approximate gamma function (real-numbered factorial)
+    //              output using the Lanczos approximation with two coefficients
+    //              calculated using Paul Godfrey's method here:
+    //                  http://my.fit.edu/~gabdo/gamma.txt
+    //              An optimal g value for s in [0, 36] is ~1.12906830989, with
+    //              a maximum relative error of 0.000463 for 2**16 equally
+    //              evals.  We could use three coeffs (0.0000346 error) without
+    //              hurting latency, but this allows more parallelism with
+    //              outside instructions.
+	const vec4 g = vec4(1.12906830989);
+	const vec4 c0 = vec4(0.8109119309638332633713423362694399653724431);
+	const vec4 c1 = vec4(0.4808354605142681877121661197951496120000040);
+	const vec4 e = vec4(2.71828182845904523536028747135266249775724709);
+	const vec4 sph = s + vec4(0.5);
+	const vec4 lanczos_sum = c0 + c1/(s + vec4(1.0));
+	const vec4 base = (sph + g)/e;  //  or (s + g + vec4(0.5))/e
+	//  gamma(s + 1) = base**sph * lanczos_sum; divide by s for gamma(s).
+	//  This has less error for small s's than (s -= 1.0) at the beginning.
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+vec3 gamma_impl(const vec3 s, const vec3 s_inv)
+{
+    //  vec3 version:
+	const vec3 g = vec3(1.12906830989);
+	const vec3 c0 = vec3(0.8109119309638332633713423362694399653724431);
+	const vec3 c1 = vec3(0.4808354605142681877121661197951496120000040);
+	const vec3 e = vec3(2.71828182845904523536028747135266249775724709);
+	const vec3 sph = s + vec3(0.5);
+	const vec3 lanczos_sum = c0 + c1/(s + vec3(1.0));
+	const vec3 base = (sph + g)/e;
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+vec2 gamma_impl(const vec2 s, const vec2 s_inv)
+{
+    //  vec2 version:
+	const vec2 g = vec2(1.12906830989);
+	const vec2 c0 = vec2(0.8109119309638332633713423362694399653724431);
+	const vec2 c1 = vec2(0.4808354605142681877121661197951496120000040);
+	const vec2 e = vec2(2.71828182845904523536028747135266249775724709);
+	const vec2 sph = s + vec2(0.5);
+	const vec2 lanczos_sum = c0 + c1/(s + vec2(1.0));
+	const vec2 base = (sph + g)/e;
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+float gamma_impl(const float s, const float s_inv)
+{
+    //  Float version:
+	const float g = 1.12906830989;
+	const float c0 = 0.8109119309638332633713423362694399653724431;
+	const float c1 = 0.4808354605142681877121661197951496120000040;
+	const float e = 2.71828182845904523536028747135266249775724709;
+	const float sph = s + 0.5;
+	const float lanczos_sum = c0 + c1/(s + 1.0);
+	const float base = (sph + g)/e;
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+vec4 gamma(const vec4 s)
+{
+    //  Requires:   s is the standard parameter to the gamma function, and it
+    //              should lie in the [0, 36] range.
+    //  Returns:    Return approximate gamma function output with a maximum
+    //              relative error of 0.000463.  See gamma_impl for details.
+	return gamma_impl(s, vec4(1.0)/s);
+}
+
+vec3 gamma(const vec3 s)
+{
+    //  vec3 version:
+	return gamma_impl(s, vec3(1.0)/s);
+}
+
+vec2 gamma(const vec2 s)
+{
+    //  vec2 version:
+	return gamma_impl(s, vec2(1.0)/s);
+}
+
+float gamma(const float s)
+{
+    //  Float version:
+	return gamma_impl(s, 1.0/s);
+}
+
+
+////////////////  INCOMPLETE GAMMA FUNCTIONS (RESTRICTED INPUT)  ///////////////
+
+//  Lower incomplete gamma function for small s and z (implementation):
+vec4 ligamma_small_z_impl(const vec4 s, const vec4 z, const vec4 s_inv)
+{
+    //  Requires:   1.) s < ~0.5
+    //              2.) z <= ~0.775075
+    //              3.) s_inv = 1.0/s (precomputed for outside reuse)
+    //  Returns:    A series representation for the lower incomplete gamma
+    //              function for small s and small z (4 terms).
+    //  The actual "rolled up" summation looks like:
+	//      last_sign = 1.0; last_pow = 1.0; last_factorial = 1.0;
+	//      sum = last_sign * last_pow / ((s + k) * last_factorial)
+	//      for(int i = 0; i < 4; ++i)
+	//      {
+	//          last_sign *= -1.0; last_pow *= z; last_factorial *= i;
+	//          sum += last_sign * last_pow / ((s + k) * last_factorial);
+	//      }
+	//  Unrolled, constant-unfolded and arranged for madds and parallelism:
+	const vec4 scale = pow(z, s);
+	vec4 sum = s_inv;  //  Summation iteration 0 result
+	//  Summation iterations 1, 2, and 3:
+	const vec4 z_sq = z*z;
+	const vec4 denom1 = s + vec4(1.0);
+	const vec4 denom2 = 2.0*s + vec4(4.0);
+	const vec4 denom3 = 6.0*s + vec4(18.0);
+	//vec4 denom4 = 24.0*s + vec4(96.0);
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	//sum += z_sq * z_sq / denom4;
+	//  Scale and return:
+	return scale * sum;
+}
+
+vec3 ligamma_small_z_impl(const vec3 s, const vec3 z, const vec3 s_inv)
+{
+    //  vec3 version:
+	const vec3 scale = pow(z, s);
+	vec3 sum = s_inv;
+	const vec3 z_sq = z*z;
+	const vec3 denom1 = s + vec3(1.0);
+	const vec3 denom2 = 2.0*s + vec3(4.0);
+	const vec3 denom3 = 6.0*s + vec3(18.0);
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	return scale * sum;
+}
+
+vec2 ligamma_small_z_impl(const vec2 s, const vec2 z, const vec2 s_inv)
+{
+    //  vec2 version:
+	const vec2 scale = pow(z, s);
+	vec2 sum = s_inv;
+	const vec2 z_sq = z*z;
+	const vec2 denom1 = s + vec2(1.0);
+	const vec2 denom2 = 2.0*s + vec2(4.0);
+	const vec2 denom3 = 6.0*s + vec2(18.0);
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	return scale * sum;
+}
+
+float ligamma_small_z_impl(const float s, const float z, const float s_inv)
+{
+    //  Float version:
+	const float scale = pow(z, s);
+	float sum = s_inv;
+	const float z_sq = z*z;
+	const float denom1 = s + 1.0;
+	const float denom2 = 2.0*s + 4.0;
+	const float denom3 = 6.0*s + 18.0;
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	return scale * sum;
+}
+
+//  Upper incomplete gamma function for small s and large z (implementation):
+vec4 uigamma_large_z_impl(const vec4 s, const vec4 z)
+{
+    //  Requires:   1.) s < ~0.5
+    //              2.) z > ~0.775075
+    //  Returns:    Gauss's continued fraction representation for the upper
+    //              incomplete gamma function (4 terms).
+	//  The "rolled up" continued fraction looks like this.  The denominator
+    //  is truncated, and it's calculated "from the bottom up:"
+	//      denom = vec4('inf');
+	//      vec4 one = vec4(1.0);
+	//      for(int i = 4; i > 0; --i)
+	//      {
+	//          denom = ((i * 2.0) - one) + z - s + (i * (s - i))/denom;
+	//      }
+	//  Unrolled and constant-unfolded for madds and parallelism:
+	const vec4 numerator = pow(z, s) * exp(-z);
+	vec4 denom = vec4(7.0) + z - s;
+	denom = vec4(5.0) + z - s + (3.0*s - vec4(9.0))/denom;
+	denom = vec4(3.0) + z - s + (2.0*s - vec4(4.0))/denom;
+	denom = vec4(1.0) + z - s + (s - vec4(1.0))/denom;
+	return numerator / denom;
+}
+
+vec3 uigamma_large_z_impl(const vec3 s, const vec3 z)
+{
+    //  vec3 version:
+	const vec3 numerator = pow(z, s) * exp(-z);
+	vec3 denom = vec3(7.0) + z - s;
+	denom = vec3(5.0) + z - s + (3.0*s - vec3(9.0))/denom;
+	denom = vec3(3.0) + z - s + (2.0*s - vec3(4.0))/denom;
+	denom = vec3(1.0) + z - s + (s - vec3(1.0))/denom;
+	return numerator / denom;
+}
+
+vec2 uigamma_large_z_impl(const vec2 s, const vec2 z)
+{
+    //  vec2 version:
+	const vec2 numerator = pow(z, s) * exp(-z);
+	vec2 denom = vec2(7.0) + z - s;
+	denom = vec2(5.0) + z - s + (3.0*s - vec2(9.0))/denom;
+	denom = vec2(3.0) + z - s + (2.0*s - vec2(4.0))/denom;
+	denom = vec2(1.0) + z - s + (s - vec2(1.0))/denom;
+	return numerator / denom;
+}
+
+float uigamma_large_z_impl(const float s, const float z)
+{
+    //  Float version:
+	const float numerator = pow(z, s) * exp(-z);
+	float denom = 7.0 + z - s;
+	denom = 5.0 + z - s + (3.0*s - 9.0)/denom;
+	denom = 3.0 + z - s + (2.0*s - 4.0)/denom;
+	denom = 1.0 + z - s + (s - 1.0)/denom;
+	return numerator / denom;
+}
+
+//  Normalized lower incomplete gamma function for small s (implementation):
+vec4 normalized_ligamma_impl(const vec4 s, const vec4 z,
+    const vec4 s_inv, const vec4 gamma_s_inv)
+{
+    //  Requires:   1.) s < ~0.5
+    //              2.) s_inv = 1/s (precomputed for outside reuse)
+    //              3.) gamma_s_inv = 1/gamma(s) (precomputed for outside reuse)
+    //  Returns:    Approximate the normalized lower incomplete gamma function
+    //              for s < 0.5.  Since we only care about s < 0.5, we only need
+    //              to evaluate two branches (not four) based on z.  Each branch
+    //              uses four terms, with a max relative error of ~0.00182.  The
+    //              branch threshold and specifics were adapted for fewer terms
+    //              from Gil/Segura/Temme's paper here:
+    //                  http://oai.cwi.nl/oai/asset/20433/20433B.pdf
+	//  Evaluate both branches: Real branches test slower even when available.
+	const vec4 thresh = vec4(0.775075);
+	const bool4 z_is_large = z > thresh;
+	const vec4 large_z = vec4(1.0) - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const vec4 small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	//  Combine the results from both branches:
+	return large_z * vec4(z_is_large) + small_z * vec4(!z_is_large);
+}
+
+vec3 normalized_ligamma_impl(const vec3 s, const vec3 z,
+    const vec3 s_inv, const vec3 gamma_s_inv)
+{
+    //  vec3 version:
+	const vec3 thresh = vec3(0.775075);
+	const bool3 z_is_large = z > thresh;
+	const vec3 large_z = vec3(1.0) - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const vec3 small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	return large_z * vec3(z_is_large) + small_z * vec3(!z_is_large);
+}
+
+vec2 normalized_ligamma_impl(const vec2 s, const vec2 z,
+    const vec2 s_inv, const vec2 gamma_s_inv)
+{
+    //  vec2 version:
+	const vec2 thresh = vec2(0.775075);
+	const bool2 z_is_large = z > thresh;
+	const vec2 large_z = vec2(1.0) - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const vec2 small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	return large_z * vec2(z_is_large) + small_z * vec2(!z_is_large);
+}
+
+float normalized_ligamma_impl(const float s, const float z,
+    const float s_inv, const float gamma_s_inv)
+{
+    //  Float version:
+	const float thresh = 0.775075;
+	const bool z_is_large = z > thresh;
+	const float large_z = 1.0 - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const float small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	return large_z * float(z_is_large) + small_z * float(!z_is_large);
+}
+
+//  Normalized lower incomplete gamma function for small s:
+vec4 normalized_ligamma(const vec4 s, const vec4 z)
+{
+    //  Requires:   s < ~0.5
+    //  Returns:    Approximate the normalized lower incomplete gamma function
+    //              for s < 0.5.  See normalized_ligamma_impl() for details.
+	const vec4 s_inv = vec4(1.0)/s;
+	const vec4 gamma_s_inv = vec4(1.0)/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}
+
+vec3 normalized_ligamma(const vec3 s, const vec3 z)
+{
+    //  vec3 version:
+	const vec3 s_inv = vec3(1.0)/s;
+	const vec3 gamma_s_inv = vec3(1.0)/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}
+
+vec2 normalized_ligamma(const vec2 s, const vec2 z)
+{
+    //  vec2 version:
+	const vec2 s_inv = vec2(1.0)/s;
+	const vec2 gamma_s_inv = vec2(1.0)/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}
+
+float normalized_ligamma(const float s, const float z)
+{
+    //  Float version:
+	const float s_inv = 1.0/s;
+	const float gamma_s_inv = 1.0/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}
+
+
+#endif  //  SPECIAL_FUNCTIONS_H
+
+

include/special-functions.h

@@ -0,0 +1,492 @@
+/////////////////////////////////  MIT LICENSE  ////////////////////////////////
+
+//  Copyright (C) 2014 TroggleMonkey
+//
+//  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.
+
+
+/////////////////////////////////  DESCRIPTION  ////////////////////////////////
+
+//  This file implements the following mathematical special functions:
+//  1.) erf() = 2/sqrt(pi) * indefinite_integral(e**(-x**2))
+//  2.) gamma(s), a real-numbered extension of the integer factorial function
+//  It also implements normalized_ligamma(s, z), a normalized lower incomplete
+//  gamma function for s < 0.5 only.  Both gamma() and normalized_ligamma() can
+//  be called with an _impl suffix to use an implementation version with a few
+//  extra precomputed parameters (which may be useful for the caller to reuse).
+//  See below for details.
+//
+//  Design Rationale:
+//  Pretty much every line of code in this file is duplicated four times for
+//  different input types (vec4/vec3/vec2/float).  This is unfortunate,
+//  but Cg doesn't allow function templates.  Macros would be far less verbose,
+//  but they would make the code harder to document and read.  I don't expect
+//  these functions will require a whole lot of maintenance changes unless
+//  someone ever has need for more robust incomplete gamma functions, so code
+//  duplication seems to be the lesser evil in this case.
+
+
+///////////////////////////  GAUSSIAN ERROR FUNCTION  //////////////////////////
+
+vec4 erf6(vec4 x)
+{
+    //  Requires:   x is the standard parameter to erf().
+    //  Returns:    Return an Abramowitz/Stegun approximation of erf(), where:
+    //                  erf(x) = 2/sqrt(pi) * integral(e**(-x**2))
+    //              This approximation has a max absolute error of 2.5*10**-5
+    //              with solid numerical robustness and efficiency.  See:
+	//                  https://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions
+	const vec4 one = vec4(1.0);
+	const vec4 sign_x = sign(x);
+	const vec4 t = one/(one + 0.47047*abs(x));
+	const vec4 result = one - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+vec3 erf6(const vec3 x)
+{
+    //  vec3 version:
+	const vec3 one = vec3(1.0);
+	const vec3 sign_x = sign(x);
+	const vec3 t = one/(one + 0.47047*abs(x));
+	const vec3 result = one - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+vec2 erf6(const vec2 x)
+{
+    //  vec2 version:
+	const vec2 one = vec2(1.0);
+	const vec2 sign_x = sign(x);
+	const vec2 t = one/(one + 0.47047*abs(x));
+	const vec2 result = one - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+float erf6(const float x)
+{
+    //  Float version:
+	const float sign_x = sign(x);
+	const float t = 1.0/(1.0 + 0.47047*abs(x));
+	const float result = 1.0 - t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*
+		exp(-(x*x));
+	return result * sign_x;
+}
+
+vec4 erft(const vec4 x)
+{
+    //  Requires:   x is the standard parameter to erf().
+    //  Returns:    Approximate erf() with the hyperbolic tangent.  The error is
+    //              visually noticeable, but it's blazing fast and perceptually
+    //              close...at least on ATI hardware.  See:
+    //                  http://www.maplesoft.com/applications/view.aspx?SID=5525&view=html
+    //  Warning:    Only use this if your hardware drivers correctly implement
+    //              tanh(): My nVidia 8800GTS returns garbage output.
+	return tanh(1.202760580 * x);
+}
+
+vec3 erft(const vec3 x)
+{
+    //  vec3 version:
+	return tanh(1.202760580 * x);
+}
+
+vec2 erft(const vec2 x)
+{
+    //  vec2 version:
+	return tanh(1.202760580 * x);
+}
+
+float erft(const float x)
+{
+    //  Float version:
+	return tanh(1.202760580 * x);
+}
+
+vec4 erf(const vec4 x)
+{
+    //  Requires:   x is the standard parameter to erf().
+    //  Returns:    Some approximation of erf(x), depending on user settings.
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+vec3 erf(const vec3 x)
+{
+    //  vec3 version:
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+vec2 erf(const vec2 x)
+{
+    //  vec2 version:
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+float erf(const float x)
+{
+    //  Float version:
+	#ifdef ERF_FAST_APPROXIMATION
+		return erft(x);
+	#else
+		return erf6(x);
+	#endif
+}
+
+///////////////////////////  COMPLETE GAMMA FUNCTION  //////////////////////////
+
+vec4 gamma_impl(const vec4 s, const vec4 s_inv)
+{
+    //  Requires:   1.) s is the standard parameter to the gamma function, and
+    //                  it should lie in the [0, 36] range.
+    //              2.) s_inv = 1.0/s.  This implementation function requires
+    //                  the caller to precompute this value, giving users the
+    //                  opportunity to reuse it.
+    //  Returns:    Return approximate gamma function (real-numbered factorial)
+    //              output using the Lanczos approximation with two coefficients
+    //              calculated using Paul Godfrey's method here:
+    //                  http://my.fit.edu/~gabdo/gamma.txt
+    //              An optimal g value for s in [0, 36] is ~1.12906830989, with
+    //              a maximum relative error of 0.000463 for 2**16 equally
+    //              evals.  We could use three coeffs (0.0000346 error) without
+    //              hurting latency, but this allows more parallelism with
+    //              outside instructions.
+	const vec4 g = vec4(1.12906830989);
+	const vec4 c0 = vec4(0.8109119309638332633713423362694399653724431);
+	const vec4 c1 = vec4(0.4808354605142681877121661197951496120000040);
+	const vec4 e = vec4(2.71828182845904523536028747135266249775724709);
+	const vec4 sph = s + vec4(0.5);
+	const vec4 lanczos_sum = c0 + c1/(s + vec4(1.0));
+	const vec4 base = (sph + g)/e;  //  or (s + g + vec4(0.5))/e
+	//  gamma(s + 1) = base**sph * lanczos_sum; divide by s for gamma(s).
+	//  This has less error for small s's than (s -= 1.0) at the beginning.
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+vec3 gamma_impl(const vec3 s, const vec3 s_inv)
+{
+    //  vec3 version:
+	const vec3 g = vec3(1.12906830989);
+	const vec3 c0 = vec3(0.8109119309638332633713423362694399653724431);
+	const vec3 c1 = vec3(0.4808354605142681877121661197951496120000040);
+	const vec3 e = vec3(2.71828182845904523536028747135266249775724709);
+	const vec3 sph = s + vec3(0.5);
+	const vec3 lanczos_sum = c0 + c1/(s + vec3(1.0));
+	const vec3 base = (sph + g)/e;
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+vec2 gamma_impl(const vec2 s, const vec2 s_inv)
+{
+    //  vec2 version:
+	const vec2 g = vec2(1.12906830989);
+	const vec2 c0 = vec2(0.8109119309638332633713423362694399653724431);
+	const vec2 c1 = vec2(0.4808354605142681877121661197951496120000040);
+	const vec2 e = vec2(2.71828182845904523536028747135266249775724709);
+	const vec2 sph = s + vec2(0.5);
+	const vec2 lanczos_sum = c0 + c1/(s + vec2(1.0));
+	const vec2 base = (sph + g)/e;
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+float gamma_impl(const float s, const float s_inv)
+{
+    //  Float version:
+	const float g = 1.12906830989;
+	const float c0 = 0.8109119309638332633713423362694399653724431;
+	const float c1 = 0.4808354605142681877121661197951496120000040;
+	const float e = 2.71828182845904523536028747135266249775724709;
+	const float sph = s + 0.5;
+	const float lanczos_sum = c0 + c1/(s + 1.0);
+	const float base = (sph + g)/e;
+	return (pow(base, sph) * lanczos_sum) * s_inv;
+}
+
+vec4 gamma(const vec4 s)
+{
+    //  Requires:   s is the standard parameter to the gamma function, and it
+    //              should lie in the [0, 36] range.
+    //  Returns:    Return approximate gamma function output with a maximum
+    //              relative error of 0.000463.  See gamma_impl for details.
+	return gamma_impl(s, vec4(1.0)/s);
+}
+
+vec3 gamma(const vec3 s)
+{
+    //  vec3 version:
+	return gamma_impl(s, vec3(1.0)/s);
+}
+
+vec2 gamma(const vec2 s)
+{
+    //  vec2 version:
+	return gamma_impl(s, vec2(1.0)/s);
+}
+
+float gamma(const float s)
+{
+    //  Float version:
+	return gamma_impl(s, 1.0/s);
+}
+
+////////////////  INCOMPLETE GAMMA FUNCTIONS (RESTRICTED INPUT)  ///////////////
+
+//  Lower incomplete gamma function for small s and z (implementation):
+vec4 ligamma_small_z_impl(const vec4 s, const vec4 z, const vec4 s_inv)
+{
+    //  Requires:   1.) s < ~0.5
+    //              2.) z <= ~0.775075
+    //              3.) s_inv = 1.0/s (precomputed for outside reuse)
+    //  Returns:    A series representation for the lower incomplete gamma
+    //              function for small s and small z (4 terms).
+    //  The actual "rolled up" summation looks like:
+	//      last_sign = 1.0; last_pow = 1.0; last_factorial = 1.0;
+	//      sum = last_sign * last_pow / ((s + k) * last_factorial)
+	//      for(int i = 0; i < 4; ++i)
+	//      {
+	//          last_sign *= -1.0; last_pow *= z; last_factorial *= i;
+	//          sum += last_sign * last_pow / ((s + k) * last_factorial);
+	//      }
+	//  Unrolled, constant-unfolded and arranged for madds and parallelism:
+	const vec4 scale = pow(z, s);
+	vec4 sum = s_inv;  //  Summation iteration 0 result
+	//  Summation iterations 1, 2, and 3:
+	const vec4 z_sq = z*z;
+	const vec4 denom1 = s + vec4(1.0);
+	const vec4 denom2 = 2.0*s + vec4(4.0);
+	const vec4 denom3 = 6.0*s + vec4(18.0);
+	//vec4 denom4 = 24.0*s + vec4(96.0);
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	//sum += z_sq * z_sq / denom4;
+	//  Scale and return:
+	return scale * sum;
+}
+
+vec3 ligamma_small_z_impl(const vec3 s, const vec3 z, const vec3 s_inv)
+{
+    //  vec3 version:
+	const vec3 scale = pow(z, s);
+	vec3 sum = s_inv;
+	const vec3 z_sq = z*z;
+	const vec3 denom1 = s + vec3(1.0);
+	const vec3 denom2 = 2.0*s + vec3(4.0);
+	const vec3 denom3 = 6.0*s + vec3(18.0);
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	return scale * sum;
+}
+
+vec2 ligamma_small_z_impl(const vec2 s, const vec2 z, const vec2 s_inv)
+{
+    //  vec2 version:
+	const vec2 scale = pow(z, s);
+	vec2 sum = s_inv;
+	const vec2 z_sq = z*z;
+	const vec2 denom1 = s + vec2(1.0);
+	const vec2 denom2 = 2.0*s + vec2(4.0);
+	const vec2 denom3 = 6.0*s + vec2(18.0);
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	return scale * sum;
+}
+
+float ligamma_small_z_impl(const float s, const float z, const float s_inv)
+{
+    //  Float version:
+	const float scale = pow(z, s);
+	float sum = s_inv;
+	const float z_sq = z*z;
+	const float denom1 = s + 1.0;
+	const float denom2 = 2.0*s + 4.0;
+	const float denom3 = 6.0*s + 18.0;
+	sum -= z/denom1;
+	sum += z_sq/denom2;
+	sum -= z * z_sq/denom3;
+	return scale * sum;
+}
+
+//  Upper incomplete gamma function for small s and large z (implementation):
+vec4 uigamma_large_z_impl(const vec4 s, const vec4 z)
+{
+    //  Requires:   1.) s < ~0.5
+    //              2.) z > ~0.775075
+    //  Returns:    Gauss's continued fraction representation for the upper
+    //              incomplete gamma function (4 terms).
+	//  The "rolled up" continued fraction looks like this.  The denominator
+    //  is truncated, and it's calculated "from the bottom up:"
+	//      denom = vec4('inf');
+	//      vec4 one = vec4(1.0);
+	//      for(int i = 4; i > 0; --i)
+	//      {
+	//          denom = ((i * 2.0) - one) + z - s + (i * (s - i))/denom;
+	//      }
+	//  Unrolled and constant-unfolded for madds and parallelism:
+	const vec4 numerator = pow(z, s) * exp(-z);
+	vec4 denom = vec4(7.0) + z - s;
+	denom = vec4(5.0) + z - s + (3.0*s - vec4(9.0))/denom;
+	denom = vec4(3.0) + z - s + (2.0*s - vec4(4.0))/denom;
+	denom = vec4(1.0) + z - s + (s - vec4(1.0))/denom;
+	return numerator / denom;
+}
+
+vec3 uigamma_large_z_impl(const vec3 s, const vec3 z)
+{
+    //  vec3 version:
+	const vec3 numerator = pow(z, s) * exp(-z);
+	vec3 denom = vec3(7.0) + z - s;
+	denom = vec3(5.0) + z - s + (3.0*s - vec3(9.0))/denom;
+	denom = vec3(3.0) + z - s + (2.0*s - vec3(4.0))/denom;
+	denom = vec3(1.0) + z - s + (s - vec3(1.0))/denom;
+	return numerator / denom;
+}
+
+vec2 uigamma_large_z_impl(const vec2 s, const vec2 z)
+{
+    //  vec2 version:
+	const vec2 numerator = pow(z, s) * exp(-z);
+	vec2 denom = vec2(7.0) + z - s;
+	denom = vec2(5.0) + z - s + (3.0*s - vec2(9.0))/denom;
+	denom = vec2(3.0) + z - s + (2.0*s - vec2(4.0))/denom;
+	denom = vec2(1.0) + z - s + (s - vec2(1.0))/denom;
+	return numerator / denom;
+}
+
+float uigamma_large_z_impl(const float s, const float z)
+{
+    //  Float version:
+	const float numerator = pow(z, s) * exp(-z);
+	float denom = 7.0 + z - s;
+	denom = 5.0 + z - s + (3.0*s - 9.0)/denom;
+	denom = 3.0 + z - s + (2.0*s - 4.0)/denom;
+	denom = 1.0 + z - s + (s - 1.0)/denom;
+	return numerator / denom;
+}
+
+//  Normalized lower incomplete gamma function for small s (implementation):
+vec4 normalized_ligamma_impl(const vec4 s, const vec4 z,
+    const vec4 s_inv, const vec4 gamma_s_inv)
+{
+    //  Requires:   1.) s < ~0.5
+    //              2.) s_inv = 1/s (precomputed for outside reuse)
+    //              3.) gamma_s_inv = 1/gamma(s) (precomputed for outside reuse)
+    //  Returns:    Approximate the normalized lower incomplete gamma function
+    //              for s < 0.5.  Since we only care about s < 0.5, we only need
+    //              to evaluate two branches (not four) based on z.  Each branch
+    //              uses four terms, with a max relative error of ~0.00182.  The
+    //              branch threshold and specifics were adapted for fewer terms
+    //              from Gil/Segura/Temme's paper here:
+    //                  http://oai.cwi.nl/oai/asset/20433/20433B.pdf
+	//  Evaluate both branches: Real branches test slower even when available.
+	const vec4 thresh = vec4(0.775075);
+	bvec4 z_is_large = greaterThan(z , thresh);
+	vec4 z_size_check = vec4(z_is_large.x ? 1.0 : 0.0, z_is_large.y ? 1.0 : 0.0, z_is_large.z ? 1.0 : 0.0, z_is_large.w ? 1.0 : 0.0);
+	const vec4 large_z = vec4(1.0) - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const vec4 small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	//  Combine the results from both branches:
+	return large_z * vec4(z_size_check) + small_z * vec4(z_size_check);
+}
+
+vec3 normalized_ligamma_impl(const vec3 s, const vec3 z,
+    const vec3 s_inv, const vec3 gamma_s_inv)
+{
+    //  vec3 version:
+	const vec3 thresh = vec3(0.775075);
+	bvec3 z_is_large = greaterThan(z , thresh);
+	vec3 z_size_check = vec3(z_is_large.x ? 1.0 : 0.0, z_is_large.y ? 1.0 : 0.0, z_is_large.z ? 1.0 : 0.0);
+	const vec3 large_z = vec3(1.0) - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const vec3 small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	return large_z * vec3(z_size_check) + small_z * vec3(z_size_check);
+}
+
+vec2 normalized_ligamma_impl(const vec2 s, const vec2 z,
+    const vec2 s_inv, const vec2 gamma_s_inv)
+{
+    //  vec2 version:
+	const vec2 thresh = vec2(0.775075);
+	bvec2 z_is_large = greaterThan(z , thresh);
+	vec2 z_size_check = vec2(z_is_large.x ? 1.0 : 0.0, z_is_large.y ? 1.0 : 0.0);
+	const vec2 large_z = vec2(1.0) - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const vec2 small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	return large_z * vec2(z_size_check) + small_z * vec2(z_size_check);
+}
+
+float normalized_ligamma_impl(const float s, const float z,
+    const float s_inv, const float gamma_s_inv)
+{
+    //  Float version:
+	const float thresh = 0.775075;
+	float z_size_check = 0.0;
+	if (z > thresh) z_size_check = 1.0;
+	const float large_z = 1.0 - uigamma_large_z_impl(s, z) * gamma_s_inv;
+	const float small_z = ligamma_small_z_impl(s, z, s_inv) * gamma_s_inv;
+	return large_z * float(z_size_check) + small_z * float(z_size_check);
+}
+
+//  Normalized lower incomplete gamma function for small s:
+vec4 normalized_ligamma(const vec4 s, const vec4 z)
+{
+    //  Requires:   s < ~0.5
+    //  Returns:    Approximate the normalized lower incomplete gamma function
+    //              for s < 0.5.  See normalized_ligamma_impl() for details.
+	const vec4 s_inv = vec4(1.0)/s;
+	const vec4 gamma_s_inv = vec4(1.0)/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}
+
+vec3 normalized_ligamma(const vec3 s, const vec3 z)
+{
+    //  vec3 version:
+	const vec3 s_inv = vec3(1.0)/s;
+	const vec3 gamma_s_inv = vec3(1.0)/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}
+
+vec2 normalized_ligamma(const vec2 s, const vec2 z)
+{
+    //  vec2 version:
+	const vec2 s_inv = vec2(1.0)/s;
+	const vec2 gamma_s_inv = vec2(1.0)/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}
+
+float normalized_ligamma(const float s, const float z)
+{
+    //  Float version:
+	const float s_inv = 1.0/s;
+	const float gamma_s_inv = 1.0/gamma_impl(s, s_inv);
+	return normalized_ligamma_impl(s, z, s_inv, gamma_s_inv);
+}