ZFX + Developia

Code: Alles auswählen

//==============================================================================================================================
// convolution 2048.hlsl
// public domain by Krishty, 2008–2011
//==============================================================================================================================
//
// !ARCHITECTURE! marks architecture-dependent optimization opportunities. They should be enclosed in #ifdef #elif #else #endif
//	blocks, so the code base remains the same for all architectures. Check for the following #defines:
//		• AMD_RADEON_HD		AMD Radeon HD 5xxx or higher architecture (4xxx doesn't offer compute shaders)
//
// Some notes on the many discrete Fourier transformations being done here:
//	• All DFTs are realized as fast Fourier transformations with different radixes. Higher radixes require less arithmetic
//		instructions and less group synchronizations, but they also require more registers and read/write scattering and they
//		are more complex. Lower-radix FFTs perform better on current GPUs (2011) for the following reasons:
//			— Compilers fail on large shaders. This applies currently (early 2011) to Microsofts HLSL shader compiler as well as
//				to AMD's shader bytecode compiler. Both produce extremely poor code for shaders with hundreds of instructions:
//				20 % of all instructions are useless; ALU utilization is below 50 %; unneeded scratch registers both in local
//				and global memory are a huge problem.
//			— Read/write scattering is awfully slow and memory bandwidth is a bottleneck. Currently, it is unknown whether this
//				is by concept or a driver bug.
//			— Simple shaders in many threads perform faster than complex shaders in not-so-many threads, even if the total
//				amount of work is bigger (probably because the scheduler can hide latencies better).
//			— Although the number of available GPU registers is quite high (AMD provides at least 32 and at most 512 registers,
//				depending on the number of threads per group), using more than a few registers (6 on AMD hardware) reduces
//				performance drastically.
//		For the future, however, we expect higher-radix FFTs to perform better:
//			— Read/write scattering performance is constantly improving (GPGPU).
//			— GPU performance on complex shaders is constantly improving (more registers, better compilers due to GPGPU).
//	• We use the Stockham auto-sort algorithm. This is an out-of-place algorithm which reads from one source, performs an FFT
//		with an arbitrary radix and writes the result to another array. It has been proposed in "High Performance Discrete
//		Fourier Transforms on Graphics Processors" (Naga K. Govindaraju, Brandon Lloyd, Yuri Dotsenko, Burton Smith, and John
//		Manferdelli, Microsoft Corporation). Basic functionality is copied from the paper.
//	• Most papers propose to invert a FFT by negating the complex exponent. This would, however, break all higher-radix
//		butterflies (because their implementation depends on the twiddle factors being used and therefore on the sign of the
//		complex exponent). It is — far — easier to swap the real and imaginary components of the operands before and after
//		applying an ordinary FFT. This is fully compatible to all FFT butterflies, requires little code change and has nearly no
//		cost when implemented as register swizzling. The final division by the number of samples remains. Source:
//		"Algorithms for Programmers" (Jörg Arndt), chapter 1.7: "Inverse FFT for free".
//	• The group-shared memory in Direct3D 11 is limited to 32 KiB; this is
//			— 8192 'float' values
//			— 4096 complex ('float2') values
//	• For input and output of the FFT algorithm, there are two possibilies:
//		— Using two group-shared arrays ("bank 0" and "bank 1") in a ping-pong pattern. Because different FFTs access different
//			elements, a group thread synchronization is necessary before the banks are swapped. This doubles the space
//			requirement of the algorithm — FFTs are limited to 2048 samples.
//		— Using only one group-shared array and local registers as a cache. This consumes less group-shared memory (now FFTs on
//			4096 samples are possible) but requires twice as many group synchronizations.
//		Although the first possibility has proven to be slightly faster (5 %) on AMD GPUs, it also takes thrice as much time to
//		compile and performs 13.8 % worse on Nvidia GPUs.
//	• The red, green and blue components (where needed) are treated sequentially. This greatly simplifies the shaders, reduces
//		read/write scattering and allows high parallelization.
//		For horizontal FFTs, the input color channels should be stacked vertically and written out horizontally (this behaviour
//		is swapped on inverse FFTs). Rationale: the Fourier transformation of a zero signal is zero, again — therefore, all FFTs
//		on padding can be omitted. As a fortunate coincidence, out-of-bounds reads from DXGI resources return 0, too — the
//		padding can be omitted completely and the texture can be cropped down to the input's actual size. This saves a lot of
//		texture space and bandwidth (73 % on vertical passes with full HD — 1080 pixels are read and written instead of 4096)
//		and the need to clear the texture with zeroes (which limits the bandwidth, too).
//
//==============================================================================================================================



//==============================================================================================================================
// GENERIC FUNCTIONS
// This code handles mostly FFTs and is the same for all convolution implementations. It can be re-used by simply adjusting the
//	"IMAGE_SIDELENGTH" constant.
//==============================================================================================================================
// #define AMD_RADEON_HD
#define IMAGE_SIDELENGTH 2048 // group size attributes require literal constants
static const uint imageSidelength = IMAGE_SIDELENGTH;



// !ARCHITECTURE! Scratch memory for the current DFT. Never read or write it directly — use "loadGroupShared()" and
//	"writeGroupShared" to read and write, respectively. This allows quick adjustion of the group-shared memory layout for
//	different architectures, e.g. to avoid bank conflicts.
#if defined(AMD_RADEON_HD)
	//	• on AMD Radeon HD architecture, storing the values as 'float2's has proven to be 5.8 % faster than seperating real and
	//		imaginary components
	//	• on AMD Radeon HD architecture, the ping-pong pattern is >5 % faster than using a single group-shared array (but it
	//		requires more than twice the time to compile) — this is probably due to the compiler having problems transferring
	//		instructions over group thread barriers (a ping-pong pattern halves the number of required synchronizations)

	groupshared float2 groupSharedValues[2][imageSidelength];
	static uint currentGSMOutBanksIndex = 1; // first writing to bank 1 saves two cycles; I don't know why

	float2 loadGroupShared(
		const uint index
	) {
		return groupSharedValues[currentGSMOutBanksIndex ^ 1][index];
	}
	//..........................................................................................................................
	void writeGroupShared(
		const uint index,
		const float2 value
	) {
		groupSharedValues[currentGSMOutBanksIndex][index] = value;
		return;
	}
	//..........................................................................................................................
	void switchGSMBank() {
		currentGSMOutBanksIndex ^= 1;
		return;
	}

