Source browser
examples/spectral.cpp
Read the complete file below or download its original bytes. Line numbers provide stable links to a declaration or processing step.
File fingerprint
1/** @file spectral.cpp
2 * Linked 2D/3D convolution and Fourier pipelines over seed-reproducible data.
3 * Filters use true convolution with clamp-to-edge. The 2D direct DFT is a small
4 * direct transform. An independent full 2D reference
5 * checks the separable GPU/CPU stages and the inverse reconstruction.
6 */
7#include "runner.hpp" // Include runner.hpp declarations.
8#include <algorithm> // Include algorithm declarations.
9#include <complex> // Include complex declarations.
10#include <random> // Include random declarations.
11#include <vector> // Include vector declarations.
12
13namespace { // Group these helpers within their namespace.
14typedef shdr::glsl::vec2 sample; // Declare shdr for the following calculation.
15typedef shdr::glsl::uvec2 size2; // Declare shdr for the following calculation.
16typedef shdr::glsl::uvec3 size3; // Declare shdr for the following calculation.
17
18int edge(int coordinate, unsigned size) { // Implement edge for the operations below.
19 return std::max(0, std::min(coordinate, int(size) - 1)); // Return the result of this operation.
20} // End this block.
21
22// Accumulate the full reference products in double precision.
23std::vector<sample> reference2d(const std::vector<sample>& input, // Declare reference2d with the parameters continued below.
24 size2 size, // Declare size for the following calculation.
25 const std::vector<float>& weights, // Declare weights for the following calculation.
26 size2 kernel_size) { // Declare kernel size for the following calculation.
27 std::vector<sample> output(input.size()); // Construct output with these initial values.
28 for (unsigned y = 0; y < size.y; ++y) // Iterate over y for this calculation.
29 for (unsigned x = 0; x < size.x; ++x) { // Iterate over x for this calculation.
30 double real = 0, imaginary = 0; // Use this floating-point value or resource in the following calculation.
31 for (unsigned ky = 0; ky < kernel_size.y; ++ky) // Iterate over ky for this calculation.
32 for (unsigned kx = 0; kx < kernel_size.x; ++kx) { // Iterate over kx for this calculation.
33 int qx = edge(int(x) + int(kernel_size.x / 2) - int(kx), size.x); // Set qx from this expression.
34 int qy = edge(int(y) + int(kernel_size.y / 2) - int(ky), size.y); // Set qy from this expression.
35 const sample value = input[unsigned(qx) + size.x * unsigned(qy)]; // Set value from this expression.
36 const double coefficient = weights[kx + kernel_size.x * ky]; // Set coefficient from this expression.
37 real += double(value.x) * coefficient; // Update real with this contribution.
38 imaginary += double(value.y) * coefficient; // Update imaginary with this contribution.
39 } // End this block.
40 output[x + size.x * y] = sample(float(real), float(imaginary)); // Set output[x + size.x * y] from this expression.
41 } // End this block.
42 return output; // Return output.
43} // End this block.
44
45std::vector<float> reference3d(const std::vector<float>& input, // Declare reference3d with the parameters continued below.
46 size3 size, // Declare size for the following calculation.
47 const std::vector<float>& weights, // Declare weights for the following calculation.
48 size3 kernel_size) { // Declare kernel size for the following calculation.
49 std::vector<float> output(input.size()); // Construct output with these initial values.
50 for (unsigned z = 0; z < size.z; ++z) // Iterate over z for this calculation.
51 for (unsigned y = 0; y < size.y; ++y) // Iterate over y for this calculation.
52 for (unsigned x = 0; x < size.x; ++x) { // Iterate over x for this calculation.
53 double total = 0; // Use this floating-point value or resource in the following calculation.
54 for (unsigned kz = 0; kz < kernel_size.z; ++kz) // Iterate over kz for this calculation.
55 for (unsigned ky = 0; ky < kernel_size.y; ++ky) // Iterate over ky for this calculation.
56 for (unsigned kx = 0; kx < kernel_size.x; ++kx) { // Iterate over kx for this calculation.
57 unsigned qx = // Set qx from this expression.
58 unsigned(edge(int(x) + int(kernel_size.x / 2) - int(kx), size.x)); // Apply unsigned to these arguments.
59 unsigned qy = // Set qy from this expression.
60 unsigned(edge(int(y) + int(kernel_size.y / 2) - int(ky), size.y)); // Apply unsigned to these arguments.
61 unsigned qz = // Set qz from this expression.
62 unsigned(edge(int(z) + int(kernel_size.z / 2) - int(kz), size.z)); // Apply unsigned to these arguments.
63 total += // Update total with this contribution.
64 double(input[qx + size.x * (qy + size.y * qz)]) * // Use this floating-point value or resource in the following calculation.
65 double(weights[kx + kernel_size.x * (ky + kernel_size.y * kz)]); // Use this floating-point value or resource in the following calculation.
66 } // End this block.
67 output[x + size.x * (y + size.y * z)] = float(total); // Set output[x + size.x * (y + size.y * z)] from this expression.
68 } // End this block.
69 return output; // Return output.
70} // End this block.
71
72// Evaluate the full two-dimensional DFT directly as an independent reference.
73std::vector<sample> reference_dft(const std::vector<sample>& input, size2 size) { // Implement reference dft for the operations below.
74 std::vector<sample> output(input.size()); // Construct output with these initial values.
75 const double tau = 6.283185307179586476925286766559; // Set tau from this expression.
76 for (unsigned ky = 0; ky < size.y; ++ky) // Iterate over ky for this calculation.
77 for (unsigned kx = 0; kx < size.x; ++kx) { // Iterate over kx for this calculation.
78 std::complex<double> total(0, 0); // Construct total with these initial values.
79 for (unsigned y = 0; y < size.y; ++y) // Iterate over y for this calculation.
80 for (unsigned x = 0; x < size.x; ++x) { // Iterate over x for this calculation.
81 sample value = input[x + size.x * y]; // Set value from this expression.
82 double angle = -tau * (double(kx) * x / size.x + double(ky) * y / size.y); // Use this floating-point value or resource in the following calculation.
83 total += std::complex<double>(value.x, value.y) * std::polar(1.0, angle); // Update total with this contribution.
84 } // End this block.
85 output[kx + size.x * ky] = sample(float(total.real()), float(total.imag())); // Set output[kx + size.x * ky] from this expression.
86 } // End this block.
87 return output; // Return output.
88} // End this block.
89
90template <class T> shdr::buffer<T> resource(const std::vector<T>& values) { // Accept the value types or dimensions used by this template.
91 shdr::buffer<T> result(values.size()); // Construct result with these initial values.
92 std::copy(values.begin(), values.end(), result.data()); // Apply copy to these arguments.
93 return result; // Return result.
94} // End this block.
95
96void check(const shdr::buffer<sample>& actual, // Declare check with the parameters continued below.
97 const std::vector<sample>& expected, // Declare expected for the following calculation.
98 const std::string& label, // Declare label for the following calculation.
99 float tolerance = 0.0002f) { // Set tolerance from this expression.
100 for (std::size_t i = 0; i < expected.size(); ++i) { // Iterate over i for this calculation.
101 shdr_examples::require_close(actual.data()[i].x, // Check the computed value against the independent reference.
102 expected[i].x, // Continue the expression using expected, i, x.
103 label + " real[" + std::to_string(i) + "]", // Apply to string to these arguments.
104 tolerance); // Continue the expression using tolerance.
105 shdr_examples::require_close(actual.data()[i].y, // Check the computed value against the independent reference.
106 expected[i].y, // Continue the expression using expected, i, y.
107 label + " imaginary[" + std::to_string(i) + "]", // Apply to string to these arguments.
108 tolerance); // Continue the expression using tolerance.
109 } // End this block.
110} // End this block.
111void check(const shdr::buffer<float>& actual, // Declare check with the parameters continued below.
112 const std::vector<float>& expected, // Declare expected for the following calculation.
113 const std::string& label) { // Declare label for the following calculation.
114 for (std::size_t i = 0; i < expected.size(); ++i) // Iterate over i for this calculation.
115 shdr_examples::require_close( // Check the computed value against the independent reference.
116 actual.data()[i], expected[i], label + "[" + std::to_string(i) + "]", 0.0002f); // Apply data to these arguments.
117} // End this block.
118
119std::unique_ptr<shdr::kernel> convolution2d(const shdr::buffer<sample>& input, // Declare convolution2d with the parameters continued below.
120 const shdr::buffer<sample>& output, // Declare output for the following calculation.
121 const shdr::buffer<float>& weights, // Declare weights for the following calculation.
122 size2 dimensions, // Declare dimensions for the following calculation.
123 size2 kernel_size) { // Declare kernel size for the following calculation.
124 std::unique_ptr<shdr::kernel> shader = shdr::make_kernel("convolution2d_kernel"); // Create the linked convolution2d kernel instance.
125 shader->bind("input_values", input).bind("output_values", output).bind("coefficients", weights); // Connect this shader port to its application resource.
126 shader->set<size2>("dimensions", dimensions).set<size2>("kernel_size", kernel_size); // Supply this shader uniform with its typed application value.
127 return shader; // Return shader.
128} // End this block.
129std::unique_ptr<shdr::kernel> convolution3d(const shdr::buffer<float>& input, // Declare convolution3d with the parameters continued below.
130 const shdr::buffer<float>& output, // Declare output for the following calculation.
131 const shdr::buffer<float>& weights, // Declare weights for the following calculation.
132 size3 dimensions, // Declare dimensions for the following calculation.
133 size3 kernel_size) { // Declare kernel size for the following calculation.
134 std::unique_ptr<shdr::kernel> shader = shdr::make_kernel("convolution3d_kernel"); // Create the linked convolution3d kernel instance.
135 shader->bind("input_values", input).bind("output_values", output).bind("coefficients", weights); // Connect this shader port to its application resource.
136 shader->set<size3>("dimensions", dimensions).set<size3>("kernel_size", kernel_size); // Supply this shader uniform with its typed application value.
137 return shader; // Return shader.
138} // End this block.
139std::unique_ptr<shdr::kernel> fourier(const char* name, // Declare fourier with the parameters continued below.
140 const shdr::buffer<sample>& input, // Declare input for the following calculation.
141 const shdr::buffer<sample>& output, // Declare output for the following calculation.
142 size2 dimensions, // Declare dimensions for the following calculation.
143 float direction, // Declare direction for the following calculation.
144 float normalization) { // Declare normalization for the following calculation.
145 std::unique_ptr<shdr::kernel> shader = shdr::make_kernel(name); // Create the linked shader instance.
146 shader->bind("input_values", input).bind("output_values", output); // Connect this shader port to its application resource.
147 shader->set<size2>("dimensions", dimensions) // Supply this shader uniform with its typed application value.
148 .set<float>("direction", direction) // Supply this shader uniform with its typed application value.
149 .set<float>("normalization", normalization); // Supply this shader uniform with its typed application value.
150 return shader; // Return shader.
151} // End this block.
152
153void image_and_fourier(shdr_examples::runner& run, std::mt19937& random) { // Implement image and fourier for the operations below.
154 const size2 dimensions(8, 8); // Construct dimensions with these initial values.
155 const unsigned count = dimensions.x * dimensions.y; // Set count from this expression.
156 std::uniform_real_distribution<float> distribution(-1, 1); // Construct distribution with these initial values.
157 std::vector<sample> pixels(count); // Construct pixels with these initial values.
158 for (unsigned i = 0; i < count; ++i) // Iterate over i for this calculation.
159 pixels[i] = sample(distribution(random), 0); // Set pixels[i] from this expression.
160 pixels[0].x += 3; // Edge and interior impulses make clamp behavior visible.
161 pixels[3 + dimensions.x * 4].x += 5; // Update pixels[3 + dimensions.x * 4].x with this contribution.
162 const std::vector<float> blur_coefficients = {0, 0.125f, 0, 0.125f, 0.5f, 0.125f, 0, 0.125f, 0}; // Set blur coefficients from this expression.
163 const std::vector<float> derivative_coefficients = {-0.5f, 0, 0.5f}; // Set derivative coefficients from this expression.
164 shdr::buffer<sample> input = resource(pixels), blurred(count), filtered(count), // Set input from this expression.
165 row_frequency(count), spectrum(count), inverse_rows(count), // Apply row frequency to these arguments.
166 reconstructed(count); // Apply reconstructed to these arguments.
167 shdr::buffer<float> blur_weights = resource(blur_coefficients), // Set blur weights from this expression.
168 derivative_weights = resource(derivative_coefficients); // Set derivative weights from this expression.
169 auto blur = convolution2d(input, blurred, blur_weights, dimensions, size2(3, 3)); // Set blur from this expression.
170 auto derivative = convolution2d(blurred, filtered, derivative_weights, dimensions, size2(3, 1)); // Set derivative from this expression.
171 auto rows = fourier("dft_rows_kernel", filtered, row_frequency, dimensions, -1, 1); // Set rows from this expression.
172 auto columns = fourier("dft_columns_kernel", row_frequency, spectrum, dimensions, -1, 1); // Set columns from this expression.
173 auto inverse_row = // Set inverse row from this expression.
174 fourier("dft_rows_kernel", spectrum, inverse_rows, dimensions, 1, 1.0f / dimensions.x); // Apply fourier to these arguments.
175 auto inverse_column = fourier( // Set inverse column from this expression.
176 "dft_columns_kernel", inverse_rows, reconstructed, dimensions, 1, 1.0f / dimensions.y); // Provide this text argument for the operation above.
177 const shdr::extent groups((dimensions.x + 7) / 8, (dimensions.y + 7) / 8); // Construct groups with these initial values.
178 shdr::pipeline pipeline = // Set pipeline from this expression.
179 shdr::dispatch(*blur, groups) >> shdr::dispatch(*derivative, groups) >> // Add a shader dispatch with the required workgroup counts.
180 shdr::dispatch(*rows, groups) >> shdr::dispatch(*columns, groups) >> // Add a shader dispatch with the required workgroup counts.
181 shdr::dispatch(*inverse_row, groups) >> shdr::dispatch(*inverse_column, groups); // Add a shader dispatch with the required workgroup counts.
182 run.run(pipeline); // Execute the configured pipeline on the requested backend.
183 const std::vector<sample> reference_blur = // Set reference blur from this expression.
184 reference2d(pixels, dimensions, blur_coefficients, size2(3, 3)); // Apply reference2d to these arguments.
185 const std::vector<sample> reference_filter = // Set reference filter from this expression.
186 reference2d(reference_blur, dimensions, derivative_coefficients, size2(3, 1)); // Apply reference2d to these arguments.
187 check(blurred, reference_blur, "2D convolution blur"); // Apply check to these arguments.
188 check(filtered, reference_filter, "2D asymmetric convolution"); // Apply check to these arguments.
189 check(spectrum, reference_dft(reference_filter, dimensions), "2D forward DFT", 0.002f); // Apply check to these arguments.
190 check(reconstructed, reference_filter, "2D inverse DFT reconstruction", 0.0003f); // Apply check to these arguments.
191 // Kernel resources are replaceable after composing the pipeline.
192 const std::vector<float> identity_coefficients = {1}; // Set identity coefficients from this expression.
193 shdr::buffer<float> identity_weights = resource(identity_coefficients); // Set identity weights from this expression.
194 derivative->bind("coefficients", identity_weights).set<size2>("kernel_size", size2(1)); // Connect this shader port to its application resource.
195 run.run(pipeline); // Execute the configured pipeline on the requested backend.
196 check(filtered, reference_blur, "rebound convolution coefficients"); // Apply check to these arguments.
197 check(reconstructed, reference_blur, "rebound inverse reconstruction", 0.0003f); // Apply check to these arguments.
198 std::cout << "8x8 2D convolution -> convolution -> DFT rows/columns -> inverse: verified\n"; // Write this part of the example result or diagnostic.
199} // End this block.
200
201void volume(shdr_examples::runner& run, std::mt19937& random) { // Implement volume for the operations below.
202 const size3 dimensions(5, 4, 3); // Construct dimensions with these initial values.
203 const unsigned count = dimensions.x * dimensions.y * dimensions.z; // Set count from this expression.
204 std::uniform_real_distribution<float> distribution(-1, 1); // Construct distribution with these initial values.
205 std::vector<float> voxels(count); // Construct voxels with these initial values.
206 for (unsigned i = 0; i < count; ++i) // Iterate over i for this calculation.
207 voxels[i] = distribution(random); // Set voxels[i] from this expression.
208 voxels[0] += 2; // Update voxels[0] with this contribution.
209 voxels[2 + dimensions.x * (1 + dimensions.y)] += 4; // Update voxels[2 + dimensions.x * (1 + dimensions.y)] with this contribution.
210 const std::vector<float> x_coefficients = {0.25f, 0.5f, 0.25f}; // Set x coefficients from this expression.
211 const std::vector<float> y_coefficients = {-0.25f, 0.75f, 0.5f}; // Set y coefficients from this expression.
212 const std::vector<float> z_coefficients = {0.2f, 0.3f, 0.5f}; // Set z coefficients from this expression.
213 shdr::buffer<float> input = resource(voxels), along_x(count), along_y(count), output(count); // Set input from this expression.
214 shdr::buffer<float> x_weights = resource(x_coefficients), y_weights = resource(y_coefficients), // Set x weights from this expression.
215 z_weights = resource(z_coefficients); // Set z weights from this expression.
216 auto x = convolution3d(input, along_x, x_weights, dimensions, size3(3, 1, 1)); // Set x from this expression.
217 auto y = convolution3d(along_x, along_y, y_weights, dimensions, size3(1, 3, 1)); // Set y from this expression.
218 auto z = convolution3d(along_y, output, z_weights, dimensions, size3(1, 1, 3)); // Set z from this expression.
219 const shdr::extent groups( // Declare groups with the parameters continued below.
220 (dimensions.x + 3) / 4, (dimensions.y + 3) / 4, (dimensions.z + 3) / 4); // Continue the expression using dimensions, x, dimensions.
221 shdr::pipeline pipeline = // Set pipeline from this expression.
222 shdr::dispatch(*x, groups) >> shdr::dispatch(*y, groups) >> shdr::dispatch(*z, groups); // Add a shader dispatch with the required workgroup counts.
223 run.run(pipeline); // Execute the configured pipeline on the requested backend.
224 const std::vector<float> reference_x = // Set reference x from this expression.
225 reference3d(voxels, dimensions, x_coefficients, size3(3, 1, 1)); // Apply reference3d to these arguments.
226 const std::vector<float> reference_y = // Set reference y from this expression.
227 reference3d(reference_x, dimensions, y_coefficients, size3(1, 3, 1)); // Apply reference3d to these arguments.
228 const std::vector<float> reference_z = // Set reference z from this expression.
229 reference3d(reference_y, dimensions, z_coefficients, size3(1, 1, 3)); // Apply reference3d to these arguments.
230 check(along_x, reference_x, "3D X convolution"); // Apply check to these arguments.
231 check(along_y, reference_y, "3D Y convolution"); // Apply check to these arguments.
232 check(output, reference_z, "3D Z convolution"); // Apply check to these arguments.
233 std::cout << "5x4x3 3D convolution X -> Y -> Z: verified\n"; // Write this part of the example result or diagnostic.
234} // End this block.
235} // namespace
236
237int main(int argc, char** argv) { // Run this example and validate its results.
238 try { // Run the operation while retaining the error handling below.
239 shdr_examples::runner run(argc, argv); // Construct run with these initial values.
240 std::mt19937 random(0x53484452u); // Construct random with these initial values.
241 image_and_fourier(run, random); // Apply image and fourier to these arguments.
242 volume(run, random); // Apply volume to these arguments.
243 std::cout << (run.gpu() ? "GPU" : "CPU") << " spectral pipelines passed\n"; // Write this part of the example result or diagnostic.
244 return 0; // Report successful completion.
245#ifdef SHDR_HAVE_EGL // Enable this branch when its feature macro is defined.
246 } catch (const shdr_examples::context_unavailable& error) { // Report the exception from the preceding operation.
247 std::cerr << "SKIP: " << error.what() << '\n'; // Write this part of the example result or diagnostic.
248 return 77; // Report an unavailable optional GPU backend.
249#endif // Finish this conditional compilation branch.
250 } catch (const std::exception& error) { // Report the exception from the preceding operation.
251 std::cerr << error.what() << '\n'; // Write this part of the example result or diagnostic.
252 return 1; // Report a failed check or operation.
253 } // End this block.
254} // End this block.