179 enum class RANGE { SINGLE, LAST, COORDS, LAST_COORDS, ALL } range_type;
186 struct r_convert_1D final : r_convert {
189 r_convert_1D(
difference_type p1,
const details::last_index&) : r_convert(r_convert::RANGE::LAST_COORDS, p1) { }
190 r_convert_1D(
const details::all_range&) : r_convert(r_convert::RANGE::ALL) { }
192 struct r_convert_2D_1 final : r_convert {
194 r_convert_2D_1(
difference_type p) : r_convert(r_convert::RANGE::SINGLE,p) { }
195 r_convert_2D_1(
const details::last_index&) : r_convert(r_convert::RANGE::LAST) { }
197 r_convert_2D_1(
difference_type p1,
const details::last_index&) : r_convert(r_convert::RANGE::LAST_COORDS, p1) { }
198 r_convert_2D_1(
const details::all_range&) : r_convert(r_convert::RANGE::ALL) { }
200 struct r_convert_2D_2 final : r_convert {
202 r_convert_2D_2(
difference_type p) : r_convert(r_convert::RANGE::SINGLE,p) { }
203 r_convert_2D_2(
const details::last_index&) : r_convert(r_convert::RANGE::LAST) { }
205 r_convert_2D_2(
difference_type p1,
const details::last_index&) : r_convert(r_convert::RANGE::LAST_COORDS, p1) { }
206 r_convert_2D_2(
const details::all_range&) : r_convert(r_convert::RANGE::ALL) { }
228 template <
typename T_container>
230 template <
typename T_container>
235 template <
typename T_container>
237 template <
typename T_container>
250 template<
typename T_output =
interpolator>
251 typename std::enable_if<std::is_floating_point<value_type>::value, T_output>::type
get_interpolator(
INTERP interp_type)
const {
253 switch (interp_type) {
266 template<
typename T_output = linsolver>
267 typename std::enable_if<std::is_floating_point<value_type>::value, T_output>::type
get_linsolver(
LINSOLVER linsolver_type)
const {
269 switch (linsolver_type) {
293 template <
typename T2>
297 friend std::ostream&
operator<<(std::ostream &os,
const Array2D &A) {
return A.this_stream(os); }
301 friend void save(
const Array2D &A,
const std::string &filename) { A.this_save(filename); }
302 friend void save(
const Array2D &A, std::ofstream &os) { A.this_save(os); }
415 bool empty()
const {
return s == 0; }
416 template <
typename T_container>
417 bool same_size(
const T_container &A)
const {
return h == A.h && w == A.w; }
423 std::string
size_2D_string()
const {
return "(" + std::to_string(h) +
"," + std::to_string(w) +
")"; }
437 template<
typename T_output = cv::Mat>
438 typename std::enable_if<std::is_arithmetic<value_type>::value, T_output>::type this_cv_img(
value_type,
value_type)
const;
439 template<
typename T_output =
void>
440 typename std::enable_if<std::is_arithmetic<value_type>::value, T_output>::type this_imshow(
difference_type)
const;
441 template<
typename T_output = std::ostream&>
442 typename std::enable_if<std::is_arithmetic<value_type>::value, T_output>::type this_stream(std::ostream&)
const;
446 template<
typename T_output =
void>
447 typename std::enable_if<std::is_arithmetic<value_type>::value, T_output>::type this_save(
const std::string&)
const;
448 template<
typename T_output =
void>
449 typename std::enable_if<std::is_arithmetic<value_type>::value, T_output>::type this_save(std::ofstream&)
const;
452 bool this_isequal(
const Array2D&)
const;
460 bool this_any_true()
const;
461 bool this_all_true()
const;
501 template<
typename T_output = Array2D&>
502 typename std::enable_if<std::is_arithmetic<value_type>::value, T_output>::type this_sqrt();
503 template<
typename T_output = Array2D&>
504 typename std::enable_if<std::is_arithmetic<value_type>::value, T_output>::type this_pow(
double);
507 template<
typename T_output = Array2D>
508 typename std::enable_if<std::is_same<value_type,double>::value, T_output>::type this_mat_mult(
const Array2D&)
const;
509 template<
typename T_output = Array2D>
510 typename std::enable_if<!std::is_same<value_type,double>::value, T_output>::type this_mat_mult(
const Array2D&)
const;
513 template<
typename T_output = Array2D>
514 typename std::enable_if<std::is_same<value_type,double>::value, T_output>::type this_conv(
const Array2D&)
const;
515 template<
typename T_output = Array2D>
516 typename std::enable_if<std::is_same<value_type,double>::value, T_output>::type this_deconv(
const Array2D&)
const;
517 template<
typename T_output = Array2D>
518 typename std::enable_if<std::is_same<value_type,double>::value, T_output>::type this_xcorr(
const Array2D&)
const;
519 enum class FFTW { CONV, DECONV, XCORR };
523 template<
typename T_output = value_type>
524 typename std::enable_if<std::is_floating_point<value_type>::value, T_output>::type this_dot(
const Array2D&)
const;
525 template<
typename T_output = Array2D>
526 typename std::enable_if<std::is_floating_point<value_type>::value, T_output>::type this_normalize();
527 template<
typename T_output = Array2D>
528 typename std::enable_if<std::is_floating_point<value_type>::value, T_output>::type this_linsolve(
const Array2D&)
const;
531 coords get_range(
const r_convert_1D&)
const;
532 coords get_range(
const r_convert_2D_1&)
const;
533 coords get_range(
const r_convert_2D_2&)
const;
537 template <
typename T_it>
539 void destroy_and_deallocate();
540 void make_null() { ptr =
nullptr; h = w = s = 0; }
543 void chk_samesize_op(
const Array2D&,
const std::string&)
const;
545 void chk_column_op(
const std::string&)
const;
546 void chk_square_op(
const std::string&)
const;
549 void chk_mult_size(
const Array2D&)
const;
550 void chk_kernel_size(
const Array2D&)
const;
563 template<
typename T_container>
564 struct container_traits {
565 typedef typename T_container::pointer nonconst_pointer;
566 typedef typename T_container::reference nonconst_reference;
567 typedef typename T_container::iterator nonconst_iterator;
568 typedef typename T_container::region nonconst_region;
569 typedef typename T_container::container nonconst_container;
571 typedef typename T_container::pointer pointer;
572 typedef typename T_container::reference reference;
573 typedef typename T_container::iterator iterator;
574 typedef typename T_container::region region;
575 typedef typename T_container::container container;
577 typedef typename T_container::const_pointer const_pointer;
578 typedef typename T_container::const_reference const_reference;
579 typedef typename T_container::const_iterator const_iterator;
580 typedef typename T_container::const_region const_region;
581 typedef typename T_container::const_container const_container;
584 template<
typename T_container>
585 struct container_traits<const T_container> {
586 typedef typename T_container::pointer nonconst_pointer;
587 typedef typename T_container::reference nonconst_reference;
588 typedef typename T_container::iterator nonconst_iterator;
589 typedef typename T_container::region nonconst_region;
590 typedef typename T_container::container nonconst_container;
592 typedef typename T_container::const_pointer pointer;
593 typedef typename T_container::const_reference reference;
594 typedef typename T_container::const_iterator iterator;
595 typedef typename T_container::const_region region;
596 typedef typename T_container::const_container container;
598 typedef typename T_container::const_pointer const_pointer;
599 typedef typename T_container::const_reference const_reference;
600 typedef typename T_container::const_iterator const_iterator;
601 typedef typename T_container::const_region const_region;
602 typedef typename T_container::const_container const_container;
606 template <
typename T_container>
614 typedef typename T_container::coords
coords;
615 typedef typename container_traits<T_container>::pointer
pointer;
616 typedef typename container_traits<T_container>::reference
reference;
618 typedef typename container_traits<T_container>::container
container;
621 template <
typename T_container2>
636 template<
typename T_container2>
648 template <
typename T_container2>
652 template <
typename T_container2>
654 return !((*this) == it);
671 template <
typename T_container>
686 template <
typename T_container2>
700 template<
typename T_container2>
713 template <
typename T_container>
728 template <
typename T_container2>
743 template <typename T_container2>
745 base_iterator<
container>(it), r_sub1_2D(it.r_sub1_2D), r_sub2_2D(it.r_sub2_2D), sub_p(it.sub_p), sub_h(it.sub_h), sub_w(it.sub_w), sub_s(it.sub_s) { }
756 void chk_valid_ranges()
const;
766 template <
typename T_container>
782 template <
typename T_container2>
797 template <typename T_container2>
811 void chk_same_size()
const;
816 template <
typename T_iterator>
825 typedef typename T_iterator::coords
coords;
832 template <
typename T_iterator2>
849 template <
typename T_iterator2>
851 ptr(it.ptr ? it.ptr->const_clone() : nullptr) { }
861 template <
typename T_iterator2>
863 return *(it.ptr) == *(ptr);
865 template <
typename T_iterator2>
867 return *(it.ptr) != *(ptr);
871 std::shared_ptr<T_iterator> ptr;
875 template <
typename T_container>
882 typedef typename T_container::coords
coords;
883 typedef typename container_traits<T_container>::pointer
pointer;
884 typedef typename container_traits<T_container>::const_pointer
const_pointer;
885 typedef typename container_traits<T_container>::reference
reference;
887 typedef typename container_traits<T_container>::iterator
iterator;
890 typedef typename container_traits<T_container>::container
container;
893 template <
typename T_container2>
908 template<
typename T_container2>
916 template<
typename T_container2>
950 template <
typename T_container>
968 template <
typename T_container2>
983 template<typename T_container2>
999 void chk_valid_range()
const;
1004 template <
typename T_container>
1022 template <
typename T_container2>
1037 template<typename T_container2>
1039 base_region<
container>(reg), r_sub1_2D(reg.r_sub1_2D), r_sub2_2D(reg.r_sub2_2D), sub_h(reg.sub_h), sub_w(reg.sub_w), sub_s(reg.sub_s) { }
1055 void chk_valid_ranges()
const;
1064 template <
typename T_container>
1083 template <
typename T_container2>
1098 template<typename T_container2>
1114 void chk_same_size()
const;
1116 std::shared_ptr<bool_container> A_bool_ptr;
1119 template <
typename T_region>
1138 template <
typename T_region2>
1158 template <
typename T_region2>
1160 ptr(reg.ptr ? reg.ptr->const_clone() : nullptr) { }
1163 template <
typename T_region2>
1164 typename std::enable_if<std::is_same<typename T_region2::nonconst_container, nonconst_container>::value,
interface_region&>::type
operator=(
const T_region2 ®) { *(ptr) = *(reg.ptr);
return *
this; }
1184 std::shared_ptr<T_region> ptr;
1188 template <
typename T_container>
1198 typedef typename container_traits<T_container>::nonconst_container
container;
1234 template <
typename T_container>
1269 bool out_of_bounds(
double p1,
double p2)
const override {
return p1 <= -0.5 || p2 <= -0.5 || p1 >= this->
A_ptr->height() - 0.5 || p2 >= this->
A_ptr->width() - 0.5; }
1272 template <
typename T_container>
1310 bool out_of_bounds(
double p1,
double p2)
const override {
return p1 < 0 || p1 >= this->
A_ptr->height() - 1 || p2 < 0 || p2 >= this->
A_ptr->width() - 1; }
1313 template <
typename T_container>
1314 class coef_mat_interp_base :
public base_interp<T_container> {
1329 coef_mat_interp_base() : order() { }
1330 coef_mat_interp_base(
const coef_mat_interp_base&) =
default;
1331 coef_mat_interp_base(coef_mat_interp_base&&) noexcept = default;
1332 coef_mat_interp_base& operator=(const coef_mat_interp_base&) = default;
1333 coef_mat_interp_base& operator=(coef_mat_interp_base&&) = default;
1334 virtual ~coef_mat_interp_base() noexcept = default;
1337 coef_mat_interp_base(const_container &A, difference_type order) :
1338 base_interp<container>(A), order(order), p1_pow_buf(order+1,1), p2_pow_buf(order+1,1), p1_pow_dp1_buf(order+1,1), p2_pow_dp2_buf(order+1,1) {
1341 throw std::invalid_argument(
"Attempted to form coef_mat_interp_base with order of: " + std::to_string(order) +
" order must be 1 or greater.");
1346 value_type operator()(
double,
double)
const override;
1347 const_container& first_order(
double,
double)
const override;
1351 virtual container& get_p_pow(container&,
double)
const;
1352 virtual container& get_dp_pow(container&, const_container&)
const;
1353 virtual value_type t_vec_mat_vec(const_container&, const_container&, const_container&)
const;
1354 virtual const_container& calc_coef_mat(container&, const_container&, difference_type, difference_type)
const = 0;
1355 virtual const_container& get_coef_mat(difference_type, difference_type)
const = 0;
1358 mutable container p1_pow_buf;
1359 mutable container p2_pow_buf;
1360 mutable container p1_pow_dp1_buf;
1361 mutable container p2_pow_dp2_buf;
1364 template <
typename T_container>
1365 class cubic_interp_base :
public coef_mat_interp_base<T_container> {
1368 typedef typename coef_mat_interp_base<T_container>::value_type value_type;
1369 typedef typename coef_mat_interp_base<T_container>::reference reference;
1370 typedef typename coef_mat_interp_base<T_container>::size_type size_type;
1371 typedef typename coef_mat_interp_base<T_container>::difference_type
difference_type;
1372 typedef typename coef_mat_interp_base<T_container>::coords coords;
1373 typedef typename coef_mat_interp_base<T_container>::container container;
1374 typedef typename coef_mat_interp_base<T_container>::const_container const_container;
1379 cubic_interp_base() =
default;
1380 cubic_interp_base(
const cubic_interp_base&) =
default;
1381 cubic_interp_base(cubic_interp_base&&) noexcept = default;
1382 cubic_interp_base& operator=(const cubic_interp_base&) = default;
1383 cubic_interp_base& operator=(cubic_interp_base&&) = default;
1384 virtual ~cubic_interp_base() noexcept = default;
1387 cubic_interp_base(const_container &A) : coef_mat_interp_base<container>(A,3) { }
1390 value_type operator()(
double,
double)
const override;
1391 const_container& first_order(
double,
double)
const override;
1395 const_container& calc_coef_mat(container&, const_container&, difference_type, difference_type)
const override;
1401 bool out_of_bounds(
double p1,
double p2)
const override {
return p1 < 1 || p1 >= this->
A_ptr->height() - 2 || p2 < 1 || p2 >= this->
A_ptr->width() - 2; }
1404 template <
typename T_container>
1408 typedef typename cubic_interp_base<T_container>::value_type
value_type;
1409 typedef typename cubic_interp_base<T_container>::reference
reference;
1410 typedef typename cubic_interp_base<T_container>::size_type
size_type;
1412 typedef typename cubic_interp_base<T_container>::coords
coords;
1413 typedef typename cubic_interp_base<T_container>::container
container;
1436 return this->calc_coef_mat(coef_mat_buf, *this->
A_ptr, p1 - 1, p2 - 1);
1442 template <
typename T_container>
1448 typedef typename cubic_interp_base<T_container>::value_type
value_type;
1449 typedef typename cubic_interp_base<T_container>::reference
reference;
1450 typedef typename cubic_interp_base<T_container>::size_type
size_type;
1452 typedef typename cubic_interp_base<T_container>::coords
coords;
1453 typedef typename cubic_interp_base<T_container>::container
container;
1477 std::shared_ptr<Array2D<container>> coef_mat_precompute_ptr;
1480 template <
typename T_container>
1481 class quintic_interp_base :
public coef_mat_interp_base<T_container> {
1484 typedef typename coef_mat_interp_base<T_container>::value_type value_type;
1485 typedef typename coef_mat_interp_base<T_container>::reference reference;
1486 typedef typename coef_mat_interp_base<T_container>::size_type size_type;
1487 typedef typename coef_mat_interp_base<T_container>::difference_type difference_type;
1488 typedef typename coef_mat_interp_base<T_container>::coords coords;
1489 typedef typename coef_mat_interp_base<T_container>::container container;
1490 typedef typename coef_mat_interp_base<T_container>::const_container const_container;
1495 quintic_interp_base() =
default;
1496 quintic_interp_base(
const quintic_interp_base&) =
default;
1497 quintic_interp_base(quintic_interp_base&&) noexcept = default;
1498 quintic_interp_base& operator=(const quintic_interp_base&) = default;
1499 quintic_interp_base& operator=(quintic_interp_base&&) = default;
1500 virtual ~quintic_interp_base() noexcept = default;
1503 quintic_interp_base(const_container &A) : coef_mat_interp_base<container>(A,5) { }
1506 value_type operator()(
double,
double)
const override;
1507 const_container& first_order(
double,
double)
const override;
1511 std::shared_ptr<container> get_bspline_mat_ptr(const_container&)
const;
1512 const_container& calc_coef_mat(container&, const_container&, difference_type, difference_type)
const override;
1518 bool out_of_bounds(
double p1,
double p2)
const override {
return p1 < 0 || p1 > this->
A_ptr->height() - 1 || p2 < 0 || p2 > this->
A_ptr->width() - 1; }
1522 difference_type bcoef_border = 20;
1525 template <
typename T_container>
1529 typedef typename quintic_interp_base<T_container>::value_type
value_type;
1530 typedef typename quintic_interp_base<T_container>::reference
reference;
1531 typedef typename quintic_interp_base<T_container>::size_type
size_type;
1533 typedef typename quintic_interp_base<T_container>::coords
coords;
1534 typedef typename quintic_interp_base<T_container>::container
container;
1557 return this->calc_coef_mat(coef_mat_buf, *bcoef_ptr, p1 + this->bcoef_border - 2, p2 + this->bcoef_border - 2);
1561 std::shared_ptr<container> bcoef_ptr;
1564 template <
typename T_container>
1570 typedef typename quintic_interp_base<T_container>::value_type
value_type;
1571 typedef typename quintic_interp_base<T_container>::reference
reference;
1572 typedef typename quintic_interp_base<T_container>::size_type
size_type;
1574 typedef typename quintic_interp_base<T_container>::coords
coords;
1575 typedef typename quintic_interp_base<T_container>::container
container;
1599 std::shared_ptr<Array2D<container>> coef_mat_precompute_ptr;
1602 template <
typename T_
interp>
1615 template <
typename T_
interp2>
1637 std::shared_ptr<T_interp> ptr;
1641 template <
typename T_container>
1651 typedef typename container_traits<T_container>::nonconst_container
container;
1670 virtual operator bool() = 0;
1685 template <
typename T_container>
1712 operator
bool()
override {
return full_rank; }
1721 std::shared_ptr<container> piv_ptr;
1725 template <
typename T_container>
1752 operator
bool()
override {
return full_rank; }
1762 std::shared_ptr<container> piv_ptr;
1763 std::shared_ptr<container> beta_ptr;
1768 template <
typename T_container>
1795 operator
bool()
override {
return pos_def; }
1807 template <
typename T_linsolver>
1820 template <
typename T_linsolver2>
1841 operator bool() {
return (*ptr); }
1847 std::shared_ptr<T_linsolver> ptr;
1856template <
typename T,
typename T_alloc>
1859 auto new_ptr = allocate_and_copy(A.s,A.ptr);
1860 destroy_and_deallocate();
1865 this->ptr = new_ptr;
1870template <
typename T,
typename T_alloc>
1874 destroy_and_deallocate();
1888template <
typename T,
typename T_alloc>
1890 chk_minsize_op(0,0,
"Array2D construction");
1892 this->ptr = allocate_and_init(this->h*this->w, val);
1895template <
typename T,
typename T_alloc>
1898 std::copy(il.begin(), il.end(), this->ptr);
1901template <
typename T,
typename T_alloc>
1906 for (
auto it_row = il_row.begin(); it_row != il_row.end(); ++it_row) {
1908 auto il_col = *it_row;
1910 if (il_col.size() != this->w) {
1916 throw std::invalid_argument(
"Initializer list for Array2D must not be jagged.");
1920 for (
auto it_col = il_col.begin(); it_col != il_col.end(); ++it_col) {
1921 (*this)(p1,p2) = *it_col;
1929template <
typename T,
typename T_alloc>
1930template <
typename T_output>
1933 std::ifstream is(filename.c_str(), std::ios::in | std::ios::binary);
1934 if (!is.is_open()) {
1935 throw std::invalid_argument(
"Could not open " + filename +
" for loading Array2D.");
1947template <
typename T,
typename T_alloc>
1948template <
typename T_output>
1954 is.read(
reinterpret_cast<char*
>(&A.h), std::streamsize(
sizeof(
difference_type)));
1955 is.read(
reinterpret_cast<char*
>(&A.w), std::streamsize(
sizeof(
difference_type)));
1959 A.ptr = A.allocate_and_init(A.s,
value_type());
1962 is.read(
reinterpret_cast<char*
>(A.ptr), std::streamsize(
sizeof(
value_type) * A.s));
1968template <
typename T,
typename T_alloc>
1970 chk_size_op(1,1,
"value_type conversion");
1976template <
typename T,
typename T_alloc>
1980 chk_minsize_op(1,1,
"1D single-element indexing");
1981 chk_in_bounds_op(p,
"1D single-element indexing");
1987template <
typename T,
typename T_alloc>
1991 chk_minsize_op(1,1,
"2D single-element indexing");
1992 chk_in_bounds_op(p1,p2,
"2D single-element indexing");
1995 return ptr[sub2ind(p1,p2)];
1999template <
typename T,
typename T_alloc>
2000template <
typename T_output>
2006 cv::Mat cv_img(h, w, cv::DataType<uchar>::type);
2007 double conversion = (std::abs(
max -
min) > std::numeric_limits<double>::epsilon()) ? std::numeric_limits<uchar>::max() / (
max -
min) : 0;
2008 for (difference_type p2 = 0; p2 < w; ++p2) {
2009 for (difference_type p1 = 0; p1 < h; ++p1) {
2011 double val = ((*this)(p1,p2) -
min) * conversion;
2013 cv_img.at<uchar>(p1,p2) = 0;
2014 }
else if (val > 255) {
2015 cv_img.at<uchar>(p1,p2) = 255;
2017 cv_img.at<uchar>(p1,p2) = val;
2026template <
typename T,
typename T_alloc>
2027template <
typename T_output>
2028typename std::enable_if<std::is_arithmetic<T>::value, T_output>::type Array2D<T,T_alloc>::this_imshow(difference_type delay)
const {
2030 chk_minsize_op(1,1,
"imshow");
2033 cv::Mat cv_img = this_cv_img(this_min(), this_max());
2034 cv::imshow(
"Array2D", cv_img);
2035 delay == -1 ? cv::waitKey() : cv::waitKey(delay);
2038template <
typename T,
typename T_alloc>
2039template <
typename T_output>
2040typename std::enable_if<std::is_arithmetic<T>::value, T_output>::type Array2D<T,T_alloc>::this_stream(std::ostream &os)
const {
2042 for (difference_type p1 = 0; p1 < h; ++p1) {
2043 for (difference_type p2 = 0; p2 < w; ++p2) {
2044 os << (*this)(p1,p2) << ((p2 < w - 1) ?
'\t' :
'\0');
2052template <
typename T,
typename T_alloc>
2053Array2D<T, T_alloc> Array2D<T, T_alloc>::this_repmat(difference_type rows, difference_type cols)
const {
2054 Array2D B(h * rows, w * cols);
2055 for (difference_type p2 = 0; p2 < cols; ++p2) {
2056 for (difference_type p1 = 0; p1 < rows; ++p1) {
2057 B({p1*h, (p1+1)*h-1},{p2*w, (p2+1)*w-1}) = (*
this);
2064template <
typename T,
typename T_alloc>
2065Array2D<T, T_alloc> Array2D<T, T_alloc>::this_pad(difference_type padding,
PAD type)
const {
2066 Array2D A(h + 2*padding, w + 2*padding);
2067 A({padding, h + padding - 1},{padding, w + padding - 1}) = *
this;
2074 A({0,padding-1},{0,padding-1}) = (*
this)(0,0);
2075 A({0,padding-1},{w+padding,w+2*padding-1}) = (*
this)(0,
last);
2076 A({h+padding,
last},{0,padding-1}) = (*
this)(
last,0);
2077 A({h+padding,
last},{w+padding,w+2*padding-1}) = (*
this)(
last,
last);
2079 A({0,padding-1},{padding,w+padding-1}) = repmat((*
this)(0,
all),padding,1);
2080 A({padding,h+padding-1},{0,padding-1}) = repmat((*
this)(
all,0),1,padding);
2081 A({padding,h+padding-1},{w+padding,
last}) = repmat((*
this)(
all,
last),1,padding);
2082 A({h+padding,
last},{padding,w+padding-1}) = repmat((*
this)(
last,
all),padding,1);
2092template <
typename T,
typename T_alloc>
2093Array2D<T, T_alloc> Array2D<T, T_alloc>::this_t()
const {
2095 for (difference_type p2 = 0; p2 < w; ++p2) {
2096 for (difference_type p1 = 0; p1 < h; ++p1) {
2097 A(p2,p1) = (*this)(p1,p2);
2104template <
typename T,
typename T_alloc>
2105template <
typename T_output>
2106typename std::enable_if<std::is_arithmetic<T>::value, T_output>::type Array2D<T, T_alloc>::this_save(
const std::string &filename)
const {
2108 std::ofstream os(filename.c_str(), std::ios::out | std::ios::binary | std::ios::trunc);
2109 if (!os.is_open()) {
2110 throw std::invalid_argument(
"Could not open " + filename +
" for saving Array2D.");
2120template <
typename T,
typename T_alloc>
2121template <
typename T_output>
2122typename std::enable_if<std::is_arithmetic<T>::value, T_output>::type Array2D<T, T_alloc>::this_save(std::ofstream &os)
const {
2124 os.write(
reinterpret_cast<const char*
>(&h), std::streamsize(
sizeof(difference_type)));
2125 os.write(
reinterpret_cast<const char*
>(&w), std::streamsize(
sizeof(difference_type)));
2126 os.write(
reinterpret_cast<const char*
>(ptr), std::streamsize(
sizeof(value_type) * s));
2131template <
typename T,
typename T_alloc>
2132bool Array2D<T, T_alloc>::this_isequal(
const Array2D &A)
const {
2133 if (h != A.h || w != A.w) {
2137 for (difference_type p = 0; p < s; ++p) {
2138 if ((*
this)(p) != A(p)) {
2146template <
typename T,
typename T_alloc>
2148 chk_samesize_op(A,
"element-wise equality");
2150 bool_container B(h,w);
2151 for (difference_type p = 0; p < s; ++p) {
2152 B(p) = (*this)(p) == A(p);
2158template <
typename T,
typename T_alloc>
2160 bool_container B(h,w);
2161 for (difference_type p = 0; p < s; ++p) {
2162 B(p) = (*this)(p) == val;
2168template <
typename T,
typename T_alloc>
2170 chk_samesize_op(A,
"element-wise inequality");
2172 bool_container B(h,w);
2173 for (difference_type p = 0; p < s; ++p) {
2174 B(p) = (*this)(p) != A(p);
2180template <
typename T,
typename T_alloc>
2182 bool_container B(h,w);
2183 for (difference_type p = 0; p < s; ++p) {
2184 B(p) = (*this)(p) != val;
2190template <
typename T,
typename T_alloc>
2192 chk_samesize_op(A,
"element-wise 'and'");
2194 bool_container B(h,w);
2195 for (difference_type p = 0; p < s; ++p) {
2196 B(p) = (*this)(p) && A(p);
2202template <
typename T,
typename T_alloc>
2204 chk_samesize_op(A,
"element-wise 'or'");
2206 bool_container B(h,w);
2207 for (difference_type p = 0; p < s; ++p) {
2208 B(p) = (*this)(p) || A(p);
2214template <
typename T,
typename T_alloc>
2216 bool_container B(h,w);
2217 for (difference_type p = 0; p < s; ++p) {
2224template <
typename T,
typename T_alloc>
2225bool Array2D<T,T_alloc>::this_any_true()
const {
2226 for (difference_type p = 0; p < s; ++p) {
2235template <
typename T,
typename T_alloc>
2236bool Array2D<T,T_alloc>::this_all_true()
const {
2237 for (difference_type p = 0; p < s; ++p) {
2246template <
typename T,
typename T_alloc>
2249 for (difference_type p = start; p < s; ++p) {
2259template <
typename T,
typename T_alloc>
2261 chk_samesize_op(A,
"element-wise greater-than");
2263 bool_container B(h,w);
2264 for (difference_type p = 0; p < s; ++p) {
2265 B(p) = (*this)(p) > A(p);
2271template <
typename T,
typename T_alloc>
2273 bool_container B(h,w);
2274 for (difference_type p = 0; p < s; ++p) {
2275 B(p) = (*this)(p) > val;
2281template <
typename T,
typename T_alloc>
2283 chk_samesize_op(A,
"element-wise greater-than or equal to");
2285 bool_container B(h,w);
2286 for (difference_type p = 0; p < s; ++p) {
2287 B(p) = (*this)(p) >= A(p);
2293template <
typename T,
typename T_alloc>
2295 bool_container B(h,w);
2296 for (difference_type p = 0; p < s; ++p) {
2297 B(p) = (*this)(p) >= val;
2303template <
typename T,
typename T_alloc>
2305 chk_samesize_op(A,
"element-wise less-than");
2307 bool_container B(h,w);
2308 for (difference_type p = 0; p < s; ++p) {
2309 B(p) = (*this)(p) < A(p);
2315template <
typename T,
typename T_alloc>
2317 bool_container B(h,w);
2318 for (difference_type p = 0; p < s; ++p) {
2319 B(p) = (*this)(p) < val;
2325template <
typename T,
typename T_alloc>
2327 chk_samesize_op(A,
"element-wise less-than or equal to");
2329 bool_container B(h,w);
2330 for (difference_type p = 0; p < s; ++p) {
2331 B(p) = (*this)(p) <= A(p);
2337template <
typename T,
typename T_alloc>
2339 bool_container B(h,w);
2340 for (difference_type p = 0; p < s; ++p) {
2341 B(p) = (*this)(p) <= val;
2348template <
typename T,
typename T_alloc>
2349Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator+=(
const Array2D &A) {
2350 chk_samesize_op(A,
"element-wise addition");
2352 for (difference_type p = 0; p < s; ++p) {
2359template <
typename T,
typename T_alloc>
2360Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator+=(const_reference val) {
2361 for (difference_type p = 0; p < s; ++p) {
2368template <
typename T,
typename T_alloc>
2369Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator-=(
const Array2D &A) {
2370 chk_samesize_op(A,
"element-wise subtraction");
2372 for (difference_type p = 0; p < s; ++p) {
2379template <
typename T,
typename T_alloc>
2380Array2D<T,T_alloc>& Array2D<T,T_alloc>::this_flipsubassign(
const Array2D &A) {
2381 chk_samesize_op(A,
"element-wise subtraction");
2383 for (difference_type p = 0; p < s; ++p) {
2384 (*this)(p) = A(p)-(*this)(p);
2390template <
typename T,
typename T_alloc>
2391Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator-=(const_reference val) {
2392 for (difference_type p = 0; p < s; ++p) {
2399template <
typename T,
typename T_alloc>
2400Array2D<T,T_alloc>& Array2D<T,T_alloc>::this_flipsubassign(const_reference val) {
2401 for (difference_type p = 0; p < s; ++p) {
2402 (*this)(p) = val-(*
this)(p);
2408template <
typename T,
typename T_alloc>
2409Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator*=(
const Array2D &A) {
2410 chk_samesize_op(A,
"element-wise multiplication");
2412 for (difference_type p = 0; p < s; ++p) {
2419template <
typename T,
typename T_alloc>
2420Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator*=(const_reference val) {
2421 for (difference_type p = 0; p < s; ++p) {
2428template <
typename T,
typename T_alloc>
2429Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator/=(
const Array2D &A) {
2430 chk_samesize_op(A,
"element-wise division");
2432 for (difference_type p = 0; p < s; ++p) {
2439template <
typename T,
typename T_alloc>
2440Array2D<T,T_alloc>& Array2D<T,T_alloc>::this_flipdivassign(
const Array2D &A) {
2441 chk_samesize_op(A,
"element-wise division");
2443 for (difference_type p = 0; p < s; ++p) {
2444 (*this)(p) = A(p)/(*this)(p);
2450template <
typename T,
typename T_alloc>
2451Array2D<T,T_alloc>& Array2D<T,T_alloc>::operator/=(const_reference val) {
2452 for (difference_type p = 0; p < s; ++p) {
2459template <
typename T,
typename T_alloc>
2460Array2D<T,T_alloc>& Array2D<T,T_alloc>::this_flipdivassign(const_reference val) {
2461 for (difference_type p = 0; p < s; ++p) {
2462 (*this)(p) = val/(*
this)(p);
2468template <
typename T,
typename T_alloc>
2473template <
typename T,
typename T_alloc>
2475 for (difference_type p = 0; p < s; ++p) {
2476 (*this)(p) = -(*
this)(p);
2482template <
typename T,
typename T_alloc>
2483Array2D<T, T_alloc>& Array2D<T,T_alloc>::this_sort() {
2484 std::sort(ptr, ptr + s);
2489template <
typename T,
typename T_alloc>
2490T Array2D<T,T_alloc>::this_sum()
const {
2491 value_type sum = value_type();
2492 for (difference_type p = 0; p < s; ++p) {
2499template <
typename T,
typename T_alloc>
2500T Array2D<T,T_alloc>::this_max()
const {
2501 chk_minsize_op(1,1,
"max");
2503 return *std::max_element(ptr, ptr + s);
2506template <
typename T,
typename T_alloc>
2507T Array2D<T,T_alloc>::this_min()
const {
2508 chk_minsize_op(1,1,
"min");
2510 return *std::min_element(ptr, ptr + s);
2513template <
typename T,
typename T_alloc>
2514T Array2D<T,T_alloc>::this_prctile(
double percent) {
2515 chk_minsize_op(1,1,
"prctile");
2516 if (percent < 0 || percent > 1) {
2517 throw std::invalid_argument(
"Input percent of: " + std::to_string(percent) +
" provided to prctile operator. Value must be between (inclusive) 0 and 1.");
2522 std::sort(ptr, ptr + s);
2524 return (*
this)((s-1) * percent);
2527template <
typename T,
typename T_alloc>
2528template<
typename T_output>
2529typename std::enable_if<std::is_arithmetic<T>::value, T_output>::type Array2D<T,T_alloc>::this_sqrt() {
2530 for (difference_type p = 0; p < s; ++p) {
2531 (*this)(p) = std::sqrt((*
this)(p));
2537template <
typename T,
typename T_alloc>
2538template<
typename T_output>
2539typename std::enable_if<std::is_arithmetic<T>::value, T_output>::type Array2D<T,T_alloc>::this_pow(
double n) {
2540 for (difference_type p = 0; p < s; ++p) {
2541 (*this)(p) = std::pow((*
this)(p),n);
2547template <
typename T,
typename T_alloc>
2548template <
typename T_output>
2549typename std::enable_if<!std::is_same<T,double>::value, T_output>::type Array2D<T,T_alloc>::this_mat_mult(
const Array2D<T,T_alloc> &A)
const {
2555 for (difference_type p = 0; p < A.w; ++p) {
2556 for (difference_type p2 = 0; p2 < w; ++p2) {
2557 for (difference_type p1 = 0; p1 < h; ++p1) {
2558 B(p1,p) += (*this)(p1,p2) * A(p2,p);
2566template <
typename T,
typename T_alloc>
2567template <
typename T_output>
2568typename std::enable_if<std::is_same<T,double>::value, T_output>::type Array2D<T,T_alloc>::this_mat_mult(
const Array2D<T,T_alloc> &A)
const {
2584 dgemm_(&TRANS, &TRANS, &M, &N, &K, &ALPHA, ptr, &M, A.ptr, &K, &BETA, B.ptr, &M);
2593 template <
typename T>
2594 class fftw_allocator final {
2596 typedef T value_type;
2598 typedef const T* const_pointer;
2599 typedef T& reference;
2600 typedef const T& const_reference;
2601 typedef std::size_t size_type;
2605 fftw_allocator() noexcept = default;
2606 fftw_allocator(const fftw_allocator&) noexcept = default;
2607 fftw_allocator(fftw_allocator&&) noexcept = default;
2608 fftw_allocator& operator=(const fftw_allocator&) = default;
2609 fftw_allocator& operator=(fftw_allocator&&) = default;
2610 ~fftw_allocator() noexcept = default;
2615 pointer allocate(difference_type s, const
void * = 0) {
2616 pointer ptr =
static_cast<pointer
>(malloc_function(s *
sizeof(value_type)));
2618 NLOG_ERROR <<
"Failed to allocate memory using allocate in fftw_allocator.";
2619 throw std::bad_alloc();
2624 void deallocate(pointer ptr, difference_type) { free_function(ptr); }
2625 void construct(pointer ptr, const_reference val = value_type()) { ::new((
void *)ptr) value_type(val); }
2626 void destroy(pointer ptr) { ptr->~value_type(); }
2627 template<
typename T2>
2633 void* malloc_function(difference_type s) {
return malloc(s); }
2634 void free_function(pointer ptr) {
return free(ptr); }
2640 inline void* fftw_allocator<double>::malloc_function(difference_type s) {
return fftw_malloc(s); }
2642 inline void fftw_allocator<double>::free_function(pointer ptr) { fftw_free(ptr); }
2648template <
typename T,
typename T_alloc>
2649template <
typename T_output>
2650typename std::enable_if<std::is_same<T,double>::value, T_output>::type Array2D<T,T_alloc>::this_conv(
const Array2D &A)
const {
return this_conv_base(A, FFTW::CONV); }
2652template <
typename T,
typename T_alloc>
2653template <
typename T_output>
2654typename std::enable_if<std::is_same<T,double>::value, T_output>::type Array2D<T,T_alloc>::this_deconv(
const Array2D &A)
const {
return this_conv_base(A, FFTW::DECONV); }
2656template <
typename T,
typename T_alloc>
2657template <
typename T_output>
2658typename std::enable_if<std::is_same<T,double>::value, T_output>::type Array2D<T,T_alloc>::this_xcorr(
const Array2D &A)
const {
return this_conv_base(A, FFTW::XCORR); }
2663template <
typename T,
typename T_alloc>
2664Array2D<T,T_alloc> Array2D<T,T_alloc>::this_conv_base(
const Array2D &kernel, FFTW fftw_type)
const {
2665 using details::fftw_allocator;
2667 chk_kernel_size(kernel);
2678 double* A_rm =
static_cast<double*
>(fftw_malloc(
sizeof(
double) * H * W));
2679 double* K_rm =
static_cast<double*
>(fftw_malloc(
sizeof(
double) * H * W));
2680 if (!A_rm || !K_rm) {
2681 if (A_rm) fftw_free(A_rm);
2682 if (K_rm) fftw_free(K_rm);
2683 throw std::bad_alloc();
2686 std::fill(A_rm, A_rm + H*W, 0.0);
2687 std::fill(K_rm, K_rm + H*W, 0.0);
2690 for (difference_type r = 0; r < H; ++r) {
2691 for (difference_type c = 0; c < W; ++c) {
2692 A_rm[r*W + c] = (*this)(r, c);
2703 for (difference_type r = 0; r <= c_kh; ++r) {
2704 for (difference_type c = 0; c <= c_kw; ++c) {
2705 K_rm[r*W + c] = kernel(c_kh + r, c_kw + c);
2710 for (difference_type r = 0; r <= c_kh; ++r) {
2711 for (difference_type c = 0; c < c_kw; ++c) {
2712 K_rm[r*W + (W - c_kw + c)] = kernel(c_kh + r, c);
2718 for (difference_type r = 0; r < c_kh; ++r) {
2719 for (difference_type c = 0; c <= c_kw; ++c) {
2720 K_rm[(H - c_kh + r)*W + c] = kernel(r, c_kh + c);
2725 if (c_kh > 0 && c_kw > 0) {
2726 for (difference_type r = 0; r < c_kh; ++r) {
2727 for (difference_type c = 0; c < c_kw; ++c) {
2728 K_rm[(H - c_kh + r)*W + (W - c_kw + c)] = kernel(r, c);
2734 fftw_complex* A_fft =
reinterpret_cast<fftw_complex*
>(fftw_malloc(
sizeof(fftw_complex) * H * halfWplus1));
2735 fftw_complex* K_fft =
reinterpret_cast<fftw_complex*
>(fftw_malloc(
sizeof(fftw_complex) * H * halfWplus1));
2736 if (!A_fft || !K_fft) {
2737 if (A_fft) fftw_free(A_fft);
2738 if (K_fft) fftw_free(K_fft);
2741 throw std::bad_alloc();
2745 fftw_plan plan_A = fftw_plan_dft_r2c_2d(
static_cast<int>(H),
static_cast<int>(W), A_rm, A_fft, FFTW_ESTIMATE);
2746 fftw_plan plan_kernel = fftw_plan_dft_r2c_2d(
static_cast<int>(H),
static_cast<int>(W), K_rm, K_fft, FFTW_ESTIMATE);
2747 fftw_plan plan_output = fftw_plan_dft_c2r_2d(
static_cast<int>(H),
static_cast<int>(W), A_fft, A_rm, FFTW_ESTIMATE);
2750 auto plan_deleter = [](fftw_plan *p) {
return fftw_destroy_plan(*p); };
2751 std::unique_ptr<fftw_plan,
decltype(plan_deleter)> plan_A_delete(&plan_A, plan_deleter);
2752 std::unique_ptr<fftw_plan,
decltype(plan_deleter)> plan_kernel_delete(&plan_kernel, plan_deleter);
2753 std::unique_ptr<fftw_plan,
decltype(plan_deleter)> plan_output_delete(&plan_output, plan_deleter);
2759 fftw_execute(plan_A);
2760 fftw_execute(plan_kernel);
2763 for (difference_type r = 0; r < H; ++r) {
2765 for (difference_type c = 0; c < halfWplus1; ++c) {
2767 double ar = A_fft[idx][0];
2768 double ai = A_fft[idx][1];
2769 double kr = K_fft[idx][0];
2770 double ki = K_fft[idx][1];
2771 switch (fftw_type) {
2773 double nr = ar*kr - ai*ki;
2774 double ni = ar*ki + ai*kr;
2779 case FFTW::DECONV: {
2780 double den = kr*kr + ki*ki;
2782 A_fft[idx][0] = 0.0;
2783 A_fft[idx][1] = 0.0;
2785 double nr = (ar*kr + ai*ki)/den;
2786 double ni = (ai*kr - ar*ki)/den;
2794 double nr = ar*kr + ai*ki;
2795 double ni = -ar*ki + ai*kr;
2805 fftw_execute(plan_output);
2812 const double scale = 1.0 /
static_cast<double>(H * W);
2813 for (difference_type r = 0; r < H; ++r) {
2814 for (difference_type c = 0; c < W; ++c) {
2815 C(r, c) = A_rm[r*W + c] * scale;
2829template <
typename T,
typename T_alloc>
2830template <
typename T_output>
2831typename std::enable_if<std::is_floating_point<T>::value, T_output>::type Array2D<T,T_alloc>::this_dot(
const Array2D &x)
const {
2832 chk_samesize_op(x,
"dot product");
2834 value_type sum = value_type();
2835 for (difference_type p = 0; p < s; ++p) {
2836 sum += (*this)(p)*x(p);
2842template <
typename T,
typename T_alloc>
2843template <
typename T_output>
2844typename std::enable_if<std::is_floating_point<T>::value, T_output>::type Array2D<T,T_alloc>::this_normalize() {
2845 value_type norm = value_type();
2846 for (difference_type p = 0; p < s; ++p) {
2847 norm += std::pow((*
this)(p),2.0);
2849 norm = std::sqrt(norm);
2851 return (*
this) /= norm;
2854template <
typename T,
typename T_alloc>
2855template <
typename T_output>
2856typename std::enable_if<std::is_floating_point<T>::value, T_output>::type Array2D<T,T_alloc>::this_linsolve(
const Array2D &b)
const {
2862 container x(w, b.width());
2863 for (difference_type p2 = 0; p2 < b.width(); ++p2) {
2864 x(
all,p2) = linsolver.solve(b(
all,p2));
2871template <
typename T,
typename T_alloc>
2876 switch (r.range_type) {
2877 case r_convert::RANGE::SINGLE :
2879 throw std::invalid_argument(
"Attempted to use 'single' element range for 1D range");
2880 case r_convert::RANGE::LAST :
2882 throw std::invalid_argument(
"Attempted to use 'last' element range for 1D range");
2883 case r_convert::RANGE::COORDS :
2887 case r_convert::RANGE::LAST_COORDS :
2891 case r_convert::RANGE::ALL :
2903 if (p1 < 0 || p2 >= s) {
2904 throw std::invalid_argument(
"Attempted to perform 1D range type indexing using: ({" + std::to_string(p1) +
"," + std::to_string(p2) +
2905 "}) which is beyond the size of the array: " + size_string());
2913template <
typename T,
typename T_alloc>
2918 switch (r.range_type) {
2919 case r_convert::RANGE::SINGLE :
2923 case r_convert::RANGE::LAST :
2927 case r_convert::RANGE::COORDS :
2931 case r_convert::RANGE::LAST_COORDS :
2935 case r_convert::RANGE::ALL :
2947 if (p1 < 0 || p2 >= h) {
2948 throw std::invalid_argument(
"Attempted to perform 2D range type indexing using: ({" + std::to_string(p1) +
"," + std::to_string(p2) +
2949 "},) which is beyond the size of the Array: " + size_2D_string() +
".");
2957template <
typename T,
typename T_alloc>
2962 switch (r.range_type) {
2963 case r_convert::RANGE::SINGLE :
2967 case r_convert::RANGE::LAST :
2971 case r_convert::RANGE::COORDS :
2975 case r_convert::RANGE::LAST_COORDS :
2979 case r_convert::RANGE::ALL :
2991 if (p1 < 0 || p2 >= w) {
2992 throw std::invalid_argument(
"Attempted to perform 2D range type indexing using: (,{" + std::to_string(p1) +
"," + std::to_string(p2) +
2993 "}) which is beyond the size of the Array: " + size_2D_string() +
".");
3001template <
typename T,
typename T_alloc>
3005 ptr_new = alloc.allocate(s_init);
3006 }
catch (std::bad_alloc &e) {
3008 NLOG_ERROR <<
"Error: failed to allocate memory for Array2D.";
3015template <
typename T,
typename T_alloc>
3017 pointer ptr_new = allocate(s_init);
3018 std::uninitialized_fill_n(ptr_new, s_init, val);
3023template <
typename T,
typename T_alloc>
3024template <
typename T_it>
3026 pointer ptr_new = allocate(s_init);
3027 std::uninitialized_copy_n(it, s_init, ptr_new);
3032template <
typename T,
typename T_alloc>
3033void Array2D<T,T_alloc>::destroy_and_deallocate() {
3039 for (pointer ptr_delete = ptr + s; ptr_delete != ptr; ) {
3040 alloc.destroy(--ptr_delete);
3042 alloc.deallocate(ptr, s);
3046template <
typename T,
typename T_alloc>
3047void Array2D<T,T_alloc>::chk_size_op(difference_type h, difference_type w,
const std::string &op)
const {
3048 if (this->h != h || this->w != w) {
3049 throw std::invalid_argument(
"Attempted to use " + op +
" operator on array of size " + size_2D_string() +
3050 ". Array must have size of (" + std::to_string(h) +
"," + std::to_string(w) +
").");
3054template <
typename T,
typename T_alloc>
3055void Array2D<T,T_alloc>::chk_samesize_op(
const Array2D &A,
const std::string &op)
const {
3056 if (!same_size(A)) {
3057 throw std::invalid_argument(
"Attempted to use " + op +
" operator on array of size " + size_2D_string() +
3058 " with array of size " + A.size_2D_string() +
". Arrays must have the same size.");
3062template <
typename T,
typename T_alloc>
3063void Array2D<T,T_alloc>::chk_minsize_op(difference_type h, difference_type w,
const std::string &op)
const {
3064 if (this->h < h || this->w < w) {
3065 throw std::invalid_argument(
"Attempted to use " + op +
" operator on array of size " + size_2D_string() +
3066 ". Array must be of size (" + std::to_string(h) +
"," + std::to_string(w) +
") or greater.");
3070template <
typename T,
typename T_alloc>
3071void Array2D<T,T_alloc>::chk_column_op(
const std::string &op)
const {
3073 throw std::invalid_argument(
"Attempted to use " + op +
" operator on array of size " + size_2D_string() +
3074 ". Array must be a column.");
3078template <
typename T,
typename T_alloc>
3079void Array2D<T,T_alloc>::chk_square_op(
const std::string &op)
const {
3080 if (this->h != this->w) {
3081 throw std::invalid_argument(
"Attempted to use " + op +
" operator on array of size " + size_2D_string() +
3082 ". Array must be square.");
3086template <
typename T,
typename T_alloc>
3087void Array2D<T,T_alloc>::chk_in_bounds_op(difference_type p,
const std::string &op)
const {
3088 if (!in_bounds(p)) {
3089 throw std::invalid_argument(
"Attempted to use " + op +
" with input of " + std::to_string(p) +
" on array of size " +
3090 size_string() +
".");
3094template <
typename T,
typename T_alloc>
3095void Array2D<T,T_alloc>::chk_in_bounds_op(difference_type p1, difference_type p2,
const std::string &op)
const {
3096 if (!in_bounds(p1,p2)) {
3097 throw std::invalid_argument(
"Attempted to use " + op +
" with input of (" + std::to_string(p1) +
"," + std::to_string(p2) +
") on array of size " +
3098 size_2D_string() +
".");
3102template <
typename T,
typename T_alloc>
3103void Array2D<T,T_alloc>::chk_mult_size(
const Array2D &A)
const {
3104 if (this->w != A.h) {
3105 throw std::invalid_argument(
"Attempted to multiply matrix of size " + size_2D_string() +
" with " + A.size_2D_string() +
3106 ". These sizes are incompatible for matrix multiplication.");
3110template <
typename T,
typename T_alloc>
3111void Array2D<T,T_alloc>::chk_kernel_size(
const Array2D &kernel)
const {
3113 if (!(kernel.h % 2) || !(kernel.w % 2) || kernel.w > this->w || kernel.h > this->h) {
3114 throw std::invalid_argument(
"Attempted to convolve matrix of size " + size_2D_string() +
" with kernel of size " + kernel.size_2D_string() +
3115 ". Kernel must have odd dimensions and be equal to or smaller than matrix.");
3121 template <
typename T_container>
3128 template <
typename T_container>
3130 if (p >= A_ptr->size()) {
3131 throw std::invalid_argument(
"Attempted to increment Array2D iterator beyond last element.");
3135 template <
typename T_container>
3138 throw std::invalid_argument(
"Attempted to decrement Array2D iterator beyond first element.");
3142 template <
typename T_container>
3144 if (!A_ptr->in_bounds(p)) {
3145 throw std::invalid_argument(
"Attempted to dereference Array2D iterator out of range.");
3150 template <
typename T_container>
3152 this->chk_valid_increment();
3159 template <
typename T_container>
3161 this->chk_valid_decrement();
3169 template <
typename T_container>
3172 r_sub1_2D(r_sub1_2D),
3173 r_sub2_2D(r_sub2_2D),
3174 sub_p(p_2D.first-r_sub1_2D.first + (p_2D.second-r_sub2_2D.first)*(r_sub1_2D.second-r_sub1_2D.first)),
3175 sub_h(r_sub1_2D.second-r_sub1_2D.first),
3176 sub_w(r_sub2_2D.second-r_sub2_2D.first),
3177 sub_s((r_sub1_2D.second-r_sub1_2D.first) * (r_sub2_2D.second-r_sub2_2D.first)) {
3182 template <
typename T_container>
3184 this->chk_valid_increment();
3188 this->p = this->A_ptr->sub2ind(sub_p % sub_h + r_sub1_2D.first, sub_p / sub_h + r_sub2_2D.first);
3193 template <
typename T_container>
3195 this->chk_valid_decrement();
3199 this->p = this->A_ptr->sub2ind(sub_p % sub_h + r_sub1_2D.first, sub_p / sub_h + r_sub2_2D.first);
3204 template <
typename T_container>
3206 if (r_sub1_2D.first > r_sub1_2D.second || r_sub1_2D.first < 0 || r_sub1_2D.second > this->A_ptr->height() ||
3207 r_sub2_2D.first > r_sub2_2D.second || r_sub2_2D.first < 0 || r_sub2_2D.second > this->A_ptr->width()) {
3208 throw std::invalid_argument(
"Range of ({" + std::to_string(r_sub1_2D.first) +
"," + std::to_string(r_sub1_2D.second-1) +
"},{" +
3209 std::to_string(r_sub2_2D.first) +
"," + std::to_string(r_sub2_2D.second-1) +
"}) is not valid for sub iterator with size of: " +
3210 this->A_ptr->size_2D_string() +
".");
3215 template <
typename T_container>
3223 if (this->p == 0 && !this->
A_ptr->empty() && !(*this->A_bool_ptr)(0)) {
3228 template <
typename T_container>
3230 this->chk_valid_increment();
3233 while (++this->p < this->A_ptr->size()) {
3234 if ((*A_bool_ptr)(this->p)) {
3241 template <
typename T_container>
3243 this->chk_valid_decrement();
3246 while (this->p > 0) {
3247 if ((*A_bool_ptr)(--this->p)) {
3255 template <
typename T_container>
3257 if (!this->A_ptr->same_size(*A_bool_ptr)) {
3258 throw std::invalid_argument(
"Attempted to use boolean Array of size: " + A_bool_ptr->size_2D_string() +
3259 " to form iterator for Array2D of size: " + this->A_ptr->size_2D_string() +
".");
3264 template <
typename T_container>
3268 " to region of size: " + region_size_2D_string() +
".");
3271 std::copy(reg.
begin(), reg.
end(), this->begin());
3276 template <
typename T_container>
3277 template <
typename T_container2>
3281 " to region of size: " + region_size_2D_string() +
".");
3284 std::copy(reg.
begin(), reg.
end(), this->begin());
3289 template <
typename T_container>
3291 if (this->region_h != A.height() || this->region_w != A.width()) {
3292 throw std::invalid_argument(
"Attempted to assign Array2D of size: " + A.size_2D_string() +
3293 " to region of size: " + region_size_2D_string() +
".");
3296 std::copy(A.begin(), A.end(), this->begin());
3301 template <
typename T_container>
3303 std::fill(this->begin(), this->end(), val);
3309 template <
typename T_container>
3315 template <
typename T_container>
3317 if (r.first > r.second || r.first < 0 || r.second > this->A_ptr->size()) {
3318 throw std::invalid_argument(
"Range of: ({" + std::to_string(r.first) +
"," + std::to_string(r.second-1) +
3319 "}) is not valid for region with max size of: " + this->A_ptr->size_string() +
".");
3324 template <
typename T_container>
3327 r_sub1_2D(r_sub1_2D),
3328 r_sub2_2D(r_sub2_2D),
3329 sub_h(r_sub1_2D.second-r_sub1_2D.first),
3330 sub_w(r_sub2_2D.second-r_sub2_2D.first),
3331 sub_s((r_sub1_2D.second-r_sub1_2D.first) * (r_sub2_2D.second-r_sub2_2D.first)) {
3336 template <
typename T_container>
3338 if (r_sub1_2D.first > r_sub1_2D.second || r_sub1_2D.first < 0 || r_sub1_2D.second > this->A_ptr->height() ||
3339 r_sub2_2D.first > r_sub2_2D.second || r_sub2_2D.first < 0 || r_sub2_2D.second > this->A_ptr->width()) {
3340 throw std::invalid_argument(
"Range of ({" + std::to_string(r_sub1_2D.first) +
"," + std::to_string(r_sub1_2D.second-1) +
"},{" +
3341 std::to_string(r_sub2_2D.first) +
"," + std::to_string(r_sub2_2D.second-1) +
"}) is not valid for sub region with max size of: " +
3342 this->A_ptr->size_2D_string() +
".");
3347 template <
typename T_container>
3353 template <
typename T_container>
3355 if (!this->A_ptr->same_size(*A_bool_ptr)) {
3356 throw std::invalid_argument(
"Attempted to use boolean Array of size: " + A_bool_ptr->size_2D_string() +
3357 " to index Array2D of size: " + this->A_ptr->size_2D_string() +
".");
3362 template <
typename T_container>
3364 if (out_of_bounds(p1, p2)) {
3365 return std::numeric_limits<value_type>::quiet_NaN();
3368 return (*this->A_ptr)(std::round(p1),std::round(p2));
3371 template <
typename T_container>
3373 if (out_of_bounds(p1, p2)) {
3374 this->first_order_buf(0) = this->first_order_buf(1) = this->first_order_buf(2) = std::numeric_limits<value_type>::quiet_NaN();
3375 return this->first_order_buf;
3378 this->first_order_buf(0) = (*this->A_ptr)(std::round(p1),std::round(p2));
3379 this->first_order_buf(1) = this->first_order_buf(2) = 0;
3381 return this->first_order_buf;
3385 template <
typename T_container>
3387 if (out_of_bounds(p1,p2)) {
3388 return std::numeric_limits<value_type>::quiet_NaN();
3391 double delta_p1 = p1 - std::floor(p1);
3392 double delta_p2 = p2 - std::floor(p2);
3394 return (*this->A_ptr)(p1,p2) * (1-delta_p1) * (1-delta_p2) +
3395 (*this->A_ptr)(p1+1,p2) * delta_p1 * (1-delta_p2) +
3396 (*this->A_ptr)(p1,p2+1) * (1-delta_p1) * delta_p2 +
3397 (*this->A_ptr)(p1+1,p2+1) * delta_p1 * delta_p2;
3400 template <
typename T_container>
3402 if (out_of_bounds(p1,p2)) {
3403 this->first_order_buf(0) = this->first_order_buf(1) = this->first_order_buf(2) = std::numeric_limits<value_type>::quiet_NaN();
3404 return this->first_order_buf;
3407 double delta_p1 = p1 - std::floor(p1);
3408 double delta_p2 = p2 - std::floor(p2);
3410 this->first_order_buf(0) = (*this->A_ptr)(p1,p2) * (1-delta_p1) * (1-delta_p2) +
3411 (*this->A_ptr)(p1+1,p2) * delta_p1 * (1-delta_p2) +
3412 (*this->A_ptr)(p1,p2+1) * (1-delta_p1) * delta_p2 +
3413 (*this->A_ptr)(p1+1,p2+1) * delta_p1 * delta_p2;
3415 this->first_order_buf(1) = (*this->A_ptr)(p1,p2) * (-1) * (1-delta_p2) +
3416 (*this->A_ptr)(p1+1,p2) * 1 * (1-delta_p2) +
3417 (*this->A_ptr)(p1,p2+1) * (-1) * delta_p2 +
3418 (*this->A_ptr)(p1+1,p2+1) * 1 * delta_p2;
3420 this->first_order_buf(2) = (*this->A_ptr)(p1,p2) * (1-delta_p1) * (-1) +
3421 (*this->A_ptr)(p1+1,p2) * delta_p1 * (-1) +
3422 (*this->A_ptr)(p1,p2+1) * (1-delta_p1) * 1 +
3423 (*this->A_ptr)(p1+1,p2+1) * delta_p1 * 1;
3425 return this->first_order_buf;
3429 template <
typename T_container>
3430 inline typename coef_mat_interp_base<T_container>::value_type coef_mat_interp_base<T_container>::operator()(
double p1,
double p2)
const {
3434 if (this->out_of_bounds(p1,p2)) {
3435 return std::numeric_limits<value_type>::quiet_NaN();
3439 p1_pow_buf = get_p_pow(p1_pow_buf, p1 - std::floor(p1));
3442 p2_pow_buf = get_p_pow(p2_pow_buf, p2 - std::floor(p2));
3445 const auto &coef_mat = this->get_coef_mat(p1, p2);
3448 return t_vec_mat_vec(p1_pow_buf, coef_mat, p2_pow_buf);
3451 template <
typename T_container>
3452 inline typename coef_mat_interp_base<T_container>::const_container& coef_mat_interp_base<T_container>::first_order(
double p1,
double p2)
const {
3456 if (this->out_of_bounds(p1,p2)) {
3457 this->first_order_buf(0) = this->first_order_buf(1) = this->first_order_buf(2) = std::numeric_limits<value_type>::quiet_NaN();
3458 return this->first_order_buf;
3462 p1_pow_buf = get_p_pow(p1_pow_buf, p1 - std::floor(p1));
3465 p1_pow_dp1_buf = get_dp_pow(p1_pow_dp1_buf, p1_pow_buf);
3468 p2_pow_buf = get_p_pow(p2_pow_buf, p2 - std::floor(p2));
3471 p2_pow_dp2_buf = get_dp_pow(p2_pow_dp2_buf, p2_pow_buf);
3474 const auto &coef_mat = this->get_coef_mat(p1, p2);
3477 this->first_order_buf(0) = t_vec_mat_vec(p1_pow_buf, coef_mat, p2_pow_buf);
3478 this->first_order_buf(1) = t_vec_mat_vec(p1_pow_dp1_buf, coef_mat, p2_pow_buf);
3479 this->first_order_buf(2) = t_vec_mat_vec(p1_pow_buf, coef_mat, p2_pow_dp2_buf);
3481 return this->first_order_buf;
3484 template <
typename T_container>
3485 inline typename coef_mat_interp_base<T_container>::container& coef_mat_interp_base<T_container>::get_p_pow(container &p_pow,
double delta_p)
const {
3490 for (difference_type i = 2; i <= this->order; ++i) {
3491 p_pow(i) = p_pow(i-1) * delta_p;
3497 template <
typename T_container>
3498 inline typename coef_mat_interp_base<T_container>::container& coef_mat_interp_base<T_container>::get_dp_pow(container &p_pow_dp, const_container &p_pow)
const {
3503 for (difference_type i = 2; i <= this->order; ++i) {
3504 p_pow_dp(i) = p_pow(i-1) * i;
3510 template <
typename T_container>
3511 inline typename coef_mat_interp_base<T_container>::value_type coef_mat_interp_base<T_container>::t_vec_mat_vec(const_container &vec1, const_container &mat, const_container &vec2)
const {
3512 return value_type( t( vec1 ) * mat * vec2 );
3516 template <
typename T_container>
3517 inline typename cubic_interp_base<T_container>::value_type cubic_interp_base<T_container>::operator()(
double p1,
double p2)
const {
3518 if (out_of_bounds(p1,p2)) {
3519 return std::numeric_limits<value_type>::quiet_NaN();
3522 double delta_p1 = p1 - std::floor(p1);
3523 double delta_p2 = p2 - std::floor(p2);
3525 this->p1_pow_buf(0) = 1.0;
3526 this->p1_pow_buf(1) = delta_p1;
3527 this->p1_pow_buf(2) = this->p1_pow_buf(1)*delta_p1;
3528 this->p1_pow_buf(3) = this->p1_pow_buf(2)*delta_p1;
3530 this->p2_pow_buf(0) = 1.0;
3531 this->p2_pow_buf(1) = delta_p2;
3532 this->p2_pow_buf(2) = this->p2_pow_buf(1)*delta_p2;
3533 this->p2_pow_buf(3) = this->p2_pow_buf(2)*delta_p2;
3535 const auto &coef_mat = this->get_coef_mat(p1, p2);
3537 return (this->p2_pow_buf(0)*coef_mat(0,0)+this->p2_pow_buf(1)*coef_mat(0,1)+this->p2_pow_buf(2)*coef_mat(0,2)+this->p2_pow_buf(3)*coef_mat(0,3))*this->p1_pow_buf(0)+
3538 (this->p2_pow_buf(0)*coef_mat(1,0)+this->p2_pow_buf(1)*coef_mat(1,1)+this->p2_pow_buf(2)*coef_mat(1,2)+this->p2_pow_buf(3)*coef_mat(1,3))*this->p1_pow_buf(1)+
3539 (this->p2_pow_buf(0)*coef_mat(2,0)+this->p2_pow_buf(1)*coef_mat(2,1)+this->p2_pow_buf(2)*coef_mat(2,2)+this->p2_pow_buf(3)*coef_mat(2,3))*this->p1_pow_buf(2)+
3540 (this->p2_pow_buf(0)*coef_mat(3,0)+this->p2_pow_buf(1)*coef_mat(3,1)+this->p2_pow_buf(2)*coef_mat(3,2)+this->p2_pow_buf(3)*coef_mat(3,3))*this->p1_pow_buf(3);
3543 template <
typename T_container>
3544 inline typename cubic_interp_base<T_container>::const_container& cubic_interp_base<T_container>::first_order(
double p1,
double p2)
const {
3545 if (out_of_bounds(p1,p2)) {
3546 this->first_order_buf(0) = this->first_order_buf(1) = this->first_order_buf(2) = std::numeric_limits<value_type>::quiet_NaN();
3547 return this->first_order_buf;
3550 double delta_p1 = p1 - std::floor(p1);
3551 double delta_p2 = p2 - std::floor(p2);
3553 this->p1_pow_buf(0) = 1.0;
3554 this->p1_pow_buf(1) = delta_p1;
3555 this->p1_pow_buf(2) = this->p1_pow_buf(1) * delta_p1;
3556 this->p1_pow_buf(3) = this->p1_pow_buf(2) * delta_p1;
3558 this->p1_pow_dp1_buf(0) = 0.0;
3559 this->p1_pow_dp1_buf(1) = 1.0;
3560 this->p1_pow_dp1_buf(2) = 2.0 * this->p1_pow_buf(1);
3561 this->p1_pow_dp1_buf(3) = 3.0 * this->p1_pow_buf(2);
3563 this->p2_pow_buf(0) = 1.0;
3564 this->p2_pow_buf(1) = delta_p2;
3565 this->p2_pow_buf(2) = this->p2_pow_buf(1) * delta_p2;
3566 this->p2_pow_buf(3) = this->p2_pow_buf(2) * delta_p2;
3568 this->p2_pow_dp2_buf(0) = 0.0;
3569 this->p2_pow_dp2_buf(1) = 1.0;
3570 this->p2_pow_dp2_buf(2) = 2.0 * this->p2_pow_buf(1);
3571 this->p2_pow_dp2_buf(3) = 3.0 * this->p2_pow_buf(2);
3573 const auto &coef_mat = this->get_coef_mat(p1, p2);
3575 this->first_order_buf(0) = (this->p2_pow_buf(0)*coef_mat(0,0)+this->p2_pow_buf(1)*coef_mat(0,1)+this->p2_pow_buf(2)*coef_mat(0,2)+this->p2_pow_buf(3)*coef_mat(0,3))*this->p1_pow_buf(0)+
3576 (this->p2_pow_buf(0)*coef_mat(1,0)+this->p2_pow_buf(1)*coef_mat(1,1)+this->p2_pow_buf(2)*coef_mat(1,2)+this->p2_pow_buf(3)*coef_mat(1,3))*this->p1_pow_buf(1)+
3577 (this->p2_pow_buf(0)*coef_mat(2,0)+this->p2_pow_buf(1)*coef_mat(2,1)+this->p2_pow_buf(2)*coef_mat(2,2)+this->p2_pow_buf(3)*coef_mat(2,3))*this->p1_pow_buf(2)+
3578 (this->p2_pow_buf(0)*coef_mat(3,0)+this->p2_pow_buf(1)*coef_mat(3,1)+this->p2_pow_buf(2)*coef_mat(3,2)+this->p2_pow_buf(3)*coef_mat(3,3))*this->p1_pow_buf(3);
3580 this->first_order_buf(1) = (this->p2_pow_buf(0)*coef_mat(0,0)+this->p2_pow_buf(1)*coef_mat(0,1)+this->p2_pow_buf(2)*coef_mat(0,2)+this->p2_pow_buf(3)*coef_mat(0,3))*this->p1_pow_dp1_buf(0)+
3581 (this->p2_pow_buf(0)*coef_mat(1,0)+this->p2_pow_buf(1)*coef_mat(1,1)+this->p2_pow_buf(2)*coef_mat(1,2)+this->p2_pow_buf(3)*coef_mat(1,3))*this->p1_pow_dp1_buf(1)+
3582 (this->p2_pow_buf(0)*coef_mat(2,0)+this->p2_pow_buf(1)*coef_mat(2,1)+this->p2_pow_buf(2)*coef_mat(2,2)+this->p2_pow_buf(3)*coef_mat(2,3))*this->p1_pow_dp1_buf(2)+
3583 (this->p2_pow_buf(0)*coef_mat(3,0)+this->p2_pow_buf(1)*coef_mat(3,1)+this->p2_pow_buf(2)*coef_mat(3,2)+this->p2_pow_buf(3)*coef_mat(3,3))*this->p1_pow_dp1_buf(3);
3585 this->first_order_buf(2) = (this->p2_pow_dp2_buf(0)*coef_mat(0,0)+this->p2_pow_dp2_buf(1)*coef_mat(0,1)+this->p2_pow_dp2_buf(2)*coef_mat(0,2)+this->p2_pow_dp2_buf(3)*coef_mat(0,3))*this->p1_pow_buf(0)+
3586 (this->p2_pow_dp2_buf(0)*coef_mat(1,0)+this->p2_pow_dp2_buf(1)*coef_mat(1,1)+this->p2_pow_dp2_buf(2)*coef_mat(1,2)+this->p2_pow_dp2_buf(3)*coef_mat(1,3))*this->p1_pow_buf(1)+
3587 (this->p2_pow_dp2_buf(0)*coef_mat(2,0)+this->p2_pow_dp2_buf(1)*coef_mat(2,1)+this->p2_pow_dp2_buf(2)*coef_mat(2,2)+this->p2_pow_dp2_buf(3)*coef_mat(2,3))*this->p1_pow_buf(2)+
3588 (this->p2_pow_dp2_buf(0)*coef_mat(3,0)+this->p2_pow_dp2_buf(1)*coef_mat(3,1)+this->p2_pow_dp2_buf(2)*coef_mat(3,2)+this->p2_pow_dp2_buf(3)*coef_mat(3,3))*this->p1_pow_buf(3);
3590 return this->first_order_buf;
3593 template <
typename T_container>
3594 inline typename cubic_interp_base<T_container>::const_container& cubic_interp_base<T_container>::calc_coef_mat(container &coef_mat_buf, const_container &A, difference_type p1, difference_type p2)
const {
3597 if (p1 < 0 || p1 + 3 >= A.height() || p2 < 0 || p2 + 3 >= A.width()) {
3598 throw std::invalid_argument(
"p1 and p2 are outside range of array for calc_coef_mat() with cubic interpolation - this is a programmer error.");
3602 coef_mat_buf(0,0) = 1.0*A(p1+1,p2+1);
3603 coef_mat_buf(1,0) = 0.5*A(p1+2,p2+1)-0.5*A(p1,p2+1);
3604 coef_mat_buf(2,0) = 1.0*A(p1,p2+1)-2.5*A(p1+1,p2+1)+2.0*A(p1+2,p2+1)-0.5*A(p1+3,p2+1);
3605 coef_mat_buf(3,0) = 1.5*A(p1+1,p2+1)-0.5*A(p1,p2+1)-1.5*A(p1+2,p2+1)+0.5*A(p1+3,p2+1);
3606 coef_mat_buf(0,1) = 0.5*A(p1+1,p2+2)-0.5*A(p1+1,p2);
3607 coef_mat_buf(1,1) = 0.25*A(p1,p2)-0.25*A(p1,p2+2)-0.25*A(p1+2,p2)+0.25*A(p1+2,p2+2);
3608 coef_mat_buf(2,1) = 0.5*A(p1,p2+2)-0.5*A(p1,p2)+1.25*A(p1+1,p2)-1.25*A(p1+1,p2+2)-1.0*A(p1+2,p2)+1.0*A(p1+2,p2+2)+0.25*A(p1+3,p2)-0.25*A(p1+3,p2+2);
3609 coef_mat_buf(3,1) = 0.25*A(p1,p2)-0.25*A(p1,p2+2)-0.75*A(p1+1,p2)+0.75*A(p1+1,p2+2)+0.75*A(p1+2,p2)-0.75*A(p1+2,p2+2)-0.25*A(p1+3,p2)+0.25*A(p1+3,p2+2);
3610 coef_mat_buf(0,2) = 1.0*A(p1+1,p2)-2.5*A(p1+1,p2+1)+2.0*A(p1+1,p2+2)-0.5*A(p1+1,p2+3);
3611 coef_mat_buf(1,2) = 1.25*A(p1,p2+1)-0.5*A(p1,p2)-1.0*A(p1,p2+2)+0.25*A(p1,p2+3)+0.5*A(p1+2,p2)-1.25*A(p1+2,p2+1)+1.0*A(p1+2,p2+2)-0.25*A(p1+2,p2+3);
3612 coef_mat_buf(2,2) = 1.0*A(p1,p2)-2.5*A(p1,p2+1)+2.0*A(p1,p2+2)-0.5*A(p1,p2+3)-2.5*A(p1+1,p2)+6.25*A(p1+1,p2+1)-5.0*A(p1+1,p2+2)+1.25*A(p1+1,p2+3)+2.0*A(p1+2,p2)-5.0*A(p1+2,p2+1)+4.0*A(p1+2,p2+2)-1.0*A(p1+2,p2+3)-0.5*A(p1+3,p2)+1.25*A(p1+3,p2+1)-1.0*A(p1+3,p2+2)+0.25*A(p1+3,p2+3);
3613 coef_mat_buf(3,2) = 1.25*A(p1,p2+1)-0.5*A(p1,p2)-1.0*A(p1,p2+2)+0.25*A(p1,p2+3)+1.5*A(p1+1,p2)-3.75*A(p1+1,p2+1)+3.0*A(p1+1,p2+2)-0.75*A(p1+1,p2+3)-1.5*A(p1+2,p2)+3.75*A(p1+2,p2+1)-3.0*A(p1+2,p2+2)+0.75*A(p1+2,p2+3)+0.5*A(p1+3,p2)-1.25*A(p1+3,p2+1)+1.0*A(p1+3,p2+2)-0.25*A(p1+3,p2+3);
3614 coef_mat_buf(0,3) = 1.5*A(p1+1,p2+1)-0.5*A(p1+1,p2)-1.5*A(p1+1,p2+2)+0.5*A(p1+1,p2+3);
3615 coef_mat_buf(1,3) = 0.25*A(p1,p2)-0.75*A(p1,p2+1)+0.75*A(p1,p2+2)-0.25*A(p1,p2+3)-0.25*A(p1+2,p2)+0.75*A(p1+2,p2+1)-0.75*A(p1+2,p2+2)+0.25*A(p1+2,p2+3);
3616 coef_mat_buf(2,3) = 1.5*A(p1,p2+1)-0.5*A(p1,p2)-1.5*A(p1,p2+2)+0.5*A(p1,p2+3)+1.25*A(p1+1,p2)-3.75*A(p1+1,p2+1)+3.75*A(p1+1,p2+2)-1.25*A(p1+1,p2+3)-1.0*A(p1+2,p2)+3.0*A(p1+2,p2+1)-3.0*A(p1+2,p2+2)+1.0*A(p1+2,p2+3)+0.25*A(p1+3,p2)-0.75*A(p1+3,p2+1)+0.75*A(p1+3,p2+2)-0.25*A(p1+3,p2+3);
3617 coef_mat_buf(3,3) = 0.25*A(p1,p2)-0.75*A(p1,p2+1)+0.75*A(p1,p2+2)-0.25*A(p1,p2+3)-0.75*A(p1+1,p2)+2.25*A(p1+1,p2+1)-2.25*A(p1+1,p2+2)+0.75*A(p1+1,p2+3)+0.75*A(p1+2,p2)-2.25*A(p1+2,p2+1)+2.25*A(p1+2,p2+2)-0.75*A(p1+2,p2+3)-0.25*A(p1+3,p2)+0.75*A(p1+3,p2+1)-0.75*A(p1+3,p2+2)+0.25*A(p1+3,p2+3);
3619 return coef_mat_buf;
3623 template <
typename T_container>
3627 coef_mat_precompute_ptr = std::make_shared<Array2D<container>>(this->
A_ptr->height(), this->
A_ptr->width(),
container(4,4));
3629 for (difference_type p2 = 1; p2 < this->
A_ptr->width() - 2; ++p2) {
3631 this->calc_coef_mat((*coef_mat_precompute_ptr)(p1,p2), *this->
A_ptr, p1 - 1, p2 - 1);
3637 template <
typename T_container>
3638 inline typename quintic_interp_base<T_container>::value_type quintic_interp_base<T_container>::operator()(
double p1,
double p2)
const {
3639 if (out_of_bounds(p1,p2)) {
3640 return std::numeric_limits<value_type>::quiet_NaN();
3643 double delta_p1 = p1 - std::floor(p1);
3644 double delta_p2 = p2 - std::floor(p2);
3646 this->p1_pow_buf(0) = 1.0;
3647 this->p1_pow_buf(1) = delta_p1;
3648 this->p1_pow_buf(2) = this->p1_pow_buf(1)*delta_p1;
3649 this->p1_pow_buf(3) = this->p1_pow_buf(2)*delta_p1;
3650 this->p1_pow_buf(4) = this->p1_pow_buf(3)*delta_p1;
3651 this->p1_pow_buf(5) = this->p1_pow_buf(4)*delta_p1;
3653 this->p2_pow_buf(0) = 1.0;
3654 this->p2_pow_buf(1) = delta_p2;
3655 this->p2_pow_buf(2) = this->p2_pow_buf(1)*delta_p2;
3656 this->p2_pow_buf(3) = this->p2_pow_buf(2)*delta_p2;
3657 this->p2_pow_buf(4) = this->p2_pow_buf(3)*delta_p2;
3658 this->p2_pow_buf(5) = this->p2_pow_buf(4)*delta_p2;
3660 const auto &coef_mat = this->get_coef_mat(p1, p2);
3662 return (this->p2_pow_buf(0)*coef_mat(0,0)+this->p2_pow_buf(1)*coef_mat(0,1)+this->p2_pow_buf(2)*coef_mat(0,2)+this->p2_pow_buf(3)*coef_mat(0,3)+this->p2_pow_buf(4)*coef_mat(0,4)+this->p2_pow_buf(5)*coef_mat(0,5))*this->p1_pow_buf(0)+
3663 (this->p2_pow_buf(0)*coef_mat(1,0)+this->p2_pow_buf(1)*coef_mat(1,1)+this->p2_pow_buf(2)*coef_mat(1,2)+this->p2_pow_buf(3)*coef_mat(1,3)+this->p2_pow_buf(4)*coef_mat(1,4)+this->p2_pow_buf(5)*coef_mat(1,5))*this->p1_pow_buf(1)+
3664 (this->p2_pow_buf(0)*coef_mat(2,0)+this->p2_pow_buf(1)*coef_mat(2,1)+this->p2_pow_buf(2)*coef_mat(2,2)+this->p2_pow_buf(3)*coef_mat(2,3)+this->p2_pow_buf(4)*coef_mat(2,4)+this->p2_pow_buf(5)*coef_mat(2,5))*this->p1_pow_buf(2)+
3665 (this->p2_pow_buf(0)*coef_mat(3,0)+this->p2_pow_buf(1)*coef_mat(3,1)+this->p2_pow_buf(2)*coef_mat(3,2)+this->p2_pow_buf(3)*coef_mat(3,3)+this->p2_pow_buf(4)*coef_mat(3,4)+this->p2_pow_buf(5)*coef_mat(3,5))*this->p1_pow_buf(3)+
3666 (this->p2_pow_buf(0)*coef_mat(4,0)+this->p2_pow_buf(1)*coef_mat(4,1)+this->p2_pow_buf(2)*coef_mat(4,2)+this->p2_pow_buf(3)*coef_mat(4,3)+this->p2_pow_buf(4)*coef_mat(4,4)+this->p2_pow_buf(5)*coef_mat(4,5))*this->p1_pow_buf(4)+
3667 (this->p2_pow_buf(0)*coef_mat(5,0)+this->p2_pow_buf(1)*coef_mat(5,1)+this->p2_pow_buf(2)*coef_mat(5,2)+this->p2_pow_buf(3)*coef_mat(5,3)+this->p2_pow_buf(4)*coef_mat(5,4)+this->p2_pow_buf(5)*coef_mat(5,5))*this->p1_pow_buf(5);
3670 template <
typename T_container>
3671 inline typename quintic_interp_base<T_container>::const_container& quintic_interp_base<T_container>::first_order(
double p1,
double p2)
const {
3672 if (out_of_bounds(p1,p2)) {
3673 this->first_order_buf(0) = this->first_order_buf(1) = this->first_order_buf(2) = std::numeric_limits<value_type>::quiet_NaN();
3674 return this->first_order_buf;
3677 double delta_p1 = p1 - std::floor(p1);
3678 double delta_p2 = p2 - std::floor(p2);
3680 this->p1_pow_buf(0) = 1.0;
3681 this->p1_pow_buf(1) = delta_p1;
3682 this->p1_pow_buf(2) = this->p1_pow_buf(1)*delta_p1;
3683 this->p1_pow_buf(3) = this->p1_pow_buf(2)*delta_p1;
3684 this->p1_pow_buf(4) = this->p1_pow_buf(3)*delta_p1;
3685 this->p1_pow_buf(5) = this->p1_pow_buf(4)*delta_p1;
3687 this->p1_pow_dp1_buf(0) = 0.0;
3688 this->p1_pow_dp1_buf(1) = 1.0;
3689 this->p1_pow_dp1_buf(2) = 2.0 * this->p1_pow_buf(1);
3690 this->p1_pow_dp1_buf(3) = 3.0 * this->p1_pow_buf(2);
3691 this->p1_pow_dp1_buf(4) = 4.0 * this->p1_pow_buf(3);
3692 this->p1_pow_dp1_buf(5) = 5.0 * this->p1_pow_buf(4);
3694 this->p2_pow_buf(0) = 1.0;
3695 this->p2_pow_buf(1) = delta_p2;
3696 this->p2_pow_buf(2) = this->p2_pow_buf(1)*delta_p2;
3697 this->p2_pow_buf(3) = this->p2_pow_buf(2)*delta_p2;
3698 this->p2_pow_buf(4) = this->p2_pow_buf(3)*delta_p2;
3699 this->p2_pow_buf(5) = this->p2_pow_buf(4)*delta_p2;
3701 this->p2_pow_dp2_buf(0) = 0.0;
3702 this->p2_pow_dp2_buf(1) = 1.0;
3703 this->p2_pow_dp2_buf(2) = 2.0 * this->p2_pow_buf(1);
3704 this->p2_pow_dp2_buf(3) = 3.0 * this->p2_pow_buf(2);
3705 this->p2_pow_dp2_buf(4) = 4.0 * this->p2_pow_buf(3);
3706 this->p2_pow_dp2_buf(5) = 5.0 * this->p2_pow_buf(4);
3708 const auto &coef_mat = this->get_coef_mat(p1, p2);
3710 this->first_order_buf(0) = (this->p2_pow_buf(0)*coef_mat(0,0)+this->p2_pow_buf(1)*coef_mat(0,1)+this->p2_pow_buf(2)*coef_mat(0,2)+this->p2_pow_buf(3)*coef_mat(0,3)+this->p2_pow_buf(4)*coef_mat(0,4)+this->p2_pow_buf(5)*coef_mat(0,5))*this->p1_pow_buf(0)+
3711 (this->p2_pow_buf(0)*coef_mat(1,0)+this->p2_pow_buf(1)*coef_mat(1,1)+this->p2_pow_buf(2)*coef_mat(1,2)+this->p2_pow_buf(3)*coef_mat(1,3)+this->p2_pow_buf(4)*coef_mat(1,4)+this->p2_pow_buf(5)*coef_mat(1,5))*this->p1_pow_buf(1)+
3712 (this->p2_pow_buf(0)*coef_mat(2,0)+this->p2_pow_buf(1)*coef_mat(2,1)+this->p2_pow_buf(2)*coef_mat(2,2)+this->p2_pow_buf(3)*coef_mat(2,3)+this->p2_pow_buf(4)*coef_mat(2,4)+this->p2_pow_buf(5)*coef_mat(2,5))*this->p1_pow_buf(2)+
3713 (this->p2_pow_buf(0)*coef_mat(3,0)+this->p2_pow_buf(1)*coef_mat(3,1)+this->p2_pow_buf(2)*coef_mat(3,2)+this->p2_pow_buf(3)*coef_mat(3,3)+this->p2_pow_buf(4)*coef_mat(3,4)+this->p2_pow_buf(5)*coef_mat(3,5))*this->p1_pow_buf(3)+
3714 (this->p2_pow_buf(0)*coef_mat(4,0)+this->p2_pow_buf(1)*coef_mat(4,1)+this->p2_pow_buf(2)*coef_mat(4,2)+this->p2_pow_buf(3)*coef_mat(4,3)+this->p2_pow_buf(4)*coef_mat(4,4)+this->p2_pow_buf(5)*coef_mat(4,5))*this->p1_pow_buf(4)+
3715 (this->p2_pow_buf(0)*coef_mat(5,0)+this->p2_pow_buf(1)*coef_mat(5,1)+this->p2_pow_buf(2)*coef_mat(5,2)+this->p2_pow_buf(3)*coef_mat(5,3)+this->p2_pow_buf(4)*coef_mat(5,4)+this->p2_pow_buf(5)*coef_mat(5,5))*this->p1_pow_buf(5);
3717 this->first_order_buf(1) = (this->p2_pow_buf(0)*coef_mat(0,0)+this->p2_pow_buf(1)*coef_mat(0,1)+this->p2_pow_buf(2)*coef_mat(0,2)+this->p2_pow_buf(3)*coef_mat(0,3)+this->p2_pow_buf(4)*coef_mat(0,4)+this->p2_pow_buf(5)*coef_mat(0,5))*this->p1_pow_dp1_buf(0)+
3718 (this->p2_pow_buf(0)*coef_mat(1,0)+this->p2_pow_buf(1)*coef_mat(1,1)+this->p2_pow_buf(2)*coef_mat(1,2)+this->p2_pow_buf(3)*coef_mat(1,3)+this->p2_pow_buf(4)*coef_mat(1,4)+this->p2_pow_buf(5)*coef_mat(1,5))*this->p1_pow_dp1_buf(1)+
3719 (this->p2_pow_buf(0)*coef_mat(2,0)+this->p2_pow_buf(1)*coef_mat(2,1)+this->p2_pow_buf(2)*coef_mat(2,2)+this->p2_pow_buf(3)*coef_mat(2,3)+this->p2_pow_buf(4)*coef_mat(2,4)+this->p2_pow_buf(5)*coef_mat(2,5))*this->p1_pow_dp1_buf(2)+
3720 (this->p2_pow_buf(0)*coef_mat(3,0)+this->p2_pow_buf(1)*coef_mat(3,1)+this->p2_pow_buf(2)*coef_mat(3,2)+this->p2_pow_buf(3)*coef_mat(3,3)+this->p2_pow_buf(4)*coef_mat(3,4)+this->p2_pow_buf(5)*coef_mat(3,5))*this->p1_pow_dp1_buf(3)+
3721 (this->p2_pow_buf(0)*coef_mat(4,0)+this->p2_pow_buf(1)*coef_mat(4,1)+this->p2_pow_buf(2)*coef_mat(4,2)+this->p2_pow_buf(3)*coef_mat(4,3)+this->p2_pow_buf(4)*coef_mat(4,4)+this->p2_pow_buf(5)*coef_mat(4,5))*this->p1_pow_dp1_buf(4)+
3722 (this->p2_pow_buf(0)*coef_mat(5,0)+this->p2_pow_buf(1)*coef_mat(5,1)+this->p2_pow_buf(2)*coef_mat(5,2)+this->p2_pow_buf(3)*coef_mat(5,3)+this->p2_pow_buf(4)*coef_mat(5,4)+this->p2_pow_buf(5)*coef_mat(5,5))*this->p1_pow_dp1_buf(5);
3724 this->first_order_buf(2) = (this->p2_pow_dp2_buf(0)*coef_mat(0,0)+this->p2_pow_dp2_buf(1)*coef_mat(0,1)+this->p2_pow_dp2_buf(2)*coef_mat(0,2)+this->p2_pow_dp2_buf(3)*coef_mat(0,3)+this->p2_pow_dp2_buf(4)*coef_mat(0,4)+this->p2_pow_dp2_buf(5)*coef_mat(0,5))*this->p1_pow_buf(0)+
3725 (this->p2_pow_dp2_buf(0)*coef_mat(1,0)+this->p2_pow_dp2_buf(1)*coef_mat(1,1)+this->p2_pow_dp2_buf(2)*coef_mat(1,2)+this->p2_pow_dp2_buf(3)*coef_mat(1,3)+this->p2_pow_dp2_buf(4)*coef_mat(1,4)+this->p2_pow_dp2_buf(5)*coef_mat(1,5))*this->p1_pow_buf(1)+
3726 (this->p2_pow_dp2_buf(0)*coef_mat(2,0)+this->p2_pow_dp2_buf(1)*coef_mat(2,1)+this->p2_pow_dp2_buf(2)*coef_mat(2,2)+this->p2_pow_dp2_buf(3)*coef_mat(2,3)+this->p2_pow_dp2_buf(4)*coef_mat(2,4)+this->p2_pow_dp2_buf(5)*coef_mat(2,5))*this->p1_pow_buf(2)+
3727 (this->p2_pow_dp2_buf(0)*coef_mat(3,0)+this->p2_pow_dp2_buf(1)*coef_mat(3,1)+this->p2_pow_dp2_buf(2)*coef_mat(3,2)+this->p2_pow_dp2_buf(3)*coef_mat(3,3)+this->p2_pow_dp2_buf(4)*coef_mat(3,4)+this->p2_pow_dp2_buf(5)*coef_mat(3,5))*this->p1_pow_buf(3)+
3728 (this->p2_pow_dp2_buf(0)*coef_mat(4,0)+this->p2_pow_dp2_buf(1)*coef_mat(4,1)+this->p2_pow_dp2_buf(2)*coef_mat(4,2)+this->p2_pow_dp2_buf(3)*coef_mat(4,3)+this->p2_pow_dp2_buf(4)*coef_mat(4,4)+this->p2_pow_dp2_buf(5)*coef_mat(4,5))*this->p1_pow_buf(4)+
3729 (this->p2_pow_dp2_buf(0)*coef_mat(5,0)+this->p2_pow_dp2_buf(1)*coef_mat(5,1)+this->p2_pow_dp2_buf(2)*coef_mat(5,2)+this->p2_pow_dp2_buf(3)*coef_mat(5,3)+this->p2_pow_dp2_buf(4)*coef_mat(5,4)+this->p2_pow_dp2_buf(5)*coef_mat(5,5))*this->p1_pow_buf(5);
3731 return this->first_order_buf;
3746 template <
typename T_storage>
3748 std::ptrdiff_t stride = 1) {
3753 constexpr double z[2] = {
3754 -0.43057534709997432,
3755 -0.04309628799923444
3757 constexpr int NB_POLES = 2;
3762 double lambda = 1.0;
3763 for (
int p = 0; p < NB_POLES; ++p) {
3764 lambda *= (1.0 - z[p]) * (1.0 - 1.0 / z[p]);
3766 for (std::ptrdiff_t i = 0; i < N; ++i) {
3767 s[i * stride] = T_storage(
double(s[i * stride]) * lambda);
3771 for (
int p = 0; p < NB_POLES; ++p) {
3772 const double zp = z[p];
3775 constexpr double tol = 1e-9;
3776 std::ptrdiff_t horizon =
static_cast<std::ptrdiff_t
>(
3777 std::ceil(std::log(tol) / std::log(std::abs(zp))));
3778 if (horizon > N) horizon = N;
3779 double sum_zk = double(s[0]);
3781 for (std::ptrdiff_t k = 1; k < horizon; ++k) {
3782 sum_zk += double(s[k * stride]) * zk;
3785 s[0] = T_storage(sum_zk);
3788 for (std::ptrdiff_t i = 1; i < N; ++i) {
3789 s[i * stride] = T_storage(
double(s[i * stride]) +
3790 zp *
double(s[(i - 1) * stride]));
3794 s[(N - 1) * stride] = T_storage(
3795 (zp / (zp * zp - 1.0)) *
3796 (double(s[(N - 1) * stride]) + zp * double(s[(N - 2) * stride])));
3799 for (std::ptrdiff_t i = N - 2; i >= 0; --i) {
3800 s[i * stride] = T_storage(
3801 zp * (
double(s[(i + 1) * stride]) -
double(s[i * stride])));
3806 template <
typename T_container>
3807 std::shared_ptr<typename quintic_interp_base<T_container>::container> quintic_interp_base<T_container>::get_bspline_mat_ptr(const_container &A)
const {
3809 if (bcoef_border < 3) {
3810 throw std::invalid_argument(
"B-coefficient border cannot be less than three when calling get_bspline_coef() - this is a programmer error.");
3821 const auto H = padded.height();
3822 const auto W = padded.width();
3823 if (H < 2 || W < 2) {
3824 return std::make_shared<container>(std::move(padded));
3830 for (
typename container::difference_type p1 = 0; p1 < H; ++p1) {
3831 quintic_bspline_recursive_1d<double>(&padded(p1, 0), W, H);
3834 for (
typename container::difference_type p2 = 0; p2 < W; ++p2) {
3835 quintic_bspline_recursive_1d<double>(&padded(0, p2), H, 1);
3838 return std::make_shared<container>(std::move(padded));
3841 template <
typename T_container>
3842 inline typename quintic_interp_base<T_container>::const_container& quintic_interp_base<T_container>::calc_coef_mat(container &coef_mat_buf, const_container &bcoef, difference_type p1, difference_type p2)
const {
3845 if (p1 < 0 || p1+5 >= bcoef.height() || p2 < 0 || p2+5 >= bcoef.width()) {
3846 throw std::invalid_argument(
"p1 and p2 are outside range of array for calc_coef_mat() with quintic interpolation - this is a programmer error.");
3850 coef_mat_buf(0,0) = 0.00006944444444444444*bcoef(p1,p2)+0.001805555555555556*bcoef(p1+1,p2)+0.001805555555555556*bcoef(p1+4,p2+1)+0.004583333333333333*bcoef(p1,p2+2)+0.1191666666666667*bcoef(p1+1,p2+2)+0.3025*bcoef(p1+2,p2+2)+0.1191666666666667*bcoef(p1+3,p2+2)+0.004583333333333333*bcoef(p1+4,p2+2)+0.001805555555555556*bcoef(p1,p2+3)+0.04694444444444444*bcoef(p1+1,p2+3)+0.004583333333333333*bcoef(p1+2,p2)+0.1191666666666667*bcoef(p1+2,p2+3)+0.04694444444444444*bcoef(p1+3,p2+3)+0.001805555555555556*bcoef(p1+4,p2+3)+0.00006944444444444444*bcoef(p1,p2+4)+0.001805555555555556*bcoef(p1+1,p2+4)+0.004583333333333333*bcoef(p1+2,p2+4)+0.001805555555555556*bcoef(p1+3,p2+4)+0.00006944444444444444*bcoef(p1+4,p2+4)+0.001805555555555556*bcoef(p1+3,p2)+0.00006944444444444444*bcoef(p1+4,p2)+0.001805555555555556*bcoef(p1,p2+1)+0.04694444444444444*bcoef(p1+1,p2+1)+0.1191666666666667*bcoef(p1+2,p2+1)+0.04694444444444444*bcoef(p1+3,p2+1);
3851 coef_mat_buf(1,0) = 0.009027777777777778*bcoef(p1+4,p2+1)-0.003472222222222222*bcoef(p1+1,p2)-0.0003472222222222222*bcoef(p1,p2)-0.02291666666666667*bcoef(p1,p2+2)-0.2291666666666667*bcoef(p1+1,p2+2)+0.2291666666666667*bcoef(p1+3,p2+2)+0.02291666666666667*bcoef(p1+4,p2+2)-0.009027777777777778*bcoef(p1,p2+3)-0.09027777777777778*bcoef(p1+1,p2+3)+0.09027777777777778*bcoef(p1+3,p2+3)+0.009027777777777778*bcoef(p1+4,p2+3)-0.0003472222222222222*bcoef(p1,p2+4)-0.003472222222222222*bcoef(p1+1,p2+4)+0.003472222222222222*bcoef(p1+3,p2+4)+0.0003472222222222222*bcoef(p1+4,p2+4)+0.003472222222222222*bcoef(p1+3,p2)+0.0003472222222222222*bcoef(p1+4,p2)-0.009027777777777778*bcoef(p1,p2+1)-0.09027777777777778*bcoef(p1+1,p2+1)+0.09027777777777778*bcoef(p1+3,p2+1);
3852 coef_mat_buf(2,0) = 0.0006944444444444444*bcoef(p1,p2)+0.001388888888888889*bcoef(p1+1,p2)+0.01805555555555556*bcoef(p1+4,p2+1)+0.04583333333333333*bcoef(p1,p2+2)+0.09166666666666667*bcoef(p1+1,p2+2)-0.275*bcoef(p1+2,p2+2)+0.09166666666666667*bcoef(p1+3,p2+2)+0.04583333333333333*bcoef(p1+4,p2+2)+0.01805555555555556*bcoef(p1,p2+3)+0.03611111111111111*bcoef(p1+1,p2+3)-0.004166666666666667*bcoef(p1+2,p2)-0.1083333333333333*bcoef(p1+2,p2+3)+0.03611111111111111*bcoef(p1+3,p2+3)+0.01805555555555556*bcoef(p1+4,p2+3)+0.0006944444444444444*bcoef(p1,p2+4)+0.001388888888888889*bcoef(p1+1,p2+4)-0.004166666666666667*bcoef(p1+2,p2+4)+0.001388888888888889*bcoef(p1+3,p2+4)+0.0006944444444444444*bcoef(p1+4,p2+4)+0.001388888888888889*bcoef(p1+3,p2)+0.0006944444444444444*bcoef(p1+4,p2)+0.01805555555555556*bcoef(p1,p2+1)+0.03611111111111111*bcoef(p1+1,p2+1)-0.1083333333333333*bcoef(p1+2,p2+1)+0.03611111111111111*bcoef(p1+3,p2+1);
3853 coef_mat_buf(3,0) = 0.001388888888888889*bcoef(p1+1,p2)-0.0006944444444444444*bcoef(p1,p2)+0.01805555555555556*bcoef(p1+4,p2+1)-0.04583333333333333*bcoef(p1,p2+2)+0.09166666666666667*bcoef(p1+1,p2+2)-0.09166666666666667*bcoef(p1+3,p2+2)+0.04583333333333333*bcoef(p1+4,p2+2)-0.01805555555555556*bcoef(p1,p2+3)+0.03611111111111111*bcoef(p1+1,p2+3)-0.03611111111111111*bcoef(p1+3,p2+3)+0.01805555555555556*bcoef(p1+4,p2+3)-0.0006944444444444444*bcoef(p1,p2+4)+0.001388888888888889*bcoef(p1+1,p2+4)-0.001388888888888889*bcoef(p1+3,p2+4)+0.0006944444444444444*bcoef(p1+4,p2+4)-0.001388888888888889*bcoef(p1+3,p2)+0.0006944444444444444*bcoef(p1+4,p2)-0.01805555555555556*bcoef(p1,p2+1)+0.03611111111111111*bcoef(p1+1,p2+1)-0.03611111111111111*bcoef(p1+3,p2+1);
3854 coef_mat_buf(4,0) = 0.0003472222222222222*bcoef(p1,p2)-0.001388888888888889*bcoef(p1+1,p2)+0.009027777777777778*bcoef(p1+4,p2+1)+0.02291666666666667*bcoef(p1,p2+2)-0.09166666666666667*bcoef(p1+1,p2+2)+0.1375*bcoef(p1+2,p2+2)-0.09166666666666667*bcoef(p1+3,p2+2)+0.02291666666666667*bcoef(p1+4,p2+2)+0.009027777777777778*bcoef(p1,p2+3)-0.03611111111111111*bcoef(p1+1,p2+3)+0.002083333333333333*bcoef(p1+2,p2)+0.05416666666666667*bcoef(p1+2,p2+3)-0.03611111111111111*bcoef(p1+3,p2+3)+0.009027777777777778*bcoef(p1+4,p2+3)+0.0003472222222222222*bcoef(p1,p2+4)-0.001388888888888889*bcoef(p1+1,p2+4)+0.002083333333333333*bcoef(p1+2,p2+4)-0.001388888888888889*bcoef(p1+3,p2+4)+0.0003472222222222222*bcoef(p1+4,p2+4)-0.001388888888888889*bcoef(p1+3,p2)+0.0003472222222222222*bcoef(p1+4,p2)+0.009027777777777778*bcoef(p1,p2+1)-0.03611111111111111*bcoef(p1+1,p2+1)+0.05416666666666667*bcoef(p1+2,p2+1)-0.03611111111111111*bcoef(p1+3,p2+1);
3855 coef_mat_buf(5,0) = 0.0003472222222222222*bcoef(p1+1,p2)-0.00006944444444444444*bcoef(p1,p2)-0.009027777777777778*bcoef(p1+4,p2+1)+0.001805555555555556*bcoef(p1+5,p2+1)-0.004583333333333333*bcoef(p1,p2+2)+0.02291666666666667*bcoef(p1+1,p2+2)-0.04583333333333333*bcoef(p1+2,p2+2)+0.04583333333333333*bcoef(p1+3,p2+2)-0.02291666666666667*bcoef(p1+4,p2+2)+0.004583333333333333*bcoef(p1+5,p2+2)-0.001805555555555556*bcoef(p1,p2+3)+0.009027777777777778*bcoef(p1+1,p2+3)-0.0006944444444444444*bcoef(p1+2,p2)-0.01805555555555556*bcoef(p1+2,p2+3)+0.01805555555555556*bcoef(p1+3,p2+3)-0.009027777777777778*bcoef(p1+4,p2+3)+0.001805555555555556*bcoef(p1+5,p2+3)-0.00006944444444444444*bcoef(p1,p2+4)+0.0003472222222222222*bcoef(p1+1,p2+4)-0.0006944444444444444*bcoef(p1+2,p2+4)+0.0006944444444444444*bcoef(p1+3,p2+4)-0.0003472222222222222*bcoef(p1+4,p2+4)+0.00006944444444444444*bcoef(p1+5,p2+4)+0.0006944444444444444*bcoef(p1+3,p2)-0.0003472222222222222*bcoef(p1+4,p2)+0.00006944444444444444*bcoef(p1+5,p2)-0.001805555555555556*bcoef(p1,p2+1)+0.009027777777777778*bcoef(p1+1,p2+1)-0.01805555555555556*bcoef(p1+2,p2+1)+0.01805555555555556*bcoef(p1+3,p2+1);
3856 coef_mat_buf(0,1) = 0.003472222222222222*bcoef(p1,p2+3)-0.009027777777777778*bcoef(p1+1,p2)-0.003472222222222222*bcoef(p1+4,p2+1)-0.0003472222222222222*bcoef(p1,p2)+0.09027777777777778*bcoef(p1+1,p2+3)-0.02291666666666667*bcoef(p1+2,p2)+0.2291666666666667*bcoef(p1+2,p2+3)+0.09027777777777778*bcoef(p1+3,p2+3)+0.003472222222222222*bcoef(p1+4,p2+3)+0.0003472222222222222*bcoef(p1,p2+4)+0.009027777777777778*bcoef(p1+1,p2+4)+0.02291666666666667*bcoef(p1+2,p2+4)+0.009027777777777778*bcoef(p1+3,p2+4)+0.0003472222222222222*bcoef(p1+4,p2+4)-0.009027777777777778*bcoef(p1+3,p2)-0.0003472222222222222*bcoef(p1+4,p2)-0.003472222222222222*bcoef(p1,p2+1)-0.09027777777777778*bcoef(p1+1,p2+1)-0.2291666666666667*bcoef(p1+2,p2+1)-0.09027777777777778*bcoef(p1+3,p2+1);
3857 coef_mat_buf(1,1) = 0.001736111111111111*bcoef(p1,p2)+0.01736111111111111*bcoef(p1+1,p2)-0.01736111111111111*bcoef(p1+4,p2+1)-0.01736111111111111*bcoef(p1,p2+3)-0.1736111111111111*bcoef(p1+1,p2+3)+0.1736111111111111*bcoef(p1+3,p2+3)+0.01736111111111111*bcoef(p1+4,p2+3)-0.001736111111111111*bcoef(p1,p2+4)-0.01736111111111111*bcoef(p1+1,p2+4)+0.01736111111111111*bcoef(p1+3,p2+4)+0.001736111111111111*bcoef(p1+4,p2+4)-0.01736111111111111*bcoef(p1+3,p2)-0.001736111111111111*bcoef(p1+4,p2)+0.01736111111111111*bcoef(p1,p2+1)+0.1736111111111111*bcoef(p1+1,p2+1)-0.1736111111111111*bcoef(p1+3,p2+1);
3858 coef_mat_buf(2,1) = 0.03472222222222222*bcoef(p1,p2+3)-0.006944444444444444*bcoef(p1+1,p2)-0.03472222222222222*bcoef(p1+4,p2+1)-0.003472222222222222*bcoef(p1,p2)+0.06944444444444444*bcoef(p1+1,p2+3)+0.02083333333333333*bcoef(p1+2,p2)-0.2083333333333333*bcoef(p1+2,p2+3)+0.06944444444444444*bcoef(p1+3,p2+3)+0.03472222222222222*bcoef(p1+4,p2+3)+0.003472222222222222*bcoef(p1,p2+4)+0.006944444444444444*bcoef(p1+1,p2+4)-0.02083333333333333*bcoef(p1+2,p2+4)+0.006944444444444444*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)-0.006944444444444444*bcoef(p1+3,p2)-0.003472222222222222*bcoef(p1+4,p2)-0.03472222222222222*bcoef(p1,p2+1)-0.06944444444444444*bcoef(p1+1,p2+1)+0.2083333333333333*bcoef(p1+2,p2+1)-0.06944444444444444*bcoef(p1+3,p2+1);
3859 coef_mat_buf(3,1) = 0.003472222222222222*bcoef(p1,p2)-0.006944444444444444*bcoef(p1+1,p2)-0.03472222222222222*bcoef(p1+4,p2+1)-0.03472222222222222*bcoef(p1,p2+3)+0.06944444444444444*bcoef(p1+1,p2+3)-0.06944444444444444*bcoef(p1+3,p2+3)+0.03472222222222222*bcoef(p1+4,p2+3)-0.003472222222222222*bcoef(p1,p2+4)+0.006944444444444444*bcoef(p1+1,p2+4)-0.006944444444444444*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)+0.006944444444444444*bcoef(p1+3,p2)-0.003472222222222222*bcoef(p1+4,p2)+0.03472222222222222*bcoef(p1,p2+1)-0.06944444444444444*bcoef(p1+1,p2+1)+0.06944444444444444*bcoef(p1+3,p2+1);
3860 coef_mat_buf(4,1) = 0.006944444444444444*bcoef(p1+1,p2)-0.001736111111111111*bcoef(p1,p2)-0.01736111111111111*bcoef(p1+4,p2+1)+0.01736111111111111*bcoef(p1,p2+3)-0.06944444444444444*bcoef(p1+1,p2+3)-0.01041666666666667*bcoef(p1+2,p2)+0.1041666666666667*bcoef(p1+2,p2+3)-0.06944444444444444*bcoef(p1+3,p2+3)+0.01736111111111111*bcoef(p1+4,p2+3)+0.001736111111111111*bcoef(p1,p2+4)-0.006944444444444444*bcoef(p1+1,p2+4)+0.01041666666666667*bcoef(p1+2,p2+4)-0.006944444444444444*bcoef(p1+3,p2+4)+0.001736111111111111*bcoef(p1+4,p2+4)+0.006944444444444444*bcoef(p1+3,p2)-0.001736111111111111*bcoef(p1+4,p2)-0.01736111111111111*bcoef(p1,p2+1)+0.06944444444444444*bcoef(p1+1,p2+1)-0.1041666666666667*bcoef(p1+2,p2+1)+0.06944444444444444*bcoef(p1+3,p2+1);
3861 coef_mat_buf(5,1) = 0.0003472222222222222*bcoef(p1,p2)-0.001736111111111111*bcoef(p1+1,p2)+0.01736111111111111*bcoef(p1+4,p2+1)-0.003472222222222222*bcoef(p1+5,p2+1)-0.003472222222222222*bcoef(p1,p2+3)+0.01736111111111111*bcoef(p1+1,p2+3)+0.003472222222222222*bcoef(p1+2,p2)-0.03472222222222222*bcoef(p1+2,p2+3)+0.03472222222222222*bcoef(p1+3,p2+3)-0.01736111111111111*bcoef(p1+4,p2+3)+0.003472222222222222*bcoef(p1+5,p2+3)-0.0003472222222222222*bcoef(p1,p2+4)+0.001736111111111111*bcoef(p1+1,p2+4)-0.003472222222222222*bcoef(p1+2,p2+4)+0.003472222222222222*bcoef(p1+3,p2+4)-0.001736111111111111*bcoef(p1+4,p2+4)+0.0003472222222222222*bcoef(p1+5,p2+4)-0.003472222222222222*bcoef(p1+3,p2)+0.001736111111111111*bcoef(p1+4,p2)-0.0003472222222222222*bcoef(p1+5,p2)+0.003472222222222222*bcoef(p1,p2+1)-0.01736111111111111*bcoef(p1+1,p2+1)+0.03472222222222222*bcoef(p1+2,p2+1)-0.03472222222222222*bcoef(p1+3,p2+1);
3862 coef_mat_buf(0,2) = 0.0006944444444444444*bcoef(p1,p2)+0.01805555555555556*bcoef(p1+1,p2)+0.001388888888888889*bcoef(p1+4,p2+1)-0.004166666666666667*bcoef(p1,p2+2)-0.1083333333333333*bcoef(p1+1,p2+2)-0.275*bcoef(p1+2,p2+2)-0.1083333333333333*bcoef(p1+3,p2+2)-0.004166666666666667*bcoef(p1+4,p2+2)+0.001388888888888889*bcoef(p1,p2+3)+0.03611111111111111*bcoef(p1+1,p2+3)+0.04583333333333333*bcoef(p1+2,p2)+0.09166666666666667*bcoef(p1+2,p2+3)+0.03611111111111111*bcoef(p1+3,p2+3)+0.001388888888888889*bcoef(p1+4,p2+3)+0.0006944444444444444*bcoef(p1,p2+4)+0.01805555555555556*bcoef(p1+1,p2+4)+0.04583333333333333*bcoef(p1+2,p2+4)+0.01805555555555556*bcoef(p1+3,p2+4)+0.0006944444444444444*bcoef(p1+4,p2+4)+0.01805555555555556*bcoef(p1+3,p2)+0.0006944444444444444*bcoef(p1+4,p2)+0.001388888888888889*bcoef(p1,p2+1)+0.03611111111111111*bcoef(p1+1,p2+1)+0.09166666666666667*bcoef(p1+2,p2+1)+0.03611111111111111*bcoef(p1+3,p2+1);
3863 coef_mat_buf(1,2) = 0.006944444444444444*bcoef(p1+4,p2+1)-0.03472222222222222*bcoef(p1+1,p2)-0.003472222222222222*bcoef(p1,p2)+0.02083333333333333*bcoef(p1,p2+2)+0.2083333333333333*bcoef(p1+1,p2+2)-0.2083333333333333*bcoef(p1+3,p2+2)-0.02083333333333333*bcoef(p1+4,p2+2)-0.006944444444444444*bcoef(p1,p2+3)-0.06944444444444444*bcoef(p1+1,p2+3)+0.06944444444444444*bcoef(p1+3,p2+3)+0.006944444444444444*bcoef(p1+4,p2+3)-0.003472222222222222*bcoef(p1,p2+4)-0.03472222222222222*bcoef(p1+1,p2+4)+0.03472222222222222*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)+0.03472222222222222*bcoef(p1+3,p2)+0.003472222222222222*bcoef(p1+4,p2)-0.006944444444444444*bcoef(p1,p2+1)-0.06944444444444444*bcoef(p1+1,p2+1)+0.06944444444444444*bcoef(p1+3,p2+1);
3864 coef_mat_buf(2,2) = 0.006944444444444444*bcoef(p1,p2)+0.01388888888888889*bcoef(p1+1,p2)+0.01388888888888889*bcoef(p1+4,p2+1)-0.04166666666666667*bcoef(p1,p2+2)-0.08333333333333333*bcoef(p1+1,p2+2)+0.25*bcoef(p1+2,p2+2)-0.08333333333333333*bcoef(p1+3,p2+2)-0.04166666666666667*bcoef(p1+4,p2+2)+0.01388888888888889*bcoef(p1,p2+3)+0.02777777777777778*bcoef(p1+1,p2+3)-0.04166666666666667*bcoef(p1+2,p2)-0.08333333333333333*bcoef(p1+2,p2+3)+0.02777777777777778*bcoef(p1+3,p2+3)+0.01388888888888889*bcoef(p1+4,p2+3)+0.006944444444444444*bcoef(p1,p2+4)+0.01388888888888889*bcoef(p1+1,p2+4)-0.04166666666666667*bcoef(p1+2,p2+4)+0.01388888888888889*bcoef(p1+3,p2+4)+0.006944444444444444*bcoef(p1+4,p2+4)+0.01388888888888889*bcoef(p1+3,p2)+0.006944444444444444*bcoef(p1+4,p2)+0.01388888888888889*bcoef(p1,p2+1)+0.02777777777777778*bcoef(p1+1,p2+1)-0.08333333333333333*bcoef(p1+2,p2+1)+0.02777777777777778*bcoef(p1+3,p2+1);
3865 coef_mat_buf(3,2) = 0.01388888888888889*bcoef(p1+1,p2)-0.006944444444444444*bcoef(p1,p2)+0.01388888888888889*bcoef(p1+4,p2+1)+0.04166666666666667*bcoef(p1,p2+2)-0.08333333333333333*bcoef(p1+1,p2+2)+0.08333333333333333*bcoef(p1+3,p2+2)-0.04166666666666667*bcoef(p1+4,p2+2)-0.01388888888888889*bcoef(p1,p2+3)+0.02777777777777778*bcoef(p1+1,p2+3)-0.02777777777777778*bcoef(p1+3,p2+3)+0.01388888888888889*bcoef(p1+4,p2+3)-0.006944444444444444*bcoef(p1,p2+4)+0.01388888888888889*bcoef(p1+1,p2+4)-0.01388888888888889*bcoef(p1+3,p2+4)+0.006944444444444444*bcoef(p1+4,p2+4)-0.01388888888888889*bcoef(p1+3,p2)+0.006944444444444444*bcoef(p1+4,p2)-0.01388888888888889*bcoef(p1,p2+1)+0.02777777777777778*bcoef(p1+1,p2+1)-0.02777777777777778*bcoef(p1+3,p2+1);
3866 coef_mat_buf(4,2) = 0.003472222222222222*bcoef(p1,p2)-0.01388888888888889*bcoef(p1+1,p2)+0.006944444444444444*bcoef(p1+4,p2+1)-0.02083333333333333*bcoef(p1,p2+2)+0.08333333333333333*bcoef(p1+1,p2+2)-0.125*bcoef(p1+2,p2+2)+0.08333333333333333*bcoef(p1+3,p2+2)-0.02083333333333333*bcoef(p1+4,p2+2)+0.006944444444444444*bcoef(p1,p2+3)-0.02777777777777778*bcoef(p1+1,p2+3)+0.02083333333333333*bcoef(p1+2,p2)+0.04166666666666667*bcoef(p1+2,p2+3)-0.02777777777777778*bcoef(p1+3,p2+3)+0.006944444444444444*bcoef(p1+4,p2+3)+0.003472222222222222*bcoef(p1,p2+4)-0.01388888888888889*bcoef(p1+1,p2+4)+0.02083333333333333*bcoef(p1+2,p2+4)-0.01388888888888889*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)-0.01388888888888889*bcoef(p1+3,p2)+0.003472222222222222*bcoef(p1+4,p2)+0.006944444444444444*bcoef(p1,p2+1)-0.02777777777777778*bcoef(p1+1,p2+1)+0.04166666666666667*bcoef(p1+2,p2+1)-0.02777777777777778*bcoef(p1+3,p2+1);
3867 coef_mat_buf(5,2) = 0.003472222222222222*bcoef(p1+1,p2)-0.0006944444444444444*bcoef(p1,p2)-0.006944444444444444*bcoef(p1+4,p2+1)+0.001388888888888889*bcoef(p1+5,p2+1)+0.004166666666666667*bcoef(p1,p2+2)-0.02083333333333333*bcoef(p1+1,p2+2)+0.04166666666666667*bcoef(p1+2,p2+2)-0.04166666666666667*bcoef(p1+3,p2+2)+0.02083333333333333*bcoef(p1+4,p2+2)-0.004166666666666667*bcoef(p1+5,p2+2)-0.001388888888888889*bcoef(p1,p2+3)+0.006944444444444444*bcoef(p1+1,p2+3)-0.006944444444444444*bcoef(p1+2,p2)-0.01388888888888889*bcoef(p1+2,p2+3)+0.01388888888888889*bcoef(p1+3,p2+3)-0.006944444444444444*bcoef(p1+4,p2+3)+0.001388888888888889*bcoef(p1+5,p2+3)-0.0006944444444444444*bcoef(p1,p2+4)+0.003472222222222222*bcoef(p1+1,p2+4)-0.006944444444444444*bcoef(p1+2,p2+4)+0.006944444444444444*bcoef(p1+3,p2+4)-0.003472222222222222*bcoef(p1+4,p2+4)+0.0006944444444444444*bcoef(p1+5,p2+4)+0.006944444444444444*bcoef(p1+3,p2)-0.003472222222222222*bcoef(p1+4,p2)+0.0006944444444444444*bcoef(p1+5,p2)-0.001388888888888889*bcoef(p1,p2+1)+0.006944444444444444*bcoef(p1+1,p2+1)-0.01388888888888889*bcoef(p1+2,p2+1)+0.01388888888888889*bcoef(p1+3,p2+1);
3868 coef_mat_buf(0,3) = 0.001388888888888889*bcoef(p1+4,p2+1)-0.01805555555555556*bcoef(p1+1,p2)-0.0006944444444444444*bcoef(p1,p2)-0.001388888888888889*bcoef(p1,p2+3)-0.03611111111111111*bcoef(p1+1,p2+3)-0.04583333333333333*bcoef(p1+2,p2)-0.09166666666666667*bcoef(p1+2,p2+3)-0.03611111111111111*bcoef(p1+3,p2+3)-0.001388888888888889*bcoef(p1+4,p2+3)+0.0006944444444444444*bcoef(p1,p2+4)+0.01805555555555556*bcoef(p1+1,p2+4)+0.04583333333333333*bcoef(p1+2,p2+4)+0.01805555555555556*bcoef(p1+3,p2+4)+0.0006944444444444444*bcoef(p1+4,p2+4)-0.01805555555555556*bcoef(p1+3,p2)-0.0006944444444444444*bcoef(p1+4,p2)+0.001388888888888889*bcoef(p1,p2+1)+0.03611111111111111*bcoef(p1+1,p2+1)+0.09166666666666667*bcoef(p1+2,p2+1)+0.03611111111111111*bcoef(p1+3,p2+1);
3869 coef_mat_buf(1,3) = 0.003472222222222222*bcoef(p1,p2)+0.03472222222222222*bcoef(p1+1,p2)+0.006944444444444444*bcoef(p1+4,p2+1)+0.006944444444444444*bcoef(p1,p2+3)+0.06944444444444444*bcoef(p1+1,p2+3)-0.06944444444444444*bcoef(p1+3,p2+3)-0.006944444444444444*bcoef(p1+4,p2+3)-0.003472222222222222*bcoef(p1,p2+4)-0.03472222222222222*bcoef(p1+1,p2+4)+0.03472222222222222*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)-0.03472222222222222*bcoef(p1+3,p2)-0.003472222222222222*bcoef(p1+4,p2)-0.006944444444444444*bcoef(p1,p2+1)-0.06944444444444444*bcoef(p1+1,p2+1)+0.06944444444444444*bcoef(p1+3,p2+1);
3870 coef_mat_buf(2,3) = 0.01388888888888889*bcoef(p1+4,p2+1)-0.01388888888888889*bcoef(p1+1,p2)-0.006944444444444444*bcoef(p1,p2)-0.01388888888888889*bcoef(p1,p2+3)-0.02777777777777778*bcoef(p1+1,p2+3)+0.04166666666666667*bcoef(p1+2,p2)+0.08333333333333333*bcoef(p1+2,p2+3)-0.02777777777777778*bcoef(p1+3,p2+3)-0.01388888888888889*bcoef(p1+4,p2+3)+0.006944444444444444*bcoef(p1,p2+4)+0.01388888888888889*bcoef(p1+1,p2+4)-0.04166666666666667*bcoef(p1+2,p2+4)+0.01388888888888889*bcoef(p1+3,p2+4)+0.006944444444444444*bcoef(p1+4,p2+4)-0.01388888888888889*bcoef(p1+3,p2)-0.006944444444444444*bcoef(p1+4,p2)+0.01388888888888889*bcoef(p1,p2+1)+0.02777777777777778*bcoef(p1+1,p2+1)-0.08333333333333333*bcoef(p1+2,p2+1)+0.02777777777777778*bcoef(p1+3,p2+1);
3871 coef_mat_buf(3,3) = 0.006944444444444444*bcoef(p1,p2)-0.01388888888888889*bcoef(p1+1,p2)+0.01388888888888889*bcoef(p1+4,p2+1)+0.01388888888888889*bcoef(p1,p2+3)-0.02777777777777778*bcoef(p1+1,p2+3)+0.02777777777777778*bcoef(p1+3,p2+3)-0.01388888888888889*bcoef(p1+4,p2+3)-0.006944444444444444*bcoef(p1,p2+4)+0.01388888888888889*bcoef(p1+1,p2+4)-0.01388888888888889*bcoef(p1+3,p2+4)+0.006944444444444444*bcoef(p1+4,p2+4)+0.01388888888888889*bcoef(p1+3,p2)-0.006944444444444444*bcoef(p1+4,p2)-0.01388888888888889*bcoef(p1,p2+1)+0.02777777777777778*bcoef(p1+1,p2+1)-0.02777777777777778*bcoef(p1+3,p2+1);
3872 coef_mat_buf(4,3) = 0.01388888888888889*bcoef(p1+1,p2)-0.003472222222222222*bcoef(p1,p2)+0.006944444444444444*bcoef(p1+4,p2+1)-0.006944444444444444*bcoef(p1,p2+3)+0.02777777777777778*bcoef(p1+1,p2+3)-0.02083333333333333*bcoef(p1+2,p2)-0.04166666666666667*bcoef(p1+2,p2+3)+0.02777777777777778*bcoef(p1+3,p2+3)-0.006944444444444444*bcoef(p1+4,p2+3)+0.003472222222222222*bcoef(p1,p2+4)-0.01388888888888889*bcoef(p1+1,p2+4)+0.02083333333333333*bcoef(p1+2,p2+4)-0.01388888888888889*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)+0.01388888888888889*bcoef(p1+3,p2)-0.003472222222222222*bcoef(p1+4,p2)+0.006944444444444444*bcoef(p1,p2+1)-0.02777777777777778*bcoef(p1+1,p2+1)+0.04166666666666667*bcoef(p1+2,p2+1)-0.02777777777777778*bcoef(p1+3,p2+1);
3873 coef_mat_buf(5,3) = 0.0006944444444444444*bcoef(p1,p2)-0.003472222222222222*bcoef(p1+1,p2)-0.006944444444444444*bcoef(p1+4,p2+1)+0.001388888888888889*bcoef(p1+5,p2+1)+0.001388888888888889*bcoef(p1,p2+3)-0.006944444444444444*bcoef(p1+1,p2+3)+0.006944444444444444*bcoef(p1+2,p2)+0.01388888888888889*bcoef(p1+2,p2+3)-0.01388888888888889*bcoef(p1+3,p2+3)+0.006944444444444444*bcoef(p1+4,p2+3)-0.001388888888888889*bcoef(p1+5,p2+3)-0.0006944444444444444*bcoef(p1,p2+4)+0.003472222222222222*bcoef(p1+1,p2+4)-0.006944444444444444*bcoef(p1+2,p2+4)+0.006944444444444444*bcoef(p1+3,p2+4)-0.003472222222222222*bcoef(p1+4,p2+4)+0.0006944444444444444*bcoef(p1+5,p2+4)-0.006944444444444444*bcoef(p1+3,p2)+0.003472222222222222*bcoef(p1+4,p2)-0.0006944444444444444*bcoef(p1+5,p2)-0.001388888888888889*bcoef(p1,p2+1)+0.006944444444444444*bcoef(p1+1,p2+1)-0.01388888888888889*bcoef(p1+2,p2+1)+0.01388888888888889*bcoef(p1+3,p2+1);
3874 coef_mat_buf(0,4) = 0.0003472222222222222*bcoef(p1,p2)+0.009027777777777778*bcoef(p1+1,p2)-0.001388888888888889*bcoef(p1+4,p2+1)+0.002083333333333333*bcoef(p1,p2+2)+0.05416666666666667*bcoef(p1+1,p2+2)+0.1375*bcoef(p1+2,p2+2)+0.05416666666666667*bcoef(p1+3,p2+2)+0.002083333333333333*bcoef(p1+4,p2+2)-0.001388888888888889*bcoef(p1,p2+3)-0.03611111111111111*bcoef(p1+1,p2+3)+0.02291666666666667*bcoef(p1+2,p2)-0.09166666666666667*bcoef(p1+2,p2+3)-0.03611111111111111*bcoef(p1+3,p2+3)-0.001388888888888889*bcoef(p1+4,p2+3)+0.0003472222222222222*bcoef(p1,p2+4)+0.009027777777777778*bcoef(p1+1,p2+4)+0.02291666666666667*bcoef(p1+2,p2+4)+0.009027777777777778*bcoef(p1+3,p2+4)+0.0003472222222222222*bcoef(p1+4,p2+4)+0.009027777777777778*bcoef(p1+3,p2)+0.0003472222222222222*bcoef(p1+4,p2)-0.001388888888888889*bcoef(p1,p2+1)-0.03611111111111111*bcoef(p1+1,p2+1)-0.09166666666666667*bcoef(p1+2,p2+1)-0.03611111111111111*bcoef(p1+3,p2+1);
3875 coef_mat_buf(1,4) = 0.1041666666666667*bcoef(p1+3,p2+2)-0.01736111111111111*bcoef(p1+1,p2)-0.006944444444444444*bcoef(p1+4,p2+1)-0.01041666666666667*bcoef(p1,p2+2)-0.1041666666666667*bcoef(p1+1,p2+2)-0.001736111111111111*bcoef(p1,p2)+0.01041666666666667*bcoef(p1+4,p2+2)+0.006944444444444444*bcoef(p1,p2+3)+0.06944444444444444*bcoef(p1+1,p2+3)-0.06944444444444444*bcoef(p1+3,p2+3)-0.006944444444444444*bcoef(p1+4,p2+3)-0.001736111111111111*bcoef(p1,p2+4)-0.01736111111111111*bcoef(p1+1,p2+4)+0.01736111111111111*bcoef(p1+3,p2+4)+0.001736111111111111*bcoef(p1+4,p2+4)+0.01736111111111111*bcoef(p1+3,p2)+0.001736111111111111*bcoef(p1+4,p2)+0.006944444444444444*bcoef(p1,p2+1)+0.06944444444444444*bcoef(p1+1,p2+1)-0.06944444444444444*bcoef(p1+3,p2+1);
3876 coef_mat_buf(2,4) = 0.003472222222222222*bcoef(p1,p2)+0.006944444444444444*bcoef(p1+1,p2)-0.01388888888888889*bcoef(p1+4,p2+1)+0.02083333333333333*bcoef(p1,p2+2)+0.04166666666666667*bcoef(p1+1,p2+2)-0.125*bcoef(p1+2,p2+2)+0.04166666666666667*bcoef(p1+3,p2+2)+0.02083333333333333*bcoef(p1+4,p2+2)-0.01388888888888889*bcoef(p1,p2+3)-0.02777777777777778*bcoef(p1+1,p2+3)-0.02083333333333333*bcoef(p1+2,p2)+0.08333333333333333*bcoef(p1+2,p2+3)-0.02777777777777778*bcoef(p1+3,p2+3)-0.01388888888888889*bcoef(p1+4,p2+3)+0.003472222222222222*bcoef(p1,p2+4)+0.006944444444444444*bcoef(p1+1,p2+4)-0.02083333333333333*bcoef(p1+2,p2+4)+0.006944444444444444*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)+0.006944444444444444*bcoef(p1+3,p2)+0.003472222222222222*bcoef(p1+4,p2)-0.01388888888888889*bcoef(p1,p2+1)-0.02777777777777778*bcoef(p1+1,p2+1)+0.08333333333333333*bcoef(p1+2,p2+1)-0.02777777777777778*bcoef(p1+3,p2+1);
3877 coef_mat_buf(3,4) = 0.006944444444444444*bcoef(p1+1,p2)-0.003472222222222222*bcoef(p1,p2)-0.01388888888888889*bcoef(p1+4,p2+1)-0.02083333333333333*bcoef(p1,p2+2)+0.04166666666666667*bcoef(p1+1,p2+2)-0.04166666666666667*bcoef(p1+3,p2+2)+0.02083333333333333*bcoef(p1+4,p2+2)+0.01388888888888889*bcoef(p1,p2+3)-0.02777777777777778*bcoef(p1+1,p2+3)+0.02777777777777778*bcoef(p1+3,p2+3)-0.01388888888888889*bcoef(p1+4,p2+3)-0.003472222222222222*bcoef(p1,p2+4)+0.006944444444444444*bcoef(p1+1,p2+4)-0.006944444444444444*bcoef(p1+3,p2+4)+0.003472222222222222*bcoef(p1+4,p2+4)-0.006944444444444444*bcoef(p1+3,p2)+0.003472222222222222*bcoef(p1+4,p2)+0.01388888888888889*bcoef(p1,p2+1)-0.02777777777777778*bcoef(p1+1,p2+1)+0.02777777777777778*bcoef(p1+3,p2+1);
3878 coef_mat_buf(4,4) = 0.001736111111111111*bcoef(p1,p2)-0.006944444444444444*bcoef(p1+1,p2)-0.006944444444444444*bcoef(p1+4,p2+1)+0.01041666666666667*bcoef(p1,p2+2)-0.04166666666666667*bcoef(p1+1,p2+2)+0.0625*bcoef(p1+2,p2+2)-0.04166666666666667*bcoef(p1+3,p2+2)+0.01041666666666667*bcoef(p1+4,p2+2)-0.006944444444444444*bcoef(p1,p2+3)+0.02777777777777778*bcoef(p1+1,p2+3)+0.01041666666666667*bcoef(p1+2,p2)-0.04166666666666667*bcoef(p1+2,p2+3)+0.02777777777777778*bcoef(p1+3,p2+3)-0.006944444444444444*bcoef(p1+4,p2+3)+0.001736111111111111*bcoef(p1,p2+4)-0.006944444444444444*bcoef(p1+1,p2+4)+0.01041666666666667*bcoef(p1+2,p2+4)-0.006944444444444444*bcoef(p1+3,p2+4)+0.001736111111111111*bcoef(p1+4,p2+4)-0.006944444444444444*bcoef(p1+3,p2)+0.001736111111111111*bcoef(p1+4,p2)-0.006944444444444444*bcoef(p1,p2+1)+0.02777777777777778*bcoef(p1+1,p2+1)-0.04166666666666667*bcoef(p1+2,p2+1)+0.02777777777777778*bcoef(p1+3,p2+1);
3879 coef_mat_buf(5,4) = 0.001736111111111111*bcoef(p1+1,p2)-0.0003472222222222222*bcoef(p1,p2)+0.006944444444444444*bcoef(p1+4,p2+1)-0.001388888888888889*bcoef(p1+5,p2+1)-0.002083333333333333*bcoef(p1,p2+2)+0.01041666666666667*bcoef(p1+1,p2+2)-0.02083333333333333*bcoef(p1+2,p2+2)+0.02083333333333333*bcoef(p1+3,p2+2)-0.01041666666666667*bcoef(p1+4,p2+2)+0.002083333333333333*bcoef(p1+5,p2+2)+0.001388888888888889*bcoef(p1,p2+3)-0.006944444444444444*bcoef(p1+1,p2+3)-0.003472222222222222*bcoef(p1+2,p2)+0.01388888888888889*bcoef(p1+2,p2+3)-0.01388888888888889*bcoef(p1+3,p2+3)+0.006944444444444444*bcoef(p1+4,p2+3)-0.001388888888888889*bcoef(p1+5,p2+3)-0.0003472222222222222*bcoef(p1,p2+4)+0.001736111111111111*bcoef(p1+1,p2+4)-0.003472222222222222*bcoef(p1+2,p2+4)+0.003472222222222222*bcoef(p1+3,p2+4)-0.001736111111111111*bcoef(p1+4,p2+4)+0.0003472222222222222*bcoef(p1+5,p2+4)+0.003472222222222222*bcoef(p1+3,p2)-0.001736111111111111*bcoef(p1+4,p2)+0.0003472222222222222*bcoef(p1+5,p2)+0.001388888888888889*bcoef(p1,p2+1)-0.006944444444444444*bcoef(p1+1,p2+1)+0.01388888888888889*bcoef(p1+2,p2+1)-0.01388888888888889*bcoef(p1+3,p2+1);
3880 coef_mat_buf(0,5) = 0.0003472222222222222*bcoef(p1+4,p2+1)-0.001805555555555556*bcoef(p1+1,p2)-0.00006944444444444444*bcoef(p1,p2)-0.0006944444444444444*bcoef(p1,p2+2)-0.01805555555555556*bcoef(p1+1,p2+2)-0.04583333333333333*bcoef(p1+2,p2+2)-0.01805555555555556*bcoef(p1+3,p2+2)-0.0006944444444444444*bcoef(p1+4,p2+2)+0.0006944444444444444*bcoef(p1,p2+3)+0.01805555555555556*bcoef(p1+1,p2+3)-0.004583333333333333*bcoef(p1+2,p2)+0.04583333333333333*bcoef(p1+2,p2+3)+0.01805555555555556*bcoef(p1+3,p2+3)+0.0006944444444444444*bcoef(p1+4,p2+3)-0.0003472222222222222*bcoef(p1,p2+4)-0.009027777777777778*bcoef(p1+1,p2+4)-0.02291666666666667*bcoef(p1+2,p2+4)-0.009027777777777778*bcoef(p1+3,p2+4)-0.0003472222222222222*bcoef(p1+4,p2+4)-0.001805555555555556*bcoef(p1+3,p2)+0.00006944444444444444*bcoef(p1,p2+5)+0.001805555555555556*bcoef(p1+1,p2+5)+0.004583333333333333*bcoef(p1+2,p2+5)+0.001805555555555556*bcoef(p1+3,p2+5)+0.00006944444444444444*bcoef(p1+4,p2+5)-0.00006944444444444444*bcoef(p1+4,p2)+0.0003472222222222222*bcoef(p1,p2+1)+0.009027777777777778*bcoef(p1+1,p2+1)+0.02291666666666667*bcoef(p1+2,p2+1)+0.009027777777777778*bcoef(p1+3,p2+1);
3881 coef_mat_buf(1,5) = 0.0003472222222222222*bcoef(p1,p2)+0.003472222222222222*bcoef(p1+1,p2)+0.001736111111111111*bcoef(p1+4,p2+1)+0.003472222222222222*bcoef(p1,p2+2)+0.03472222222222222*bcoef(p1+1,p2+2)-0.03472222222222222*bcoef(p1+3,p2+2)-0.003472222222222222*bcoef(p1+4,p2+2)-0.003472222222222222*bcoef(p1,p2+3)-0.03472222222222222*bcoef(p1+1,p2+3)+0.03472222222222222*bcoef(p1+3,p2+3)+0.003472222222222222*bcoef(p1+4,p2+3)+0.001736111111111111*bcoef(p1,p2+4)+0.01736111111111111*bcoef(p1+1,p2+4)-0.01736111111111111*bcoef(p1+3,p2+4)-0.001736111111111111*bcoef(p1+4,p2+4)-0.003472222222222222*bcoef(p1+3,p2)-0.0003472222222222222*bcoef(p1,p2+5)-0.003472222222222222*bcoef(p1+1,p2+5)+0.003472222222222222*bcoef(p1+3,p2+5)+0.0003472222222222222*bcoef(p1+4,p2+5)-0.0003472222222222222*bcoef(p1+4,p2)-0.001736111111111111*bcoef(p1,p2+1)-0.01736111111111111*bcoef(p1+1,p2+1)+0.01736111111111111*bcoef(p1+3,p2+1);
3882 coef_mat_buf(2,5) = 0.003472222222222222*bcoef(p1+4,p2+1)-0.001388888888888889*bcoef(p1+1,p2)-0.0006944444444444444*bcoef(p1,p2)-0.006944444444444444*bcoef(p1,p2+2)-0.01388888888888889*bcoef(p1+1,p2+2)+0.04166666666666667*bcoef(p1+2,p2+2)-0.01388888888888889*bcoef(p1+3,p2+2)-0.006944444444444444*bcoef(p1+4,p2+2)+0.006944444444444444*bcoef(p1,p2+3)+0.01388888888888889*bcoef(p1+1,p2+3)+0.004166666666666667*bcoef(p1+2,p2)-0.04166666666666667*bcoef(p1+2,p2+3)+0.01388888888888889*bcoef(p1+3,p2+3)+0.006944444444444444*bcoef(p1+4,p2+3)-0.003472222222222222*bcoef(p1,p2+4)-0.006944444444444444*bcoef(p1+1,p2+4)+0.02083333333333333*bcoef(p1+2,p2+4)-0.006944444444444444*bcoef(p1+3,p2+4)-0.003472222222222222*bcoef(p1+4,p2+4)-0.001388888888888889*bcoef(p1+3,p2)+0.0006944444444444444*bcoef(p1,p2+5)+0.001388888888888889*bcoef(p1+1,p2+5)-0.004166666666666667*bcoef(p1+2,p2+5)+0.001388888888888889*bcoef(p1+3,p2+5)+0.0006944444444444444*bcoef(p1+4,p2+5)-0.0006944444444444444*bcoef(p1+4,p2)+0.003472222222222222*bcoef(p1,p2+1)+0.006944444444444444*bcoef(p1+1,p2+1)-0.02083333333333333*bcoef(p1+2,p2+1)+0.006944444444444444*bcoef(p1+3,p2+1);
3883 coef_mat_buf(3,5) = 0.0006944444444444444*bcoef(p1,p2)-0.001388888888888889*bcoef(p1+1,p2)+0.003472222222222222*bcoef(p1+4,p2+1)+0.006944444444444444*bcoef(p1,p2+2)-0.01388888888888889*bcoef(p1+1,p2+2)+0.01388888888888889*bcoef(p1+3,p2+2)-0.006944444444444444*bcoef(p1+4,p2+2)-0.006944444444444444*bcoef(p1,p2+3)+0.01388888888888889*bcoef(p1+1,p2+3)-0.01388888888888889*bcoef(p1+3,p2+3)+0.006944444444444444*bcoef(p1+4,p2+3)+0.003472222222222222*bcoef(p1,p2+4)-0.006944444444444444*bcoef(p1+1,p2+4)+0.006944444444444444*bcoef(p1+3,p2+4)-0.003472222222222222*bcoef(p1+4,p2+4)+0.001388888888888889*bcoef(p1+3,p2)-0.0006944444444444444*bcoef(p1,p2+5)+0.001388888888888889*bcoef(p1+1,p2+5)-0.001388888888888889*bcoef(p1+3,p2+5)+0.0006944444444444444*bcoef(p1+4,p2+5)-0.0006944444444444444*bcoef(p1+4,p2)-0.003472222222222222*bcoef(p1,p2+1)+0.006944444444444444*bcoef(p1+1,p2+1)-0.006944444444444444*bcoef(p1+3,p2+1);
3884 coef_mat_buf(4,5) = 0.001388888888888889*bcoef(p1+1,p2)-0.0003472222222222222*bcoef(p1,p2)+0.001736111111111111*bcoef(p1+4,p2+1)-0.003472222222222222*bcoef(p1,p2+2)+0.01388888888888889*bcoef(p1+1,p2+2)-0.02083333333333333*bcoef(p1+2,p2+2)+0.01388888888888889*bcoef(p1+3,p2+2)-0.003472222222222222*bcoef(p1+4,p2+2)+0.003472222222222222*bcoef(p1,p2+3)-0.01388888888888889*bcoef(p1+1,p2+3)-0.002083333333333333*bcoef(p1+2,p2)+0.02083333333333333*bcoef(p1+2,p2+3)-0.01388888888888889*bcoef(p1+3,p2+3)+0.003472222222222222*bcoef(p1+4,p2+3)-0.001736111111111111*bcoef(p1,p2+4)+0.006944444444444444*bcoef(p1+1,p2+4)-0.01041666666666667*bcoef(p1+2,p2+4)+0.006944444444444444*bcoef(p1+3,p2+4)-0.001736111111111111*bcoef(p1+4,p2+4)+0.001388888888888889*bcoef(p1+3,p2)+0.0003472222222222222*bcoef(p1,p2+5)-0.001388888888888889*bcoef(p1+1,p2+5)+0.002083333333333333*bcoef(p1+2,p2+5)-0.001388888888888889*bcoef(p1+3,p2+5)+0.0003472222222222222*bcoef(p1+4,p2+5)-0.0003472222222222222*bcoef(p1+4,p2)+0.001736111111111111*bcoef(p1,p2+1)-0.006944444444444444*bcoef(p1+1,p2+1)+0.01041666666666667*bcoef(p1+2,p2+1)-0.006944444444444444*bcoef(p1+3,p2+1);
3885 coef_mat_buf(5,5) = 0.00006944444444444444*bcoef(p1,p2)-0.0003472222222222222*bcoef(p1+1,p2)-0.001736111111111111*bcoef(p1+4,p2+1)+0.0003472222222222222*bcoef(p1+5,p2+1)+0.0006944444444444444*bcoef(p1,p2+2)-0.003472222222222222*bcoef(p1+1,p2+2)+0.006944444444444444*bcoef(p1+2,p2+2)-0.006944444444444444*bcoef(p1+3,p2+2)+0.003472222222222222*bcoef(p1+4,p2+2)-0.0006944444444444444*bcoef(p1+5,p2+2)-0.0006944444444444444*bcoef(p1,p2+3)+0.003472222222222222*bcoef(p1+1,p2+3)+0.0006944444444444444*bcoef(p1+2,p2)-0.006944444444444444*bcoef(p1+2,p2+3)+0.006944444444444444*bcoef(p1+3,p2+3)-0.003472222222222222*bcoef(p1+4,p2+3)+0.0006944444444444444*bcoef(p1+5,p2+3)+0.0003472222222222222*bcoef(p1,p2+4)-0.001736111111111111*bcoef(p1+1,p2+4)+0.003472222222222222*bcoef(p1+2,p2+4)-0.003472222222222222*bcoef(p1+3,p2+4)+0.001736111111111111*bcoef(p1+4,p2+4)-0.0003472222222222222*bcoef(p1+5,p2+4)-0.0006944444444444444*bcoef(p1+3,p2)-0.00006944444444444444*bcoef(p1,p2+5)+0.0003472222222222222*bcoef(p1+1,p2+5)-0.0006944444444444444*bcoef(p1+2,p2+5)+0.0006944444444444444*bcoef(p1+3,p2+5)-0.0003472222222222222*bcoef(p1+4,p2+5)+0.00006944444444444444*bcoef(p1+5,p2+5)+0.0003472222222222222*bcoef(p1+4,p2)-0.00006944444444444444*bcoef(p1+5,p2)-0.0003472222222222222*bcoef(p1,p2+1)+0.001736111111111111*bcoef(p1+1,p2+1)-0.003472222222222222*bcoef(p1+2,p2+1)+0.003472222222222222*bcoef(p1+3,p2+1);
3887 return coef_mat_buf;
3891 template <
typename T_container>
3895 coef_mat_precompute_ptr = std::make_shared<Array2D<container>>(this->
A_ptr->height(), this->
A_ptr->width(),
container(6,6));
3897 auto bcoef = this->get_bspline_mat_ptr(A);
3900 this->calc_coef_mat((*coef_mat_precompute_ptr)(p1,p2), *bcoef, p1 + this->bcoef_border - 2, p2 + this->bcoef_border - 2);
3906 template <
typename T_container>
3908 if (A.height() != A.width()) {
3909 throw std::invalid_argument(
"Attempted to perform backward substitution using an Array with size " + A.size_2D_string() +
3910 ". Array must be square.");
3912 if (A.height() != b.height() || b.width() != 1) {
3913 throw std::invalid_argument(
"Attempted to perform backward substitution using an A with size " + A.size_2D_string() +
3914 " and a b with size " + b.size_2D_string() +
".");
3919 std::copy(b.get_pointer(), b.get_pointer() + b.size(), x_buf.get_pointer());
3921 if ((std::abs((*A_factored_ptr)(p,p)) < std::numeric_limits<value_type>::epsilon())) {
3925 x_buf(p) /= (*A_factored_ptr)(p,p);
3929 x_buf(p1) -= x_buf(p) * (*A_factored_ptr)(p1,p);
3933 if ((std::abs((*A_factored_ptr)(0,0)) < std::numeric_limits<value_type>::epsilon())) {
3936 x_buf(0) /= (*A_factored_ptr)(0,0);
3942 template <
typename T_container>
3944 if (A.height() != A.width()) {
3945 throw std::invalid_argument(
"Attempted to perform forward substitution using an Array with size " + A.size_2D_string() +
3946 ". Array must be square.");
3948 if (A.height() != b.height() || b.width() != 1) {
3949 throw std::invalid_argument(
"Attempted to perform forward substitution using an A with size " + A.size_2D_string() +
3950 " and a b with size " + b.size_2D_string() +
".");
3955 std::copy(b.get_pointer(), b.get_pointer() + b.size(), x_buf.get_pointer());
3957 if ((std::abs((*A_factored_ptr)(p,p)) < std::numeric_limits<value_type>::epsilon())) {
3961 x_buf(p) /= (*A_factored_ptr)(p,p);
3964 for (
difference_type p1 = p + 1; p1 < A_factored_ptr->height(); ++p1) {
3965 x_buf(p1) -= x_buf(p) * (*A_factored_ptr)(p1,p);
3969 if ((std::abs((*A_factored_ptr)(
last,
last)) < std::numeric_limits<value_type>::epsilon())) {
3979 template <
typename T_container>
3981 if (A.height() != A.width()) {
3982 throw std::invalid_argument(
"Attempted to perform LU decomposition on array of size " + A.size_2D_string() +
3983 ". Array must be square.");
3991 (*this->piv_ptr)(p) = mu;
4001 if (std::abs((*this->
A_factored_ptr)(p,p)) < std::numeric_limits<value_type>::epsilon()) {
4022 template <
typename T_container>
4024 if ((*this->A_factored_ptr).height() != b.height() || b.width() != 1) {
4025 throw std::invalid_argument(
"Attempted to solve LU decomposition using b of size " + b.size_2D_string() +
4026 " on an Array of size " + (*this->A_factored_ptr).size_2D_string() +
".");
4035 std::copy(b.get_pointer(), b.get_pointer() + b.size(), this->x_buf.get_pointer());
4039 this->x_buf((*piv_ptr)(p)) = this->x_buf(p);
4040 this->x_buf(p) = buf;
4044 this->x_buf(p1) -= this->x_buf(p) * (*this->A_factored_ptr)(p1,p);
4048 return this->backward_sub((*this->A_factored_ptr), this->x_buf);
4052 template <
typename T_container>
4062 c(p2) = dot(A(
all,p2),A(
all,p2));
4067 while (std::abs(tau) >= std::numeric_limits<value_type>::epsilon() && tau > 0) {
4069 (*piv_ptr)(rank-1) = k;
4083 (*beta_ptr)(rank-1) = h_pair.second;
4086 (*this->
A_factored_ptr)({rank,
last},rank-1) = h_pair.first({1,(*this->A_factored_ptr).height()-rank});
4092 if (rank < min_dim) {
4094 k = find(c == tau, rank);
4100 full_rank = (rank == (*this->A_factored_ptr).width());
4103 template <
typename T_container>
4105 if (x.width() != 1) {
4106 throw std::invalid_argument(
"Attempted to get household vector for vector of size " + x.size_2D_string() +
4107 ". Vector must be a column vector.");
4114 if (std::abs(sigma) > std::numeric_limits<value_type>::epsilon()) {
4115 value_type mu = std::sqrt(std::pow(x(0),2) + sigma);
4116 if (std::abs(x(0)) < std::numeric_limits<value_type>::epsilon() || x(0) < 0) {
4119 v(0) = -sigma / (x(0) + mu);
4121 beta = 2 * std::pow(v(0),2) / (sigma + std::pow(v(0),2));
4125 return {std::move(v),std::move(beta)};
4128 template <
typename T_container>
4130 if ((*this->A_factored_ptr).height() != b.height() || b.width() != 1) {
4131 throw std::invalid_argument(
"Attempted to solve QR decomposition using b of size " + b.size_2D_string() +
4132 " on an Array of size " + (*this->A_factored_ptr).size_2D_string() +
".");
4144 v({p+1,
last}) = (*this->A_factored_ptr)({p+1,
last},p);
4149 this->backward_sub((*this->A_factored_ptr)({0,rank-1},{0,rank-1}),y({0,rank-1}));
4152 this->x_buf({rank,
last}) = 0;
4154 value_type buf = this->x_buf((*piv_ptr)(p-1));
4155 this->x_buf((*piv_ptr)(p-1)) = this->x_buf(p-1);
4156 this->x_buf(p-1) = buf;
4163 template <
typename T_container>
4165 if (A.height() != A.width()) {
4166 throw std::invalid_argument(
"Attempted to perform Cholesky decomposition on array of size " + A.size_2D_string() +
4167 ". Array must be square.");
4185 if ((*this->
A_factored_ptr)(p,p) >= std::numeric_limits<value_type>::epsilon()) {
4204 template <
typename T_container>
4207 throw std::invalid_argument(
"Attempted to solve Cholesky decomposition with a non positive definite matrix.");
4209 if ((*this->A_factored_ptr).height() != b.height() || b.width() != 1) {
4210 throw std::invalid_argument(
"Attempted to solve Cholesky decomposition using b of size " + b.size_2D_string() +
4211 " on an Array of size " + (*this->A_factored_ptr).size_2D_string() +
".");
4219 return this->backward_sub((*this->A_factored_ptr), this->forward_sub((*this->A_factored_ptr), b));
4224template <
typename T =
double,
typename T_alloc = std::allocator<T>>
4229 throw std::invalid_argument(
"Attempted to create identity matrix with a size of: " + std::to_string(n) +
".");
4233 for (difference_type p = 0; p < n; ++p) {