#else // not AMD Radeon HD:
	//	• on Nvidia GeForce GT architecture, a single group-shared array is 13.8 % faster than the ping-pong pattern — this is
	//		probably due to a group-shared memory size limitation there; halving the GSM consumption doubles the effective
	//		wavefront size

	groupshared float2 groupSharedValues[imageSidelength];

	float2 loadGroupShared(
		const uint index
	) {
		return groupSharedValues[index];
	}
	//..........................................................................................................................
	void writeGroupShared(
		const uint index,
		const float2 value
	) {
		groupSharedValues[index] = value;
		return;
	}

#endif // default hardware



// For better readability of function parameters.
static const bool invert = true;
static const bool dontInvert = false;
static const bool horizontally = false;
static const bool vertically = true;



// The input and output buffers. Their registers overlap because all shaders expect only one input and write to one output (with
//	the exception of the convolution itself, which reads from the image's DFT as well as from the kernel's DFT).

// Stores the glare with its three color channels vertically stacked. At the beginning of the glare operator, this is the
//	resolved (and possibly downsampled) scene.
RWTexture2DArray<float>		glareRW			: register(u0);
Texture2DArray<float>		glareRO			: register(t0);

// Stores the point spread function with its three color channels vertically stacked.
Texture2DArray<float>		PSFRO			: register(t0);

// Stores the DFT of the aperture's point spread function. The real and imaginary components are stored in the X and Y component
//	and its three color channels are vertically stacked.
RWTexture2DArray<float2>	kernelsDFTRW	: register(u0);
Texture2DArray<float2>		kernelsDFTRO	: register(t1);

// Stores the DFT of the scene in the same layout as the aperture PSF DFT's texture.
RWTexture2DArray<float2>	imagesDFTRW		: register(u0);
Texture2DArray<float2>		imagesDFTRO		: register(t0);



//------------------------------------------------------------------------------------------------------------------------------
// Swaps the real and imaginary components of the given complex number if the given boolean expression evaluates TRUE.
// Allows re-use of FFT functions for iFFT by just switching a boolean expression.
//------------------------------------------------------------------------------------------------------------------------------
float2 swapIf(
	const bool		doOrDont,
	const float2	value
) {
	return doOrDont ? value.yx : value.xy;
}



//------------------------------------------------------------------------------------------------------------------------------
// Multiplies the two given complex numbers.
//------------------------------------------------------------------------------------------------------------------------------
float2 multiplyComplex(
	const float2 a,
	const float2 b
) {
	// !ARCHITECTURE! There are several ways to express a complex multiplication:
	return float2(

#		if defined(AMD_RADEON_HD)
			// Saves one out of 100 instructions and — sometimes — one register. Since the main problem on AMD Radeon HD
			//	hardware is register pressure, this performs best.
			dot(float2(a.x, -a.y), float2(b.x, b.y)),
			dot(float2(a.x, a.y), float2(b.y, b.x))
#		else
			// Saves another instruction, but does not lower register pressure. Decide after profiling.
			mad(a.x, b.x, -a.y * b.y),
			mad(a.x, b.y, a.y * b.x)
			// The default procedure. Source: http://en.wikipedia.org/wiki/Complex_numbers#Multiplication_and_division.
			//	a.x * b.x - a.y * b.y,
			//	a.x * b.y + a.y * b.x
#		endif
	);
}



//------------------------------------------------------------------------------------------------------------------------------
// Computes the base position for an FFT's output.
//------------------------------------------------------------------------------------------------------------------------------
uint baseIndexFor(
	const uint radixsIndex,
	const uint stepSizeInSamples,
	const uint radix
) {
	// Source: "High Performance Discrete Fourier Transforms on Graphics Processors" (Naga K. Govindaraju, Brandon Lloyd, Yuri
	//	Dotsenko, Burton Smith, and John Manferdelli; Microsoft Corporation), fig. 2.
	return (radixsIndex / stepSizeInSamples) * stepSizeInSamples * radix + (radixsIndex % stepSizeInSamples);
}



