PURIFY
Next-generation radio interferometric imaging
Functions
convolution.cc File Reference
#include "purify/config.h"
#include "purify/types.h"
#include <iostream>
#include "catch2/catch_all.hpp"
#include "purify/directories.h"
#include "purify/logging.h"
#include "purify/convolution.h"
+ Include dependency graph for convolution.cc:

Go to the source code of this file.

Functions

 TEST_CASE ("1d_zeropad")
 
 TEST_CASE ("1d_convolution")
 
 TEST_CASE ("2d_convolution")
 
 TEST_CASE ("2d_convolution_functions")
 

Function Documentation

◆ TEST_CASE() [1/4]

TEST_CASE ( "1d_convolution"  )

Definition at line 23 of file convolution.cc.

23  {
24  Vector<t_int> tri = Vector<t_int>::Zero(3);
25  tri << 1, 2, 1;
26  Vector<t_int> signal = convol::zero_pad<t_int>(Vector<t_int>::Random(4), tri.size());
27  const Vector<t_int> output = convol::linear_convol_1d(tri, signal);
28  CHECK(output.size() == signal.size() - tri.size() + 1);
29  for (int i = 0; i < signal.size() - tri.size() - 1; i++)
30  CHECK(output(i) == tri(0) * signal(output.size() - i - 1) +
31  tri(1) * signal(output.size() - i - 1 + 1) +
32  tri(2) * signal(output.size() - i - 1 + 2));
33 }
#define CHECK(CONDITION, ERROR)
Definition: casa.cc:6
Vector< T > linear_convol_1d(const Vector< T > &kernelf, const Vector< T > &kernelg)
1d linear convoluiton of entire signal
Definition: convolution.h:43

References CHECK, and purify::convol::linear_convol_1d().

◆ TEST_CASE() [2/4]

TEST_CASE ( "1d_zeropad"  )

Definition at line 12 of file convolution.cc.

12  {
13  const Vector<t_int> signal = Vector<t_int>::Random(3);
14  for (int i = 1; i < 5; i++) {
15  const Vector<t_int> output = convol::zero_pad(signal, i);
16  CHECK(output.size() == signal.size() + 2 * i);
17  for (int j = 0; j < output.size(); j++) {
18  if (j < i) CHECK(output(j) == 0);
19  if (j > (i + signal.size())) CHECK(output(j) == 0);
20  }
21  }
22 }
Vector< T > zero_pad(const Vector< T > &input, const t_int padding)
zero pad a vector by a given amount
Definition: convolution.h:28

References CHECK, and purify::convol::zero_pad().

◆ TEST_CASE() [3/4]

TEST_CASE ( "2d_convolution"  )

Definition at line 34 of file convolution.cc.

34  {
35  const Vector<t_real> kernelu = Vector<t_real>::Random(3);
36  const Vector<t_real> kernelv = Vector<t_real>::Random(2);
37  const Matrix<t_real> signal = Matrix<t_real>::Random(2, 2);
38  const Matrix<t_real> output = convol::linear_convol_2d(kernelu, kernelv, signal);
39  CHECK(output.cols() == signal.cols() + kernelu.size() - 1);
40  CHECK(output.rows() == signal.rows() + kernelv.size() - 1);
41  const Matrix<t_real> test_buff =
42  convol::zero_pad<t_real>(signal, kernelv.size() - 1, kernelu.size() - 1);
43  // checking for signs that there is a boundary of zero, suggesting an offset
44  CHECK(output(0, 0) != output(1, 0));
45  CHECK(output(output.rows() - 1, output.cols() - 1) !=
46  output(output.rows() - 1, output.cols() - 2));
47  for (int i = 0; i < output.cols(); i++) {
48  for (int j = 0; j < output.rows(); j++) {
49  CHECK(std::abs(
50  output(j, i) -
51  kernelv(0) *
52  (kernelu(0) * test_buff(output.rows() - j - 1, output.cols() - i - 1) +
53  kernelu(1) * test_buff(output.rows() - j - 1, output.cols() - i + 1 - 1) +
54  kernelu(2) * test_buff(output.rows() - j - 1, output.cols() - i - 1 + 2)) -
55  kernelv(1) *
56  (kernelu(0) * test_buff(output.rows() - 1 - j + 1, output.cols() - i - 1) +
57  kernelu(1) * test_buff(output.rows() - 1 - j + 1, output.cols() - i - 1 + 1) +
58  kernelu(2) *
59  test_buff(output.rows() - 1 - j + 1, output.cols() - i - 1 + 2))) < 1e-7);
60  }
61  }
62 }
Matrix< T > linear_convol_2d(const Vector< T > &kernelfu, const Vector< T > &kernelfv, const Matrix< T > &kernelg)
perform linear convolution between two separable kernels and a 2d kernel (vectors)
Definition: convolution.h:55

References CHECK, and purify::convol::linear_convol_2d().

◆ TEST_CASE() [4/4]

TEST_CASE ( "2d_convolution_functions"  )

Definition at line 63 of file convolution.cc.

