// Adaptive Sharpen
// Based on: https://github.com/libretro/slang-shaders/blob/master/sharpen/shaders/adaptive-sharpen.slang
//
// Adaptive sharpen - version 2015-05-15 - (requires ps >= 3.0)
// Designed for use post-resize, expects full range gamma light

//!BGFX EFFECT
//!VERSION 1
//!NAME Adaptive Sharpen
//!CATEGORY Sharpening
//!DESCRIPTION Content-aware sharpening algorithm that adapts to local image features. Preserves smooth gradients while enhancing edges with minimal halo artifacts.


//!PARAMETER
//!LABEL Sharpness
//!DESC Main sharpening strength control. Values between 0.3 and 2.0 work best for most content.
//!DEFAULT 0.8
//!MIN 0
//!MAX 2
//!STEP 0.01
float curveHeight;

//!TEXTURE
Texture2D INPUT;

//!TEXTURE
//!WIDTH INPUT_WIDTH
//!HEIGHT INPUT_HEIGHT
Texture2D OUTPUT;

//!SAMPLER
//!FILTER POINT
SamplerState sam;


//!PASS 1
//!IN INPUT
//!OUT OUTPUT
//!BLOCK_SIZE 16
//!NUM_THREADS 64

// Algorithm constants - these are optimized values, do not modify unless you understand the algorithm
#define curveslope      (curveHeight*1.5f)   // Sharpening curve slope for edge regions
#define D_overshoot     0.016f                // Max dark overshoot threshold before compression
#define D_comp_ratio    0.250f                // Max compression ratio for dark overshoot (1/0.25=4x)
#define L_overshoot     0.004f                // Max light overshoot threshold before compression
#define L_comp_ratio    0.167f                // Max compression ratio for light overshoot (1/0.167=6x)
#define max_scale_lim   10.0f                 // Absolute change limit before max compression (1/10=10%)

// Fast approximate gamma-aware luminance calculation
float ComputeLuminance(float3 RGB) {
	return sqrt((1.0f / 3.0f) * ((RGB * RGB).r + (RGB * RGB).g + (RGB * RGB).b));
}