//------------------------------------------------------------------------------------------------------------------------------
// Performs a fast Fourier transformation on group-shared memory.
//	• Radix 2 or 4.
//	• Can be used for an inverse FFT.
// This function can be executed in parallel — subsequent calls with the same step size but different radix indices will not
//	overlap. It reads and writes out-of-place — after the input has been written, the group-shared memory must be
//	synchronized. This is automatically done before this routine returns.
//------------------------------------------------------------------------------------------------------------------------------
void FFTOnGroupSharedMemory(
	const uint radix,
	const bool inverse,
	const uint radixsIndex,
	const uint stepSizeInSamples
) {
	// Source: "High Performance Discrete Fourier Transforms on Graphics Processors" (Naga K. Govindaraju, Brandon Lloyd, Yuri
	//	Dotsenko, Burton Smith, and John Manferdelli, Microsoft Corporation), fig. 2.
	const float angle = -6.2831853071795865f / float(stepSizeInSamples * radix) * float(radixsIndex % stepSizeInSamples);
	const uint baseIndex = baseIndexFor(radixsIndex, stepSizeInSamples, radix);

#	if defined(AMD_RADEON_HD)
		// On AMD Radeon HD hardware, a ping-pong pattern is faster.

		if(2 == radix) {
			static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// Load two samples into a local cache and multiply them with their twiddle factors.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates
				//	to (1, 0) and the complex multiplication yields the identity.
				swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex))),
					float2(cos(angle), sin(angle))
				)
			};

			// Use a radix-2 butterfly to transform the samples. Source: http://en.wikipedia.org/wiki/Butterfly_diagram.
			const float2 value0 = cache[0];
			const float2 value1 = cache[1];
			cache[0] = value0 + value1;
			cache[1] = value0 - value1;

			// Write the result out-of-order back.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), swapIf(inverse, cache[0]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), swapIf(inverse, cache[1]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();
			switchGSMBank();

		} else if(4 == radix) {
			static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// Load two samples into a local cache and multiply them with their twiddle factors.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates
				//	to (1, 0) and the complex multiplication yields the identity.
					swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex))),
					float2(cos(1.0f * angle), sin(1.0f * angle))
				),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(2, imageSidelength / radix, radixsIndex))),
					float2(cos(2.0f * angle), sin(2.0f * angle))
				),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(3, imageSidelength / radix, radixsIndex))),
					float2(cos(3.0f * angle), sin(3.0f * angle))
				)
			};

			// Use a radix-4 butterfly to transform the samples. Source: unknown; found it somewhere on the internet.
			const float2 value0 = cache[0];
			const float2 value1 = cache[1];
			const float2 value2 = cache[2];
			const float2 value3 = cache[3];

			cache[0] = float2(
				value0.x + value1.x + value2.x + value3.x,
				value0.y + value1.y + value2.y + value3.y
			);
			cache[1] = float2(
				value0.x + value1.y - value2.x - value3.y,
				value0.y - value1.x - value2.y + value3.x
			);
			cache[2] = float2(
				value0.x - value1.x + value2.x - value3.x,
				value0.y - value1.y + value2.y - value3.y
			);
			cache[3] = float2(
				value0.x - value1.y - value2.x + value3.y,
				value0.y + value1.x - value2.y - value3.x
			);

			// Write the result out-of-order back.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), swapIf(inverse, cache[0]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), swapIf(inverse, cache[1]));
			writeGroupShared(mad(2, stepSizeInSamples, baseIndex), swapIf(inverse, cache[2]));
			writeGroupShared(mad(3, stepSizeInSamples, baseIndex), swapIf(inverse, cache[3]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();
			switchGSMBank();

		}

#	else // default hardware:
		// Cache all group-shared values in registers before performing the transformation.

		if(2 == radix) {
			static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// Load two samples into a local cache and multiply them with their twiddle factors.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates
				//	to (1, 0) and the complex multiplication yields the identity.
				swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex))),
					float2(cos(angle), sin(angle))
				)
			};

			// Wait for all threads to complete their reading before the result is written back.
			GroupMemoryBarrierWithGroupSync();

			// Use a radix-2 butterfly to transform the samples. Source: http://en.wikipedia.org/wiki/Butterfly_diagram.
			const float2 value0 = cache[0];
			const float2 value1 = cache[1];
			cache[0] = value0 + value1;
			cache[1] = value0 - value1;

			// Write the result out-of-order back.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), swapIf(inverse, cache[0]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), swapIf(inverse, cache[1]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();

		} else if(4 == radix) {
			static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// Load two samples into a local cache and multiply them with their twiddle factors.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates
				//	to (1, 0) and the complex multiplication yields the identity.
					swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex))),
					float2(cos(1.0f * angle), sin(1.0f * angle))
				),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(2, imageSidelength / radix, radixsIndex))),
					float2(cos(2.0f * angle), sin(2.0f * angle))
				),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(3, imageSidelength / radix, radixsIndex))),
					float2(cos(3.0f * angle), sin(3.0f * angle))
				)
			};

			// Wait for all threads to complete their reading before the result is written back.
			GroupMemoryBarrierWithGroupSync();

			// Use a radix-4 butterfly to transform the samples. Source: unknown; found it somewhere on the internet.
			const float2 value0 = cache[0];
			const float2 value1 = cache[1];
			const float2 value2 = cache[2];
			const float2 value3 = cache[3];

			cache[0] = float2(
				value0.x + value1.x + value2.x + value3.x,
				value0.y + value1.y + value2.y + value3.y
			);
			cache[1] = float2(
				value0.x + value1.y - value2.x - value3.y,
				value0.y - value1.x - value2.y + value3.x
			);
			cache[2] = float2(
				value0.x - value1.x + value2.x - value3.x,
				value0.y - value1.y + value2.y - value3.y
			);
			cache[3] = float2(
				value0.x - value1.y - value2.x + value3.y,
				value0.y + value1.x - value2.y - value3.x
			);

			// Write the result out-of-order back.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), swapIf(inverse, cache[0]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), swapIf(inverse, cache[1]));
			writeGroupShared(mad(2, stepSizeInSamples, baseIndex), swapIf(inverse, cache[2]));
			writeGroupShared(mad(3, stepSizeInSamples, baseIndex), swapIf(inverse, cache[3]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();

		}

#	endif // default hardware

	return;
}



