Bicubic heightmap interpolation in ISPC

SampleHeightmap_Bicubic.ispc

// MIT License
// 
// Copyright (c) 2023 Voxel Plugin
// 
// 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.

// Useful references:
// https://ennogames.com/blog/bilinear-vs-bicubic-interpolation-and-how-to-get-height-of-terrain-in-unity-from-a-heightmap
// https://www.paulinternet.nl/?page=bicubic

inline float CubicValue(
	const float p0,
	const float p1,
	const float p2,
	const float p3,
	const float X)
{
    return p1 + 0.5f * X * (p2 - p0 + X * (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3 + X * (3.0f * (p1 - p2) + p3 - p0)));
}

export void SampleHeightmap_Bicubic(
	const uniform float ArrayPositionX[],
	const uniform bool bConstPositionX,
	const uniform float ArrayPositionY[],
	const uniform bool bConstPositionY,
	const uniform float BrushScaleXY,
	const uniform float ScaleZ,
	const uniform float OffsetZ,
	const uniform int32 SizeX,
	const uniform int32 SizeY,
	const uniform uint16 Heightmap[],
	const uniform int32 Num,
	uniform float Heights[])
{
	const uniform float InvBrushScaleXY = 1.f / BrushScaleXY;
	const uniform float HalfSizeX = SizeX / 2.f;
	const uniform float HalfSizeY = SizeY / 2.f;
	const uniform float SizeXMinus2 = SizeX - 2.f;
	const uniform float SizeYMinus2 = SizeY - 2.f;

	const uniform int32 NumBlocks = Num / programCount;

	const uniform float p0A = +0.0f, p0B = -0.5f, p0C = +1.0f, p0D = -0.5f;
	const uniform float p1A = +1.0f, p1B = +0.0f, p1C = -2.5f, p1D = +1.5f;
	const uniform float p2A = +0.0f, p2B = +0.5f, p2C = +2.0f, p2D = -1.5f;
	const uniform float p3A = +0.0f, p3B = +0.0f, p3C = -0.5f, p3D = +0.5f;

	for (uniform int32 BlockIndex = 0; BlockIndex < NumBlocks; BlockIndex++)
	{
		const varying int32 Index = programCount * BlockIndex + programIndex;

		varying int32 BaseIndexVector;
		varying float AlphaXVector;
		varying float p0MultiplierVector;
		varying float p1MultiplierVector;
		varying float p2MultiplierVector;
		varying float p3MultiplierVector;

		{
			varying float PositionX = bConstPositionX ? ArrayPositionX[0] : ArrayPositionX[Index];
			varying float PositionY = bConstPositionY ? ArrayPositionY[0] : ArrayPositionY[Index];

			PositionX *= InvBrushScaleXY;
			PositionY *= InvBrushScaleXY;

			PositionX += HalfSizeX;
			PositionY += HalfSizeY;

			const varying float MinXf = clamp(floor(PositionX), 1.f, SizeXMinus2);
			const varying float MinYf = clamp(floor(PositionY), 1.f, SizeYMinus2);

			const varying int32 MinX = (int32)MinXf;
			const varying int32 MinY = (int32)MinYf;

			BaseIndexVector = (MinX - 1) + SizeX * (MinY - 1);
			AlphaXVector = clamp(PositionX - MinXf, 0.f, 1.f);

			const varying float AlphaY = clamp(PositionY - MinYf, 0.f, 1.f);

			p0MultiplierVector = p0A + AlphaY * (p0B + AlphaY * (p0C + AlphaY * p0D));
			p1MultiplierVector = p1A + AlphaY * (p1B + AlphaY * (p1C + AlphaY * p1D));
			p2MultiplierVector = p2A + AlphaY * (p2B + AlphaY * (p2C + AlphaY * p2D));
			p3MultiplierVector = p3A + AlphaY * (p3B + AlphaY * (p3C + AlphaY * p3D));
		}

		varying float Height;

        #pragma unroll
		for (uniform int32 Lane = 0; Lane < programCount; Lane++)
		{
			const uniform int32 BaseIndex = extract(BaseIndexVector, Lane);
			const uniform float p0Multiplier = extract(p0MultiplierVector, Lane);
			const uniform float p1Multiplier = extract(p1MultiplierVector, Lane);
			const uniform float p2Multiplier = extract(p2MultiplierVector, Lane);
			const uniform float p3Multiplier = extract(p3MultiplierVector, Lane);

			const varying float p0 = Heightmap[BaseIndex + SizeX * 0 + programIndex] * p0Multiplier;
			const varying float p1 = Heightmap[BaseIndex + SizeX * 1 + programIndex] * p1Multiplier;
			const varying float p2 = Heightmap[BaseIndex + SizeX * 2 + programIndex] * p2Multiplier;
			const varying float p3 = Heightmap[BaseIndex + SizeX * 3 + programIndex] * p3Multiplier;

#if TARGET_WIDTH == 8
			const varying float A = { p0A, p1A, p2A, p3A, 0.f, 0.f, 0.f, 0.f };
			const varying float B = { p0B, p1B, p2B, p3B, 0.f, 0.f, 0.f, 0.f };
			const varying float C = { p0C, p1C, p2C, p3C, 0.f, 0.f, 0.f, 0.f };
			const varying float D = { p0D, p1D, p2D, p3D, 0.f, 0.f, 0.f, 0.f };
#elif TARGET_WIDTH == 4
			const varying float A = { p0A, p1A, p2A, p3A };
			const varying float B = { p0B, p1B, p2B, p3B };
			const varying float C = { p0C, p1C, p2C, p3C };
			const varying float D = { p0D, p1D, p2D, p3D };
#else
#error Unsupported TARGET_WIDTH
#endif

			const uniform float AlphaX = extract(AlphaXVector, Lane);
			const varying float MultiplierXVector = A + AlphaX * (B + AlphaX * (C + AlphaX * D));

			varying float Value = (p0 + p1 + p2 + p3) * MultiplierXVector;
			Value += rotate(Value, 1);
			Value += rotate(Value, 2);

			Height = insert(Height, Lane, extract(Value, 0));
		}

		Heights[Index] = Height * ScaleZ + OffsetZ;
	}

	foreach(Index = programCount * NumBlocks ... Num)
	{
		varying float PositionX = bConstPositionX ? ArrayPositionX[0] : ArrayPositionX[Index];
		varying float PositionY = bConstPositionY ? ArrayPositionY[0] : ArrayPositionY[Index];

		PositionX *= InvBrushScaleXY;
		PositionY *= InvBrushScaleXY;

		PositionX += HalfSizeX;
		PositionY += HalfSizeY;

		const varying float MinXf = clamp(floor(PositionX), 1.f, SizeXMinus2);
		const varying float MinYf = clamp(floor(PositionY), 1.f, SizeYMinus2);

		const varying float AlphaX = clamp(PositionX - MinXf, 0.f, 1.f);
		const varying float AlphaY = clamp(PositionY - MinYf, 0.f, 1.f);

		const varying int32 MinX = (int32)MinXf;
		const varying int32 MinY = (int32)MinYf;

		#pragma ignore warning(perf)
		const varying float p00 = Heightmap[(MinX - 1) + SizeX * (MinY - 1)];
		#pragma ignore warning(perf)
		const varying float p01 = Heightmap[(MinX + 0) + SizeX * (MinY - 1)];
		#pragma ignore warning(perf)
		const varying float p02 = Heightmap[(MinX + 1) + SizeX * (MinY - 1)];
		#pragma ignore warning(perf)
		const varying float p03 = Heightmap[(MinX + 2) + SizeX * (MinY - 1)];

		#pragma ignore warning(perf)
		const varying float p10 = Heightmap[(MinX - 1) + SizeX * (MinY + 0)];
		#pragma ignore warning(perf)
		const varying float p11 = Heightmap[(MinX + 0) + SizeX * (MinY + 0)];
		#pragma ignore warning(perf)
		const varying float p12 = Heightmap[(MinX + 1) + SizeX * (MinY + 0)];
		#pragma ignore warning(perf)
		const varying float p13 = Heightmap[(MinX + 2) + SizeX * (MinY + 0)];

		#pragma ignore warning(perf)
		const varying float p20 = Heightmap[(MinX - 1) + SizeX * (MinY + 1)];
		#pragma ignore warning(perf)
		const varying float p21 = Heightmap[(MinX + 0) + SizeX * (MinY + 1)];
		#pragma ignore warning(perf)
		const varying float p22 = Heightmap[(MinX + 1) + SizeX * (MinY + 1)];
		#pragma ignore warning(perf)
		const varying float p23 = Heightmap[(MinX + 2) + SizeX * (MinY + 1)];

		#pragma ignore warning(perf)
		const varying float p30 = Heightmap[(MinX - 1) + SizeX * (MinY + 2)];
		#pragma ignore warning(perf)
		const varying float p31 = Heightmap[(MinX + 0) + SizeX * (MinY + 2)];
		#pragma ignore warning(perf)
		const varying float p32 = Heightmap[(MinX + 1) + SizeX * (MinY + 2)];
		#pragma ignore warning(perf)
		const varying float p33 = Heightmap[(MinX + 2) + SizeX * (MinY + 2)];

		const varying float Height = CubicValue(
			CubicValue(p00, p01, p02, p03, AlphaX),
			CubicValue(p10, p11, p12, p13, AlphaX),
			CubicValue(p20, p21, p22, p23, AlphaX),
			CubicValue(p30, p31, p32, p33, AlphaX),
			AlphaY);

		Heights[Index] = Height * ScaleZ + OffsetZ;
	}
}

Godbolt: https://ispc.godbolt.org/z/76h6qKnba