63  {
64  const t_int Ju = 3;
65  const t_int Jv = 4;
66  SECTION("box and delta") {
67  const t_int Jw = 5;
68  auto convol_kernelu = [=](const t_int ju) -> t_complex { return 1.; };
69  auto convol_kernelv = [=](const t_int jv) -> t_complex { return 1.; };
70  auto convol_kernelw = [=](const t_int jv, const t_int ju) -> t_complex {
71  return ((std::floor(ju - (Jw - 1) * 0.5) == 0) and (std::floor(jv - (Jw - 1) * 0.5) == 0))
72  ? 1.
73  : 0.;
74  };
75  const Matrix<t_complex> projection_kernel = convol::linear_convol_2d<t_complex>(
76  convol_kernelu, convol_kernelv, convol_kernelw, static_cast<t_int>(Ju),
77  static_cast<t_int>(Jv), static_cast<t_int>(Jw), static_cast<t_int>(Jw));
78  CAPTURE(projection_kernel);
79  CHECK(projection_kernel.cols() == (Ju + Jw - 1));
80  CHECK(projection_kernel.rows() == (Jv + Jw - 1));
81  CHECK(projection_kernel.block(2, 2, Jv, Ju)
82  .isApprox(Matrix<t_complex>::Constant(Jv, Ju, 1.), 1e-12));
83  }
84  SECTION("delta and delta") {
85  const t_int Jw = 5;
86  auto convol_kernelu = [=](const t_int ju) -> t_complex {
87  return (std::floor(ju - (Ju - 1) * 0.5) == 0) ? 1. : 0.;
88  };
89  auto convol_kernelv = [=](const t_int jv) -> t_complex {
90  return (std::floor(jv - (Jv - 1) * 0.5) == 0) ? 1. : 0.;
91  };
92  auto convol_kernelw = [=](const t_int jv, const t_int ju) -> t_complex {
93  return ((std::floor(ju - (Jw - 1) * 0.5) == 0) and (std::floor(jv - (Jw - 1) * 0.5) == 0))
94  ? 1.
95  : 0.;
96  };
97  const Matrix<t_complex> projection_kernel = convol::linear_convol_2d<t_complex>(
98  convol_kernelu, convol_kernelv, convol_kernelw, static_cast<t_int>(Ju),
99  static_cast<t_int>(Jv), static_cast<t_int>(Jw), static_cast<t_int>(Jw));
100  CAPTURE(projection_kernel);
101  CHECK(projection_kernel.cols() == (Ju + Jw - 1));
102  CHECK(projection_kernel.rows() == (Jv + Jw - 1));
103  CHECK(projection_kernel(4, 3) == 1.);
104  }
105  SECTION("delta_size_1") {
106  const t_int Jw = 1;
107  auto convol_kernelu = [=](const t_int ju) -> t_complex {
108  return (std::floor(ju - (Ju - 1) * 0.5) == 0) ? 1. : 0.;
109  };
110  auto convol_kernelv = [=](const t_int jv) -> t_complex {
111  return (std::floor(jv - (Jv - 1) * 0.5) == 0) ? 1. : 0.;
112  };
113  auto convol_kernelw = [=](const t_int jv, const t_int ju) -> t_complex {
114  return ((std::floor(ju - (Jw - 1) * 0.5) == 0) and (std::floor(jv - (Jw - 1) * 0.5) == 0))
115  ? 1.
116  : 0.;
117  };
118  const Matrix<t_complex> projection_kernel = convol::linear_convol_2d<t_complex>(
119  convol_kernelu, convol_kernelv, convol_kernelw, static_cast<t_int>(Ju),
120  static_cast<t_int>(Jv), static_cast<t_int>(Jw), static_cast<t_int>(Jw));
121  CAPTURE(projection_kernel);
122  CHECK(projection_kernel.cols() == (Ju + Jw - 1));
123  CHECK(projection_kernel.rows() == (Jv + Jw - 1));
124  CHECK(projection_kernel(2, 1) == 1.);
125  }
126  SECTION("box_and_delta_size_1") {
127  const t_int Jw = 1;
128  auto convol_kernelu = [=](const t_int ju) -> t_complex { return ju * 1.; };
129  auto convol_kernelv = [=](const t_int jv) -> t_complex { return jv * 1.; };
130  auto convol_kernelw = [=](const t_int jv, const t_int ju) -> t_complex {
131  return ((std::floor(ju - (Jw - 1) * 0.5) == 0) and (std::floor(jv - (Jw - 1) * 0.5) == 0))
132  ? 1.
133  : 0.;
134  };
135  const Matrix<t_complex> projection_kernel = convol::linear_convol_2d<t_complex>(
136  convol_kernelu, convol_kernelv, convol_kernelw, static_cast<t_int>(Ju),
137  static_cast<t_int>(Jv), static_cast<t_int>(Jw), static_cast<t_int>(Jw));
138  CAPTURE(projection_kernel);
139  CHECK(projection_kernel.cols() == (Ju + Jw - 1));
140  CHECK(projection_kernel.rows() == (Jv + Jw - 1));
141  for (int i = 0; i < projection_kernel.cols(); i++) {
142  for (int j = 0; j < projection_kernel.rows(); j++) {
143  CHECK(std::abs(static_cast<t_real>(i * j) - projection_kernel(j, i)) < 1e-12);
144  }
145  }
146  }
147 }

References CHECK.