//------------------------------------------------------------------------------------------------------------------------------
// Performs two parallel FFTs on group-shared memory.
// This function is necessary because some powers of two can not be expressed by a single radix — e.g. 512 must be transformed
//	with radix 4-4-4-4-2. Since the shader cannot switch its thread group size while in execution, two lower-radix FFTs must be
//	performed. This would, however, enforce two additional group synchronizations (although both lower-radix FFTs do not
//	overlap). This function's purpose is to offer two parallel lower-radix FFTs without unnecessary synchronizations.
//	• Radix 2.
//	• Can be used for an inverse FFT.
//------------------------------------------------------------------------------------------------------------------------------
void TwoFFTsOnGroupSharedMemory(
	const uint radix,
	const bool inverse,
	const uint radixsSuperiorIndex,
	const uint stepSizeInSamples
) {
	const uint2 radixsIndex = uint2(
		mad(2, radixsSuperiorIndex, 0),
		mad(2, radixsSuperiorIndex, 1)
	);
	const float2 angle = float2(
		-6.2831853071795865f / float(stepSizeInSamples * radix) * float(radixsIndex[0] % stepSizeInSamples),
		-6.2831853071795865f / float(stepSizeInSamples * radix) * float(radixsIndex[1] % stepSizeInSamples)
	);

	const uint2 baseIndex = uint2(
		baseIndexFor(radixsIndex[0], stepSizeInSamples, radix),
		baseIndexFor(radixsIndex[1], stepSizeInSamples, radix)
	);

#	if defined(AMD_RADEON_HD)
		// On AMD Radeon HD hardware, a ping-pong pattern is faster.

		if(2 == radix) {
			static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// Load two samples into a local cache and multiply them with their twiddle factors.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
			//	(1, 0) and the complex multiplication yields the identity.
			const float2 cacheA[radix] = {
					swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex[0]))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex[0]))),
					float2(cos(angle[0]), sin(angle[0]))
				)
			};
			const float2 cacheB[radix] = {
					swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex[1]))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex[1]))),
					float2(cos(angle[1]), sin(angle[1]))
				)
			};

			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex[0]), swapIf(inverse, cacheA[0] + cacheA[1]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex[0]), swapIf(inverse, cacheA[0] - cacheA[1]));
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex[1]), swapIf(inverse, cacheB[0] + cacheB[1]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex[1]), swapIf(inverse, cacheB[0] - cacheB[1]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();
			switchGSMBank();

		}

#	else // default hardware:
		// Cache all group-shared values in registers before performing the transformation.

		if(2 == radix) {
			static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// Load two samples into a local cache and multiply them with their twiddle factors.
			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
			//	(1, 0) and the complex multiplication yields the identity.
			const float2 cacheA[radix] = {
					swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex[0]))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex[0]))),
					float2(cos(angle[0]), sin(angle[0]))
				)
			};
			const float2 cacheB[radix] = {
					swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex[1]))),
				multiplyComplex(
					swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex[1]))),
					float2(cos(angle[1]), sin(angle[1]))
				)
			};

			// Wait for all threads to complete their reading before the result is written back.
			GroupMemoryBarrierWithGroupSync();

			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex[0]), swapIf(inverse, cacheA[0] + cacheA[1]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex[0]), swapIf(inverse, cacheA[0] - cacheA[1]));
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex[1]), swapIf(inverse, cacheB[0] + cacheB[1]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex[1]), swapIf(inverse, cacheB[0] - cacheB[1]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();

		}

#	endif // default hardware

	return;
}



//------------------------------------------------------------------------------------------------------------------------------
// Performs a fast Fourier transformation from a complex texture to group-shared memory.
//	• Radix 2 or 4.
//	• Can be used for an inverse FFT.
//	• The input is read as complex numbers from a channel of the given source texture (either horizontally or vertically).
//------------------------------------------------------------------------------------------------------------------------------
void FFTFromComplexTexture(
	const uint				radix,
	const bool				inverse,
	Texture2DArray<float2>	source,
	const uint				channelsIndex,
	const bool				vertical,
	const uint				signalsIndex,
	const uint				radixsIndex
) {
	const uint3 basePosition = vertical
		? uint3(signalsIndex, radixsIndex, channelsIndex)
		: uint3(radixsIndex, signalsIndex, channelsIndex);
	const uint3 delta = vertical
		? uint3(0, imageSidelength / radix, 0)
		: uint3(imageSidelength / radix, 0, 0);

	if(2 == radix) {
		static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// With a step size of 1, the twiddle factors can be omitted completely — a multiplication with the complex number
		//		cos(-2×Pi) + i × sin(-2×Pi)
		//	yields the identity.
		// (Unroll manually — loop expressions with global memory access are not optimized well.)
		const float2 cache[radix] = {
			swapIf(inverse, source[mad(0, delta, basePosition)]),
			swapIf(inverse, source[mad(1, delta, basePosition)])
		};

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		writeGroupShared(mad(radix, radixsIndex, 0), swapIf(inverse, cache[0] + cache[1]));
		writeGroupShared(mad(radix, radixsIndex, 1), swapIf(inverse, cache[0] - cache[1]));

		// Wait for all threads to complete their writing.
		GroupMemoryBarrierWithGroupSync();

	} else if(4 == radix) {
		static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// With a step size of 1, the twiddle factors can be omitted completely — a multiplication with the complex number
		//		cos(-2×Pi) + i × sin(-2×Pi)
		//	yields the identity.
		// (Unroll manually — loop expressions with global memory access are not optimized well.)
		const float2 cache[radix] = {
			swapIf(inverse, source[mad(0, delta, basePosition)]),
			swapIf(inverse, source[mad(1, delta, basePosition)]),
			swapIf(inverse, source[mad(2, delta, basePosition)]),
			swapIf(inverse, source[mad(3, delta, basePosition)])
		};

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		writeGroupShared(mad(radix, radixsIndex, 0), swapIf(inverse, float2(
			cache[0].x + cache[1].x + cache[2].x + cache[3].x,
			cache[0].y + cache[1].y + cache[2].y + cache[3].y
		)));
		writeGroupShared(mad(radix, radixsIndex, 1), swapIf(inverse, float2(
			cache[0].x + cache[1].y - cache[2].x - cache[3].y,
			cache[0].y - cache[1].x - cache[2].y + cache[3].x
		)));
		writeGroupShared(mad(radix, radixsIndex, 2), swapIf(inverse, float2(
			cache[0].x - cache[1].x + cache[2].x - cache[3].x,
			cache[0].y - cache[1].y + cache[2].y - cache[3].y
		)));
		writeGroupShared(mad(radix, radixsIndex, 3), swapIf(inverse, float2(
			cache[0].x - cache[1].y - cache[2].x + cache[3].y,
			cache[0].y + cache[1].x - cache[2].y - cache[3].x
		)));

		// Wait for all threads to complete their writing.
		GroupMemoryBarrierWithGroupSync();

	}

#	if defined(AMD_RADEON_HD)
		switchGSMBank(); // on AMD Radeon HD hardware, a ping-pong pattern is faster
#	endif
	return;
}



