13 #include <boost/numeric/ublas/matrix.hpp>
20 template <
class Matrix>
21 bool in_bounds(
const Matrix& mat,
size_t row,
size_t col) {
22 return (row < mat.size1()) && (col < mat.size2());
26 bool read_matrix(
const char *filename, boost::numeric::ublas::matrix<double>& mat);
29 template <
class Matrix>
30 void output(
const Matrix& mat, std::ostream& out = std::cout) {
31 const size_t width = 12;
33 for (
size_t i=0; i < mat.size1(); i++) {
35 out << std::setw(width) << mat(i,0);
38 for (
size_t j=1; j < mat.size2(); j++) {
39 out <<
" " << std::setw(width) << mat(i,j);
61 template <
class Matrix>
63 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
64 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max()) {
65 assert(orig.size1() == repro.size1());
66 assert(orig.size2() == repro.size2());
68 if (row_end > orig.size1()) row_end = orig.size1();
69 if (col_end > orig.size2()) col_end = orig.size2();
73 for (
size_t i=row_start; i < row_end; i++) {
74 for (
size_t j=col_start; j < col_end; j++) {
75 double diff = repro(i,j) - orig(i,j);
91 template <
class Matrix>
92 double rmse(
const Matrix& orig,
const Matrix& repro,
93 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
94 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max()) {
95 assert(orig.size1() == repro.size1());
96 assert(orig.size2() == repro.size2());
98 if (row_end > orig.size1()) row_end = orig.size1();
99 if (col_end > orig.size2()) col_end = orig.size2();
102 size_t rows = (row_end - row_start);
103 size_t cols = (col_end - col_start);
109 template <
class Matrix>
110 double nrmse(
const Matrix& orig,
const Matrix& repro,
111 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
112 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max()) {
113 assert(orig.size1() == repro.size1());
114 assert(orig.size2() == repro.size2());
116 if (row_end > orig.size1()) row_end = orig.size1();
117 if (col_end > orig.size2()) col_end = orig.size2();
121 size_t rows = (row_end - row_start);
122 size_t cols = (col_end - col_start);
123 return sqrt(summary.
sum_squares / (rows * cols)) / range;
128 template <
class Matrix>
129 double psnr(
const Matrix& orig,
const Matrix& repro) {
130 assert(orig.size1() == repro.size1());
131 assert(orig.size2() == repro.size2());
134 double range = summary.orig_max - summary.orig_min;
135 return 20 * log10(range / sqrt(summary.sum_squares / (orig.size1() * orig.size2())));
140 template <
class Matrix>
142 assert(orig.size1() == repro.size1());
143 assert(orig.size2() == repro.size2());
146 double range = summary.both_max - summary.both_min;
147 return sqrt(summary.sum_squares / (orig.size1() * orig.size2())) / range;
151 template <
class Matrix>
153 size_t n = mat.size1() * mat.size2();
156 for (
size_t i=0; i < mat.size1(); i++) {
157 for (
size_t j=0; j < mat.size2(); j++) {
159 sum2 += mat(i,j) * mat(i,j);
162 double mean = sum / n;
163 double stdDevInv = 1/sqrt(sum2/n - mean*mean);
165 for (
size_t i=0; i < mat.size1(); i++) {
166 for (
size_t j=0; j < mat.size2(); j++) {
167 mat(i,j) = (mat(i,j) - mean) * stdDevInv;
174 template <
typename T>
194 template <
class Matrix>
195 typename Matrix::value_type
sum(
const Matrix& mat,
196 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
197 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max())
199 if (row_end > mat.size1()) row_end = mat.size1();
200 if (col_end > mat.size2()) col_end = mat.size2();
202 typename Matrix::value_type total = 0;
203 for (
size_t i = row_start; i < row_end; i++) {
204 for (
size_t j = col_start; j < col_end; j++) {
212 template <
class Matrix>
214 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
215 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max())
217 if (row_end > mat.size1()) row_end = mat.size1();
218 if (col_end > mat.size2()) col_end = mat.size2();
220 size_t row_len = (row_end - row_start);
221 size_t col_len = (col_end - col_start);
222 return sum(mat, row_start, row_end, col_start, col_end)
223 / ((double)row_len * col_len);
227 template <
class Matrix>
228 typename Matrix::value_type
max_val(
const Matrix& mat,
229 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
230 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max())
232 if (row_end > mat.size1()) row_end = mat.size1();
233 if (col_end > mat.size2()) col_end = mat.size2();
234 if (row_end <= row_start || col_end <= col_start)
return 0;
236 typename Matrix::value_type max = mat(row_start, col_start);
237 for (
size_t i = row_start; i < row_end; i++) {
238 for (
size_t j = col_start; j < col_end; j++) {
239 max = mat(i,j) > max ? mat(i,j) : max;
246 template <
class Matrix>
247 typename Matrix::value_type
min_val(
const Matrix& mat,
248 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
249 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max())
251 if (row_end > mat.size1()) row_end = mat.size1();
252 if (col_end > mat.size2()) col_end = mat.size2();
253 if (row_end <= row_start || col_end <= col_start)
return 0;
255 typename Matrix::value_type min = mat(row_start, col_start);
256 for (
size_t i = row_start; i < row_end; i++) {
257 for (
size_t j = col_start; j < col_end; j++) {
258 min = mat(i,j) < min ? mat(i,j) : min;
264 template <
class Matrix>
266 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
267 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max())
269 if (row_end > mat.size1()) row_end = mat.size1();
270 if (col_end > mat.size2()) col_end = mat.size2();
272 typename Matrix::value_type amax = 0;
273 for (
size_t i = row_start; i < row_end; i++) {
274 for (
size_t j = col_start; j < col_end; j++) {
275 typename Matrix::value_type tmp =
abs_val(mat(i,j));
276 amax = (tmp > amax) ? tmp : amax;
282 template <
class Matrix,
typename V>
284 size_t row_start = 0,
size_t row_end = std::numeric_limits<size_t>::max(),
285 size_t col_start = 0,
size_t col_end = std::numeric_limits<size_t>::max())
287 if (row_end > mat.size1()) row_end = mat.size1();
288 if (col_end > mat.size2()) col_end = mat.size2();
290 for (
size_t i = row_start; i < row_end; i++) {
291 for (
size_t j = col_start; j < col_end; j++) {
298 template <
typename Iterator1,
typename Iterator2>
301 Iterator1 i = first1;
302 Iterator2 j = first2;
312 template <
typename Iterator1,
typename Iterator2>
315 Iterator1 i = first1;
316 Iterator2 j = first2;
318 double diff = (*i - *j);
330 template <
class Matrix>
334 }
else if (x > (mat.size1() - 1)) {
340 }
else if (y > (mat.size2() - 1)) {
345 size_t x1 = (size_t)floor(x);
346 size_t x2 = (size_t)ceil(x);
348 size_t y1 = (size_t)floor(y);
349 size_t y2 = (size_t)ceil(y);
355 if (xd == 0 && yd == 0) {
358 }
else if (yd == 0) {
359 value = (x2 - x) * mat(x1,y1) + (x - x1) * mat(x2, y2);
361 }
else if (xd == 0) {
362 value = (y2 - y) * mat(x1,y1) + (y - y1) * mat(x2, y2);
367 mat(x1,y1)/d * ((x2-x) * (y2-y)) +
368 mat(x2,y1)/d * ((x-x1) * (y2-y)) +
369 mat(x1,y2)/d * ((x2-x) * (y-y1)) +
370 mat(x2,y2)/d * ((x-x1) * (y-y1));
376 #endif // MATRIX_UTILS_H
void output(const Matrix &mat, std::ostream &out=std::cout)
Matrix::value_type max_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
double euclidean_distance(Iterator1 first1, Iterator1 last1, Iterator2 first2)
double psnr(const Matrix &orig, const Matrix &repro)
Peak Signal to Noise Ratio.
double interp_bilinear(const Matrix &mat, double x, double y)
Interpolates a value for point (x,y) based on values in the Matrix.
Matrix::value_type sum(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Matrix::value_type min_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
double rmse(const Matrix &orig, const Matrix &repro, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
ms_summary get_summary(const Matrix &orig, const Matrix &repro, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
bool read_matrix(const char *filename, boost::numeric::ublas::matrix< double > &mat)
bool in_bounds(const Matrix &mat, size_t row, size_t col)
void set_all(const Matrix &mat, V value, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
void standardize(Matrix &mat)
double mean_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Template class for generic, type-inferred absolute value.
double similarity(const Matrix &orig, const Matrix &repro)
Similarity (NRMSE defined to be symmetric, with max and min taken from both matrices) ...
T abs_val(T num)
Generic, type-inferred absolute value function.
double manhattan_distance(Iterator1 first1, Iterator1 last1, Iterator2 first2)
ms_summary(double ss, double ma, double mi)
Matrix::value_type abs_max_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
double nrmse(const Matrix &orig, const Matrix &repro, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Normalized rms error.
bool isDivisibleBy2(size_t n, int level)
True if and only if n is divisible by 2 <level> times.