void Pass1(uint2 blockStart, uint3 threadId) {
	uint2 gxy = (TileSwizzle8x8(threadId.x) << 1) + blockStart;

	const uint2 outputSize = GetOutputSize();
	if (gxy.x >= outputSize.x || gxy.y >= outputSize.y) {
		return;
	}

	float2 inputPt = GetInputPt();
	int i, j;

	// Load 8x8 source texel neighborhood with luminance pre-computed
	float4 src[8][8];
	[unroll]
	for (i = 0; i <= 6; i += 2) {
		[unroll]
		for (j = 0; j <= 6; j += 2) {
			// Skip corner texels that aren't needed (16 texels total)
			if ((i == 0 && j == 0) || (i == 6 && j == 0) || (i == 0 && j == 6) || (i == 6 && j == 6)) {
				continue;
			}

			float2 tpos = ((int2)gxy + int2(i, j) - 2) * inputPt;
			const float4 sr = INPUT.GatherRed(sam, tpos);
			const float4 sg = INPUT.GatherGreen(sam, tpos);
			const float4 sb = INPUT.GatherBlue(sam, tpos);

			// GatherRed returns:
			// w z
			// x y
			src[i][j].rgb = float3(sr.w, sg.w, sb.w);
			src[i][j].w = ComputeLuminance(src[i][j].rgb);
			src[i][j + 1].rgb = float3(sr.x, sg.x, sb.x);
			src[i][j + 1].w = ComputeLuminance(src[i][j + 1].rgb);
			src[i + 1][j].rgb = float3(sr.z, sg.z, sb.z);
			src[i + 1][j].w = ComputeLuminance(src[i + 1][j].rgb);
			src[i + 1][j + 1].rgb = float3(sr.y, sg.y, sb.y);
			src[i + 1][j + 1].w = ComputeLuminance(src[i + 1][j + 1].rgb);
		}
	}

	// Process 2x2 output block
	[unroll]
	for (i = 0; i <= 1; ++i) {
		[unroll]
		for (j = 0; j <= 1; ++j) {
			const uint2 destPos = gxy + uint2(i, j);

			float2 pos = (destPos + 0.5f) * inputPt;

			// Sample pattern (saturate out of range values BTB & WTW):
			// [                c22               ]
			// [           c24, c9,  c23          ]
			// [      c21, c1,  c2,  c3, c18      ]
			// [ c19, c10, c4,  c0,  c5, c11, c16 ]
			// [      c20, c6,  c7,  c8, c17      ]
			// [           c15, c12, c14          ]
			// [                c13               ]

			// Compute 3x3 Gaussian blur
			float3 blur = (2 * (src[i + 3][j + 2].rgb + src[i + 2][j + 3].rgb + src[i + 4][j + 3].rgb + src[i + 3][j + 4].rgb) +
			               (src[i + 2][j + 2].rgb + src[i + 4][j + 2].rgb + src[i + 2][j + 4].rgb + src[i + 4][j + 4].rgb) +
			               4 * src[i + 3][j + 3].rgb) / 16;
			float blur_Y = (blur.r * (1.0 / 3.0) + blur.g * (1.0 / 3.0) + blur.b * (1.0 / 3.0));

			// Edge detection using weighted matrix:
			// [           1          ]
			// [       4,  4,  4      ]
			// [   1,  4,  4,  4,  1  ]
			// [       4,  4,  4      ]
			// [           1          ]
			float edge = length(abs(blur - src[i + 3][j + 3].rgb) + abs(blur - src[i + 2][j + 2].rgb) +
			                    abs(blur - src[i + 3][j + 2].rgb) + abs(blur - src[i + 4][j + 2].rgb) +
			                    abs(blur - src[i + 2][j + 3].rgb) + abs(blur - src[i + 4][j + 3].rgb) +
			                    abs(blur - src[i + 2][j + 4].rgb) + abs(blur - src[i + 3][j + 4].rgb) +
			                    abs(blur - src[i + 4][j + 4].rgb) +
			                    0.25 * (abs(blur - src[i + 3][j + 1].rgb) + abs(blur - src[i + 1][j + 3].rgb) +
			                            abs(blur - src[i + 5][j + 3].rgb) + abs(blur - src[i + 3][j + 5].rgb))) * (1.0 / 3.0);

			// Edge detection contrast compression, center = 0.5
			edge *= min((0.8 + 2.7 * pow(2, (-7.4 * blur_Y))), 3.2);

			// Get center pixel luminance
			float c0_Y = src[i + 3][j + 3].w;

			// Build kernel array for local statistics
			float kernel[25] = {
				c0_Y, src[i + 2][j + 2].w, src[i + 3][j + 2].w, src[i + 4][j + 2].w,
				src[i + 2][j + 3].w, src[i + 4][j + 3].w, src[i + 2][j + 4].w, src[i + 3][j + 4].w, src[i + 4][j + 4].w,
				src[i + 3][j + 1].w, src[i + 1][j + 3].w, src[i + 5][j + 3].w, src[i + 3][j + 5].w, src[i + 3][j + 6].w,
				src[i + 4][j + 5].w, src[i + 2][j + 5].w, src[i + 6][j + 3].w, src[i + 5][j + 4].w, src[i + 5][j + 2].w,
				src[i][j + 3].w, src[i + 1][j + 4].w, src[i + 1][j + 2].w, src[i + 3][j].w, src[i + 4][j + 1].w, src[i + 2][j + 1].w
			};

			// Partial Laplacian outer pixel weighting scheme
			float mdiff_c0 = 0.03 + 4 * (abs(kernel[0] - kernel[2]) + abs(kernel[0] - kernel[4]) +
			                              abs(kernel[0] - kernel[5]) + abs(kernel[0] - kernel[7]) +
			                              0.25 * (abs(kernel[0] - kernel[1]) + abs(kernel[0] - kernel[3]) +
			                                      abs(kernel[0] - kernel[6]) + abs(kernel[0] - kernel[8])));

			float mdiff_c9 = (abs(kernel[9] - kernel[2]) + abs(kernel[9] - kernel[24]) +
			                  abs(kernel[9] - kernel[23]) + abs(kernel[9] - kernel[22]) +
			                  0.5 * (abs(kernel[9] - kernel[1]) + abs(kernel[9] - kernel[3])));

			float mdiff_c10 = (abs(kernel[10] - kernel[20]) + abs(kernel[10] - kernel[19]) +
			                   abs(kernel[10] - kernel[21]) + abs(kernel[10] - kernel[4]) +
			                   0.5 * (abs(kernel[10] - kernel[1]) + abs(kernel[10] - kernel[6])));

			float mdiff_c11 = (abs(kernel[11] - kernel[17]) + abs(kernel[11] - kernel[5]) +
			                   abs(kernel[11] - kernel[18]) + abs(kernel[11] - kernel[16]) +
			                   0.5 * (abs(kernel[11] - kernel[3]) + abs(kernel[11] - kernel[8])));

			float mdiff_c12 = (abs(kernel[12] - kernel[13]) + abs(kernel[12] - kernel[15]) +
			                   abs(kernel[12] - kernel[7]) + abs(kernel[12] - kernel[14]) +
			                   0.5 * (abs(kernel[12] - kernel[6]) + abs(kernel[12] - kernel[8])));

			float4 weights = float4(
				min((mdiff_c0 / mdiff_c9), 2.0),
				min((mdiff_c0 / mdiff_c10), 2.0),
				min((mdiff_c0 / mdiff_c11), 2.0),
				min((mdiff_c0 / mdiff_c12), 2.0)
			);

			// Negative Laplacian matrix with variable weights:
			// [          8*         ]
			// [      4,  1,  4      ]
			// [  8*, 1,      1,  8* ]
			// [      4,  1,  4      ]
			// [          8*         ]
			float neg_laplace = (0.25 * (kernel[2] + kernel[4] + kernel[5] + kernel[7]) +
			                     (kernel[1] + kernel[3] + kernel[6] + kernel[8]) +
			                     ((kernel[9] * weights.x) + (kernel[10] * weights.y) +
			                      (kernel[11] * weights.z) + (kernel[12] * weights.w))) /
			                    (5 + weights.x + weights.y + weights.z + weights.w);

			// Compute sharpening magnitude based on edge strength and Laplacian
			float sharpen_val = 0.01 + (curveHeight / (curveslope * pow(edge, 3.5) + 0.5)) -
			                    (curveHeight / (8192 * pow((edge * 2.2), 4.5) + 0.5));

			// Calculate sharpening difference and scale
			float sharpdiff = (c0_Y - neg_laplace) * (sharpen_val * 0.8);

			// Calculate local min/max using partial cocktail sort (no branching)
			[unroll]
			for (int k = 0; k < 2; ++k) {
				[unroll]
				for (int i1 = 1 + k; i1 < 25 - k; ++i1) {
					float temp = kernel[i1 - 1];
					kernel[i1 - 1] = min(kernel[i1 - 1], kernel[i1]);
					kernel[i1] = max(temp, kernel[i1]);
				}

				[unroll]
				for (int i2 = 23 - k; i2 > k; --i2) {
					float temp = kernel[i2 - 1];
					kernel[i2 - 1] = min(kernel[i2 - 1], kernel[i2]);
					kernel[i2] = max(temp, kernel[i2]);
				}
			}

			float nmax = max(((kernel[23] + kernel[24]) / 2), c0_Y);
			float nmin = min(((kernel[0] + kernel[1]) / 2), c0_Y);

			// Calculate tanh scale factors for positive/negative overshoot
			float nmax_scale = max((1 / ((nmax - c0_Y) + L_overshoot)), max_scale_lim);
			float nmin_scale = max((1 / ((c0_Y - nmin) + D_overshoot)), max_scale_lim);

			// Soft limit sharpening with tanh, mix to control maximum compression
			sharpdiff = lerp((tanh((max(sharpdiff, 0.0)) * nmax_scale) / nmax_scale), (max(sharpdiff, 0.0)), L_comp_ratio) +
			            lerp((tanh((min(sharpdiff, 0.0)) * nmin_scale) / nmin_scale), (min(sharpdiff, 0.0)), D_comp_ratio);

			OUTPUT[destPos] = float4(src[i + 3][j + 3].rgb + sharpdiff, 1);
		}
	}
}