//------------------------------------------------------------------------------------------------------------------------------
// Performs a fast Fourier transformation from a packed real texture to group-shared memory.
// Used to perform DFTs on temporary results which had previously been stored in textures.
//	• Radix 2 or 4.
//	• The input is read as real numbers from a channel of the given texture (either horizontally or vertically).
//------------------------------------------------------------------------------------------------------------------------------
void FFTFromRealTexture(
	const uint				radix,
	Texture2DArray<float>	source,
	const uint				channelsIndex,
	const bool				vertical,
	const uint				signalsIndex,
	const uint				radixsIndex
) {
	const uint3 basePosition = vertical
		? uint3(signalsIndex, radixsIndex, channelsIndex)
		: uint3(radixsIndex, signalsIndex, channelsIndex);
	const uint3 delta = vertical
		? uint3(0, imageSidelength / radix, 0)
		: uint3(imageSidelength / radix, 0, 0);

	if(2 == radix) {
		static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// With a step size of 1, the twiddle factors can be omitted completely — a multiplication with the complex number
		//		cos(-2×Pi) + i × sin(-2×Pi)
		//	yields the identity.
		// (Unroll manually — loop expressions with global memory access are not optimized well.)
		const float cache[radix] = {
			source[mad(0, delta, basePosition)],
			source[mad(1, delta, basePosition)]
		};

		// Save many arithmetic instructions through complex multiplication with pure real numbers.
		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		writeGroupShared(mad(radix, radixsIndex, 0), float2(cache[0] + cache[1], 0.0f));
		writeGroupShared(mad(radix, radixsIndex, 1), float2(cache[0] - cache[1], 0.0f));

		// Wait for all threads to complete their writing.
		GroupMemoryBarrierWithGroupSync();

	} else if(4 == radix) {
		static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// With a step size of 1, the twiddle factors can be omitted completely — a multiplication with the complex number
		//		cos(-2×Pi) + i × sin(-2×Pi)
		//	yields the identity.
		// (Unroll manually — loop expressions with global memory access are not optimized well.)
		const float cache[radix] = {
			source[mad(0, delta, basePosition)],
			source[mad(1, delta, basePosition)],
			source[mad(2, delta, basePosition)],
			source[mad(3, delta, basePosition)]
		};

		// Save many arithmetic instructions through complex multiplication with pure real numbers.
		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		writeGroupShared(mad(radix, radixsIndex, 0), float2(cache[0] + cache[1] + cache[2] + cache[3], 0.0f));
		writeGroupShared(mad(radix, radixsIndex, 1), float2(cache[0] - cache[2], -cache[1] + cache[3]));
		writeGroupShared(mad(radix, radixsIndex, 2), float2(cache[0] - cache[1] + cache[2] - cache[3], 0.0f));
		writeGroupShared(mad(radix, radixsIndex, 3), float2(cache[0] - cache[2], cache[1] - cache[3]));

		// Wait for all threads to complete their writing.
		GroupMemoryBarrierWithGroupSync();

	}

#	if defined(AMD_RADEON_HD)
		switchGSMBank(); // on AMD Radeon HD hardware, a ping-pong pattern is faster
#	endif
	return;
}



//------------------------------------------------------------------------------------------------------------------------------
// Performs a fast Fourier transformation from group-shared memory to a complex texture.
// Used to store temporal results like horizontal DFTs.
//	• Radix 2 or 4.
//	• The input is read from the "inBank" bank of group-shared memory.
//	• The output is written as complex numbers to a channel of the given texture (either horizontally or vertically).
//------------------------------------------------------------------------------------------------------------------------------
void FFTToComplexTexture(
	const uint					radix,
	const bool					inverse,
	RWTexture2DArray<float2>	destination,
	const uint					channelsIndex,
	const bool					vertical,
	const uint					signalsIndex,
	const uint					radixsIndex,
	const uint					stepSizeInSamples
) {
	const float angle = -6.2831853071795865f / float(stepSizeInSamples * radix) * float(radixsIndex % stepSizeInSamples);
	const uint baseIndex = baseIndexFor(radixsIndex, stepSizeInSamples, radix); // will be optimized away
	const uint3 basePosition = vertical
		? uint3(signalsIndex, baseIndex, channelsIndex)
		: uint3(baseIndex, signalsIndex, channelsIndex);
	const uint3 delta = vertical
		? uint3(0, stepSizeInSamples, 0)
		: uint3(stepSizeInSamples, 0, 0);

	if(2 == radix) {
		static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		const float2 cache[radix] = {
			// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
			//	(1, 0) and the complex multiplication yields the identity.
				swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex))),
			multiplyComplex(
				swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex))),
				float2(cos(1.0f * angle), sin(1.0f * angle))
			)
		};

		// (Unroll manually — loop expressions with global memory access are not optimized well.)
		destination[mad(0, delta, basePosition)] = swapIf(inverse, cache[0] + cache[1]);
		destination[mad(1, delta, basePosition)] = swapIf(inverse, cache[0] - cache[1]);

	} else if(4 == radix) {
		static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		const float2 cache[radix] = {
			// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
			//	(1, 0) and the complex multiplication yields the identity.
				swapIf(inverse, loadGroupShared(mad(0, imageSidelength / radix, radixsIndex))),
			multiplyComplex(
				swapIf(inverse, loadGroupShared(mad(1, imageSidelength / radix, radixsIndex))),
				float2(cos(1.0f * angle), sin(1.0f * angle))
			),
			multiplyComplex(
				swapIf(inverse, loadGroupShared(mad(2, imageSidelength / radix, radixsIndex))),
				float2(cos(2.0f * angle), sin(2.0f * angle))
			),
			multiplyComplex(
				swapIf(inverse, loadGroupShared(mad(3, imageSidelength / radix, radixsIndex))),
				float2(cos(3.0f * angle), sin(3.0f * angle))
			)
		};

		// (Unroll manually — loop expressions with global memory access are not optimized well.)
		destination[mad(0, delta, basePosition)] = swapIf(inverse, float2(
														cache[0].x + cache[1].x + cache[2].x + cache[3].x,
														cache[0].y + cache[1].y + cache[2].y + cache[3].y
													));
		destination[mad(1, delta, basePosition)] = swapIf(inverse, float2(
														cache[0].x + cache[1].y - cache[2].x - cache[3].y,
														cache[0].y - cache[1].x - cache[2].y + cache[3].x
													));
		destination[mad(2, delta, basePosition)] = swapIf(inverse, float2(
														cache[0].x - cache[1].x + cache[2].x - cache[3].x,
														cache[0].y - cache[1].y + cache[2].y - cache[3].y
													));
		destination[mad(3, delta, basePosition)] = swapIf(inverse, float2(
														cache[0].x - cache[1].y - cache[2].x + cache[3].y,
														cache[0].y + cache[1].x - cache[2].y - cache[3].x
													));

	}

	return;
}



