Skip to content

Commit e4da85b

Browse files
krzysgkrzysg
committed
Gradient tests added
1 parent c365186 commit e4da85b

File tree

3 files changed

+184
-55
lines changed

3 files changed

+184
-55
lines changed

‎src/algorithm/ComputeGradient.hpp‎

Lines changed: 48 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -51,9 +51,9 @@ class ComputeGradient {
5151

5252
// Gradient computation
5353

54-
template<typename T, typename S>
54+
template<typename S>
5555
void
56-
calc_bspline_fd_ds_mag(MeshData<T> &input, MeshData<S> &grad, const float hx, const float hy, const float hz);
56+
calc_bspline_fd_ds_mag(const MeshData<S> &input, MeshData<S> &grad, const float hx, const float hy, const float hz);
5757

5858
template<typename T,typename S>
5959
void mask_gradient(MeshData<T>& grad_ds,MeshData<S>& temp_ds,MeshData<T>& temp_full,APRParameters& par);
@@ -721,13 +721,16 @@ void ComputeGradient::get_smooth_bspline_3D(MeshData<T>& input,APRParameters& pa
721721
spline_timer.stop_timer();
722722
}
723723

724-
template<typename T,typename S>
725-
void ComputeGradient::calc_bspline_fd_ds_mag(MeshData<T> &input, MeshData<S> &grad, const float hx, const float hy,const float hz) {
726-
//
727-
// Bevan Cheeseman 2016
728-
//
729-
// Calculate fd filt, for xgrad with bsplines
730-
724+
/**
725+
* Calculates downsampled gradient (maximum magnitude) with 'replicate' boundary approach (nearest border value)
726+
* @param input - input mesh
727+
* @param grad - output gradient (must be initialized)
728+
* @param hx - step in x dir
729+
* @param hy - step in y dir
730+
* @param hz - step in z dir
731+
*/
732+
template<typename S>
733+
void ComputeGradient::calc_bspline_fd_ds_mag(const MeshData<S> &input, MeshData<S> &grad, const float hx, const float hy,const float hz) {
731734
const size_t z_num = input.z_num;
732735
const size_t x_num = input.x_num;
733736
const size_t y_num = input.y_num;
@@ -738,67 +741,65 @@ void ComputeGradient::calc_bspline_fd_ds_mag(MeshData<T> &input, MeshData<S> &gr
738741
std::vector<S> temp(y_num, 0);
739742
const size_t xnumynum = x_num * y_num;
740743

741-
// 4 4
742-
// 2 1,3 ... 1 2 3 ...
743-
// 5 5
744744
#ifdef HAVE_OPENMP
745-
#pragma omp parallel for default(shared) firstprivate(temp)
745+
#pragma omp parallel for default(shared) firstprivate(temp)
746746
#endif
747-
for (size_t j = 0; j < z_num; ++j) {
748-
S *left = input.mesh.begin() + j * xnumynum + 1 * y_num;
749-
S *center = input.mesh.begin() + j * xnumynum;
747+
for (size_t z = 0; z < z_num; ++z) {
748+
// Belows pointers up, down... are forming stencil in X (left <-> right) and Z ( up <-> down) direction and
749+
// are pointing to whole Y column. If out of bounds then 'replicate' (nearest array border value) approach is used.
750+
//
751+
// up
752+
// ... left center right ...
753+
// down
754+
const S *left = input.mesh.begin() + z * xnumynum + 0 * y_num; // boundary value is chosen
755+
const S *center = input.mesh.begin() + z * xnumynum + 0 * y_num;
750756

751757
//LHS boundary condition is accounted for wiht this initialization
752-
const size_t j_m = j > 0 ? j - 1 : 0;
753-
const size_t j_p = std::min(z_num - 1, j + 1);
754-
float g[x_num][y_num];
755-
for (size_t i = 0; i < x_num - 1; ++i) {
756-
S *up = input.mesh.begin() + j_m * xnumynum + i * y_num;
757-
S *down = input.mesh.begin() + j_p * xnumynum + i * y_num;
758-
S *right = input.mesh.begin() + j * xnumynum + (i + 1) * y_num;
758+
const size_t zMinus = z > 0 ? z - 1 : 0 /* boundary */;
759+
const size_t zPlus = std::min(z + 1, z_num - 1 /* boundary */);
760+
761+
for (size_t x = 0; x < x_num; ++x) {
762+
const S *up = input.mesh.begin() + zMinus * xnumynum + x * y_num;
763+
const S *down = input.mesh.begin() + zPlus * xnumynum + x * y_num;
764+
const size_t xPlus = std::min(x + 1, x_num - 1 /* boundary */);
765+
const S *right = input.mesh.begin() + z * xnumynum + xPlus * y_num;
759766

760767
//compute the boundary values
761-
temp[0] = sqrt(pow((right[0] - left[0]) / (2 * hx), 2.0) + pow((down[0] - up[0]) / (2 * hz), 2.0));
762-
g[0][i] = temp[0];
763-
//do the y gradient
764-
#ifdef HAVE_OPENMP
765-
#pragma omp simd
766-
#endif
767-
for (size_t k = 1; k < y_num - 1; ++k) {
768-
temp[k] = sqrt(pow((right[k] - left[k]) / (2 * hx), 2.0) + pow((down[k] - up[k]) / (2 * hz), 2.0) +
769-
pow((center[k + 1] - center[k - 1]) / (2 * hy), 2.0));
770-
g[k][i] = temp[k];
768+
if (y_num >= 2) {
769+
temp[0] = sqrt(pow((right[0] - left[0]) / (2 * hx), 2.0) + pow((down[0] - up[0]) / (2 * hz), 2.0) + pow((center[1] - center[0 /* boundary */]) / (2 * hy), 2.0));
770+
temp[y_num - 1] = sqrt(pow((right[y_num - 1] - left[y_num - 1]) / (2 * hx), 2.0) + pow((down[y_num - 1] - up[y_num - 1]) / (2 * hz), 2.0) + pow((center[y_num - 1 /* boundary */] - center[y_num - 2]) / (2 * hy), 2.0));
771+
}
772+
else {
773+
temp[0] = 0; // same values minus same values in x/y/z
771774
}
772775

773-
temp[y_num - 1] = sqrt(pow((right[y_num - 1] - left[y_num - 1]) / (2 * hx), 2.0) +
774-
pow((down[y_num - 1] - up[y_num - 1]) / (2 * hz), 2.0));
775-
g[y_num - 1][i] = temp[y_num - 1];
776+
//do the y gradient in range 1..y_num-2
777+
#ifdef HAVE_OPENMP
778+
#pragma omp simd
779+
#endif
780+
for (size_t y = 1; y < y_num - 1; ++y) {
781+
temp[y] = sqrt(pow((right[y] - left[y]) / (2 * hx), 2.0) + pow((down[y] - up[y]) / (2 * hz), 2.0) + pow((center[y + 1] - center[y - 1]) / (2 * hy), 2.0));
782+
}
776783

777-
int64_t j_2 = j / 2;
778-
int64_t i_2 = i / 2;
784+
// Set as a downsampled gradient maximum from 2x2x2 gradient cubes
785+
int64_t z_2 = z / 2;
786+
int64_t x_2 = x / 2;
779787
for (size_t k = 0; k < y_num_ds; ++k) {
780788
size_t k_s = std::min(2 * k + 1, y_num - 1);
781-
const size_t idx = j_2 * x_num_ds * y_num_ds + i_2 * y_num_ds + k;
789+
const size_t idx = z_2 * x_num_ds * y_num_ds + x_2 * y_num_ds + k;
782790
grad.mesh[idx] = std::max(temp[2 * k], std::max(temp[k_s], grad.mesh[idx]));
783791
}
784792

785793
// move left, center to current center, right (both +1 to right)
786794
std::swap(left, center);
787795
std::swap(center, right);
788796
}
789-
for (int y = 0; y < y_num; ++y) {
790-
for (int x = 0; x < x_num; ++x) {
791-
std::cout << g[y][x] << " ";
792-
}
793-
std::cout << "\n";
794-
}
795797
}
796798
}
797799

798-
/*
800+
/**
799801
* Caclulation of signal value from B-Spline co-efficients
800802
*/
801-
802803
template<typename T>
803804
void ComputeGradient::calc_inv_bspline_y(MeshData<T>& input){
804805
//

‎src/data_structures/Mesh/MeshData.hpp‎

Lines changed: 24 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@
1717
#include <cmath>
1818
#include <memory>
1919
#include <sstream>
20+
#include <iostream>
21+
#include <iomanip>
2022

2123
#include "src/misc/APRTimer.hpp"
2224

@@ -396,18 +398,32 @@ public :
396398
}
397399

398400
/**
399-
* Prints X-Y planes of mesh (for debug/test purposses - use only on small meshes)
401+
* Prints X-Y or X-Z planes of mesh (for debug/test purposses - use only on small meshes)
400402
*/
401-
void printMesh() const {
402-
for (size_t z = 0; z < z_num; ++z) {
403-
std::cout << "z=" << z << "\n";
403+
void printMesh(int aColumnWidth, bool aXYplanes = true) const {
404+
if (aXYplanes) {
405+
for (size_t z = 0; z < z_num; ++z) {
406+
std::cout << "z=" << z << "\n";
407+
for (size_t y = 0; y < y_num; ++y) {
408+
for (size_t x = 0; x < x_num; ++x) {
409+
std::cout << std::setw(aColumnWidth) << at(y, x, z) << " ";
410+
}
411+
std::cout << "\n";
412+
}
413+
std::cout << std::endl;
414+
}
415+
}
416+
else { // X-Z planes
404417
for (size_t y = 0; y < y_num; ++y) {
405-
for (size_t x = 0; x < x_num; ++x) {
406-
std::cout << at(y, x, z) << " ";
418+
std::cout << "y=" << y << "\n";
419+
for (size_t z = 0; z < z_num; ++z) {
420+
for (size_t x = 0; x < x_num; ++x) {
421+
std::cout << std::setw(aColumnWidth) << at(y, x, z) << " ";
422+
}
423+
std::cout << "\n";
407424
}
408-
std::cout << "\n";
425+
std::cout << std::endl;
409426
}
410-
std::cout << std::endl;
411427
}
412428
}
413429

‎test/ComputeGradientTest.cpp‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
/*
2+
* Created by Krzysztof Gonciarz 2018
3+
*/
4+
5+
#include <gtest/gtest.h>
6+
#include "src/data_structures/Mesh/MeshData.hpp"
7+
#include "src/algorithm/ComputeGradient.hpp"
8+
9+
namespace {
10+
/**
11+
* Compares mesh with provided data
12+
* @param mesh
13+
* @param data - data with [Z][Y][X] structure
14+
* @return true if same
15+
*/
16+
template <typename T>
17+
bool compare(MeshData<T> &mesh, const float *data, const float epsilon) {
18+
size_t dataIdx = 0;
19+
for (size_t z = 0; z < mesh.z_num; ++z) {
20+
for (size_t y = 0; y < mesh.y_num; ++y) {
21+
for (size_t x = 0; x < mesh.x_num; ++x) {
22+
bool v = std::abs(mesh(y, x, z) - data[dataIdx]) < epsilon;
23+
if (v == false) {
24+
std::cerr << "Mesh and expected data differ. First place at (Y, X, Z) = " << y << ", " << x << ", " << z << ") " << mesh(y, x, z) << " vs " << data[dataIdx] << std::endl;
25+
return false;
26+
}
27+
++dataIdx;
28+
}
29+
}
30+
}
31+
return true;
32+
}
33+
34+
TEST(ComputeGradientTest, 2D_XY) {
35+
{ // Corner points
36+
MeshData<float> m(6, 6, 1, 0);
37+
// expect gradient is 3x3 X/Y plane
38+
float expect[] = {1.41, 0, 4.24,
39+
0, 0, 0,
40+
2.82, 0, 5.65};
41+
// put values in corners
42+
m(0, 0, 0) = 2;
43+
m(5, 0, 0) = 4;
44+
m(0, 5, 0) = 6;
45+
m(5, 5, 0) = 8;
46+
MeshData<float> grad;
47+
grad.initDownsampled(m, 0);
48+
ComputeGradient cg;
49+
cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
50+
ASSERT_TRUE(compare(grad, expect, 0.01));
51+
}
52+
{ // In the middle
53+
MeshData<float> m(6, 6, 1, 0);
54+
// expect gradient is 3x3 X/Y plane
55+
float expect[] = {1, 1, 0,
56+
1, 0, 0,
57+
0, 0, 0};
58+
// put values in corners
59+
m(1, 1, 0) = 2;
60+
MeshData<float> grad;
61+
grad.initDownsampled(m, 0);
62+
ComputeGradient cg;
63+
cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
64+
ASSERT_TRUE(compare(grad, expect, 0.01));
65+
}
66+
{ // One pixel image 1x1x1
67+
MeshData<float> m(1, 1, 1, 0);
68+
// expect gradient is 3x3 X/Y plane
69+
float expect[] = {0};
70+
// put values in corners
71+
m(0, 0, 0) = 2;
72+
MeshData<float> grad;
73+
grad.initDownsampled(m, 0);
74+
ComputeGradient cg;
75+
cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
76+
ASSERT_TRUE(compare(grad, expect, 0.01));
77+
}
78+
79+
}
80+
81+
TEST(ComputeGradientTest, Corners3D) {
82+
MeshData<float> m(6, 6, 4, 0);
83+
// expect gradient is 3x3x2 X/Y/Z plane
84+
float expect[] = { 1.73, 0, 5.19,
85+
0, 0, 0,
86+
3.46, 0, 6.92,
87+
88+
8.66, 0, 12.12,
89+
0, 0, 0,
90+
10.39, 0, 13.85 };
91+
// put values in corners
92+
m(0, 0, 0) = 2;
93+
m(5, 0, 0) = 4;
94+
m(0, 5, 0) = 6;
95+
m(5, 5, 0) = 8;
96+
m(0, 0, 3) = 10;
97+
m(5, 0, 3) = 12;
98+
m(0, 5, 3) = 14;
99+
m(5, 5, 3) = 16;
100+
101+
MeshData<float> grad;
102+
grad.initDownsampled(m, 0);
103+
ComputeGradient cg;
104+
cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
105+
ASSERT_TRUE(compare(grad, expect, 0.01));
106+
}
107+
}
108+
109+
int main(int argc, char **argv) {
110+
testing::InitGoogleTest(&argc, argv);
111+
return RUN_ALL_TESTS();
112+
}

0 commit comments

Comments
 (0)