//------------------------------------------------------------------------------------------------------------------------------
// Performs a fast Fourier transformation on group-shared memory and multiplies the result with complex numbers read from the
//	PSF DFT's texture (vertically).
// Used to combine the scene's DFT with the PSF's DFT before transforming back.
//	• Radix 2 or 4.
//	• The input is read from group-shared memory.
//	• The multiplicands are read vertically from the PSF DFT's texture according to the radix's index and the index in the
//		signal.
//------------------------------------------------------------------------------------------------------------------------------
void FFTOnGroupSharedMemoryWithVerticalPSFMultiplication(
	const uint			radix,
	const uint			channelsIndex,
	const uint			signalsIndex,
	const uint			radixsIndex,
	const uint			stepSizeInSamples
) {
	const float angle = -6.2831853071795865f / float(stepSizeInSamples * radix) * float(radixsIndex % stepSizeInSamples);
	const uint baseIndex = baseIndexFor(radixsIndex, stepSizeInSamples, radix); // will be optimized away
	const uint3 basePosition = uint3(signalsIndex, baseIndex, channelsIndex);
	const uint3 delta = uint3(0, stepSizeInSamples, 0);

#	if defined(AMD_RADEON_HD)
		// On AMD Radeon HD hardware, a ping-pong pattern is faster.

		if(2 == radix) {
			static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
				//	(1, 0) and the complex multiplication yields the identity.
					loadGroupShared(mad(0, imageSidelength / radix, radixsIndex)),
				multiplyComplex(
					loadGroupShared(mad(1, imageSidelength / radix, radixsIndex)),
					float2(cos(1.0f * angle), sin(1.0f * angle))
				)
			};

			// Multiply with the PSF's DFT and write to the scene's DFT texture.
			// (Unroll manually — loop expressions with global memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), multiplyComplex(
															cache[0] + cache[1],
															kernelsDFTRO[mad(0, delta, basePosition)]
														));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), multiplyComplex(
															cache[0] - cache[1],
															kernelsDFTRO[mad(1, delta, basePosition)]
														));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();
			switchGSMBank();

		} else if(4 == radix) {
			static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			const float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
				//	(1, 0) and the complex multiplication yields the identity.
					loadGroupShared(mad(0, imageSidelength / radix, radixsIndex)),
				multiplyComplex(
					loadGroupShared(mad(1, imageSidelength / radix, radixsIndex)),
					float2(cos(1.0f * angle), sin(1.0f * angle))
				),
				multiplyComplex(
					loadGroupShared(mad(2, imageSidelength / radix, radixsIndex)),
					float2(cos(2.0f * angle), sin(2.0f * angle))
				),
				multiplyComplex(
					loadGroupShared(mad(3, imageSidelength / radix, radixsIndex)),
					float2(cos(3.0f * angle), sin(3.0f * angle))
				)
			};

			// Multiply with the PSF's DFT and write to the scene's DFT texture.
			// (Unroll manually — loop expressions with global memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x + cache[1].x + cache[2].x + cache[3].x,
															cache[0].y + cache[1].y + cache[2].y + cache[3].y
														), kernelsDFTRO[mad(0, delta, basePosition)]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x + cache[1].y - cache[2].x - cache[3].y,
															cache[0].y - cache[1].x - cache[2].y + cache[3].x
														), kernelsDFTRO[mad(1, delta, basePosition)]));
			writeGroupShared(mad(2, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x - cache[1].x + cache[2].x - cache[3].x,
															cache[0].y - cache[1].y + cache[2].y - cache[3].y
														), kernelsDFTRO[mad(2, delta, basePosition)]));
			writeGroupShared(mad(3, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x - cache[1].y - cache[2].x + cache[3].y,
															cache[0].y + cache[1].x - cache[2].y - cache[3].x
														), kernelsDFTRO[mad(3, delta, basePosition)]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();
			switchGSMBank();

		}

#	else // default hardware:
		// Cache all values in registers before performing the transformation

		if(2 == radix) {
			static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
				//	(1, 0) and the complex multiplication yields the identity.
					loadGroupShared(mad(0, imageSidelength / radix, radixsIndex)),
				multiplyComplex(
					loadGroupShared(mad(1, imageSidelength / radix, radixsIndex)),
					float2(cos(1.0f * angle), sin(1.0f * angle))
				)
			};

			// Wait for all threads to complete their reading.
			GroupMemoryBarrierWithGroupSync();

			// Multiply with the PSF's DFT and write to the scene's DFT texture.
			// (Unroll manually — loop expressions with global memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), multiplyComplex(
															cache[0] + cache[1],
															kernelsDFTRO[mad(0, delta, basePosition)]
														));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), multiplyComplex(
															cache[0] - cache[1],
															kernelsDFTRO[mad(1, delta, basePosition)]
														));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();

		} else if(4 == radix) {
			static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

			// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
			const float2 cache[radix] = {
				// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
				//	(1, 0) and the complex multiplication yields the identity.
					loadGroupShared(mad(0, imageSidelength / radix, radixsIndex)),
				multiplyComplex(
					loadGroupShared(mad(1, imageSidelength / radix, radixsIndex)),
					float2(cos(1.0f * angle), sin(1.0f * angle))
				),
				multiplyComplex(
					loadGroupShared(mad(2, imageSidelength / radix, radixsIndex)),
					float2(cos(2.0f * angle), sin(2.0f * angle))
				),
				multiplyComplex(
					loadGroupShared(mad(3, imageSidelength / radix, radixsIndex)),
					float2(cos(3.0f * angle), sin(3.0f * angle))
				)
			};

			// Wait for all threads to complete their reading.
			GroupMemoryBarrierWithGroupSync();

			// Multiply with the PSF's DFT and write to the scene's DFT texture.
			// (Unroll manually — loop expressions with global memory access are not optimized well.)
			writeGroupShared(mad(0, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x + cache[1].x + cache[2].x + cache[3].x,
															cache[0].y + cache[1].y + cache[2].y + cache[3].y
														), kernelsDFTRO[mad(0, delta, basePosition)]));
			writeGroupShared(mad(1, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x + cache[1].y - cache[2].x - cache[3].y,
															cache[0].y - cache[1].x - cache[2].y + cache[3].x
														), kernelsDFTRO[mad(1, delta, basePosition)]));
			writeGroupShared(mad(2, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x - cache[1].x + cache[2].x - cache[3].x,
															cache[0].y - cache[1].y + cache[2].y - cache[3].y
														), kernelsDFTRO[mad(2, delta, basePosition)]));
			writeGroupShared(mad(3, stepSizeInSamples, baseIndex), multiplyComplex(float2(
															cache[0].x - cache[1].y - cache[2].x + cache[3].y,
															cache[0].y + cache[1].x - cache[2].y - cache[3].x
														), kernelsDFTRO[mad(3, delta, basePosition)]));

			// Wait for all threads to complete their writing.
			GroupMemoryBarrierWithGroupSync();

		}
#	endif // default hardware

	return;
}



//------------------------------------------------------------------------------------------------------------------------------
// Performs an inverse fast Fourier transformation from values in group-shared memory and writes the result as packed reals into
//	a texture.
// Used to extract the real values from the inverse DFT and to finalize the inverse transformation.
//	• Radix 2 or 4.
//	• The input is read from the group-shared memory.
//	• The result is written as a packed real number (using the real value of the complex result) to the given texture.
//------------------------------------------------------------------------------------------------------------------------------
void iFFTToRealTexture(
	const uint				radix,
	RWTexture2DArray<float>	destination,
	const uint				channelsIndex,
	const bool				vertical,
	const uint				signalsIndex,
	const uint				radixsIndex,
	const uint				stepSizeInSamples
) {
	// Source: "High Performance Discrete Fourier Transforms on Graphics Processors" (Naga K. Govindaraju, Brandon Lloyd, Yuri
	//	Dotsenko, Burton Smith, and John Manferdelli, Microsoft Corporation), fig. 2.
	const float angle = -6.2831853071795865f / float(stepSizeInSamples * radix) * float(radixsIndex % stepSizeInSamples);
	const uint baseIndex = baseIndexFor(radixsIndex, stepSizeInSamples, radix); // will be optimized away
	const uint3 basePosition = vertical
		? uint3(signalsIndex, baseIndex, channelsIndex)
		: uint3(baseIndex, signalsIndex, channelsIndex);
	const uint3 delta = vertical
		? uint3(0, stepSizeInSamples, 0)
		: uint3(stepSizeInSamples, 0, 0);

	if(2 == radix) {
		static const uint radix = 2; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		float2 cache[radix] = {
			// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
			//	(1, 0) and the complex multiplication yields the identity.
				loadGroupShared(mad(0, imageSidelength / radix, radixsIndex)).yx, // swap because inverse
			multiplyComplex(
				loadGroupShared(mad(1, imageSidelength / radix, radixsIndex)).yx, // swap because inverse
				float2(cos(angle), sin(angle))
			)
		};

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		destination[mad(0, delta, basePosition)] = (cache[0].y + cache[1].y) / imageSidelength;
		destination[mad(1, delta, basePosition)] = (cache[0].y - cache[1].y) / imageSidelength;

	} else if(4 == radix) {
		static const uint radix = 4; // redeclare 'static' — HLSL doesn't accept parameters as array dimensions

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		float2 cache[radix] = {
			// The first complex multiplication can be omitted: The angle is always zero; the complex exponent evaluates to
			//	(1, 0) and the complex multiplication yields the identity.
				loadGroupShared(mad(0, imageSidelength / radix, radixsIndex)).yx, // swap because inverse
			multiplyComplex(
				loadGroupShared(mad(1, imageSidelength / radix, radixsIndex)).yx, // swap because inverse
				float2(cos(1.0f * angle), sin(1.0f * angle))
			),
			multiplyComplex(
				loadGroupShared(mad(2, imageSidelength / radix, radixsIndex)).yx, // swap because inverse
				float2(cos(2.0f * angle), sin(2.0f * angle))
			),
			multiplyComplex(
				loadGroupShared(mad(3, imageSidelength / radix, radixsIndex)).yx, // swap because inverse
				float2(cos(3.0f * angle), sin(3.0f * angle))
			)
		};

		// (Unroll manually — loop expressions with group-shared memory access are not optimized well.)
		destination[mad(0, delta, basePosition)] = (cache[0].y + cache[1].y + cache[2].y + cache[3].y) / (imageSidelength * imageSidelength);
		destination[mad(1, delta, basePosition)] = (cache[0].y - cache[1].x - cache[2].y + cache[3].x) / (imageSidelength * imageSidelength);
		destination[mad(2, delta, basePosition)] = (cache[0].y - cache[1].y + cache[2].y - cache[3].y) / (imageSidelength * imageSidelength);
		destination[mad(3, delta, basePosition)] = (cache[0].y + cache[1].x - cache[2].y - cache[3].x) / (imageSidelength * imageSidelength);

	}

}











[numthreads(IMAGE_SIDELENGTH / 4, 1, 1)] // 4 pixels per thread — best for radix 4
void kernelHorizontal(
	const uint3 indexInSignal : SV_DispatchThreadID
) {
	const uint currentChannelsIndex = indexInSignal.z;

	FFTFromRealTexture			(4, PSFRO, currentChannelsIndex, horizontally, indexInSignal.y, indexInSignal.x);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.x,   4);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.x,  16);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.x,  64);
	TwoFFTsOnGroupSharedMemory	(2, dontInvert, indexInSignal.x, 256);
	FFTToComplexTexture			(4, dontInvert, kernelsDFTRW, currentChannelsIndex, horizontally, indexInSignal.y, indexInSignal.x, 512);

	return;
}

[numthreads(1, IMAGE_SIDELENGTH / 4, 1)] // 4 pixels per thread — best for radix 4
void kernelVertical(
	const uint3 indexInSignal : SV_DispatchThreadID
) {
	const uint currentChannelsIndex = indexInSignal.z;

	FFTFromComplexTexture		(4, dontInvert, kernelsDFTRO, currentChannelsIndex, vertically, indexInSignal.x, indexInSignal.y);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.y,   4);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.y,  16);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.y,  64);
	TwoFFTsOnGroupSharedMemory	(2, dontInvert, indexInSignal.y, 256);
	FFTToComplexTexture			(4, dontInvert, kernelsDFTRW, currentChannelsIndex, vertically, indexInSignal.x, indexInSignal.y, 512);

	return;
}

[numthreads(IMAGE_SIDELENGTH / 4, 1, 1)] // 4 pixels per thread — best for radix 4
void imageHorizontal(
	const uint3 indexInSignal : SV_DispatchThreadID
) {
	// !ARCHITECTURE!
	//	• AMD Radeon HD: Performing the FFT with radix 4-2-4-4-4-4 saves one register.
	const uint currentChannelsIndex = indexInSignal.z;

	FFTFromRealTexture			(4, glareRO, currentChannelsIndex, horizontally, indexInSignal.y, indexInSignal.x);
	TwoFFTsOnGroupSharedMemory	(2, dontInvert, indexInSignal.x,   4);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.x,   8);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.x,  32);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.x, 128);
	FFTToComplexTexture			(4, dontInvert, imagesDFTRW, currentChannelsIndex, horizontally, indexInSignal.y, indexInSignal.x, 512);

	return;
}

[numthreads(1, IMAGE_SIDELENGTH / 4, 1)] // 4 pixels per thread — best for radix 4
void imageVertical(
	const uint3 indexInSignal : SV_DispatchThreadID
) {
	// !ARCHITECTURE!
	//	• AMD Radeon HD: Performing the FFT with radix 4-4-4-4-2-4 and the iFFT with radix 2-4-4-4-4-4 saves ten registers and
	//		37 scratch registers.
	const uint currentChannelsIndex = indexInSignal.z;

	FFTFromComplexTexture		(4, dontInvert, imagesDFTRO, currentChannelsIndex, vertically, indexInSignal.x, indexInSignal.y);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.y,   4);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.y,  16);
	FFTOnGroupSharedMemory		(4, dontInvert, indexInSignal.y,  64);
	TwoFFTsOnGroupSharedMemory	(2, dontInvert, indexInSignal.y, 256);
	FFTOnGroupSharedMemoryWithVerticalPSFMultiplication(4, currentChannelsIndex, indexInSignal.x, indexInSignal.y, 512);

	TwoFFTsOnGroupSharedMemory	(2, invert, indexInSignal.y,   1);
	FFTOnGroupSharedMemory		(4, invert, indexInSignal.y,   2);
	FFTOnGroupSharedMemory		(4, invert, indexInSignal.y,   8);
	FFTOnGroupSharedMemory		(4, invert, indexInSignal.y,  32);
	FFTOnGroupSharedMemory		(4, invert, indexInSignal.y, 128);
	FFTToComplexTexture			(4, invert, imagesDFTRW, currentChannelsIndex, vertically, indexInSignal.x, indexInSignal.y, 512);

	return;
}

[numthreads(IMAGE_SIDELENGTH / 4, 1, 1)] // 4 pixels per thread — best for radix 4
void imageInverseHorizontal(
	const uint3 indexInSignal : SV_DispatchThreadID
) {
	// !ARCHITECTURE!
	//	• AMD Radeon HD: Performing the FFT with radix 4-2-4-4-4-4 saves one register.
	const uint currentChannelsIndex = indexInSignal.z;

	FFTFromComplexTexture		(4, invert, imagesDFTRO, currentChannelsIndex, horizontally, indexInSignal.y, indexInSignal.x);
	TwoFFTsOnGroupSharedMemory	(2, invert, indexInSignal.x,   4);
	FFTOnGroupSharedMemory		(4, invert, indexInSignal.x,   8);
	FFTOnGroupSharedMemory		(4, invert, indexInSignal.x,  32);
	FFTOnGroupSharedMemory		(4, invert, indexInSignal.x, 128);
	iFFTToRealTexture			(4, glareRW, currentChannelsIndex, horizontally, indexInSignal.y, indexInSignal.x, 512);

	return;
}