/******************************************************** * ██████╗ ██████╗████████╗██╗ * ██╔════╝ ██╔════╝╚══██╔══╝██║ * ██║ ███╗██║ ██║ ██║ * ██║ ██║██║ ██║ ██║ * ╚██████╔╝╚██████╗ ██║ ███████╗ * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ * Geophysical Computational Tools & Library (GCTL) * * Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn) * * GCTL is distributed under a dual licensing scheme. You can redistribute * it and/or modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation, either version 2 * of the License, or (at your option) any later version. You should have * received a copy of the GNU Lesser General Public License along with this * program. If not, see . * * If the terms and conditions of the LGPL v.2. would prevent you from using * the GCTL, please consider the option to obtain a commercial license for a * fee. These licenses are offered by the GCTL's original author. As a rule, * licenses are provided "as-is", unlimited in time for a one time fee. Please * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget * to include some description of your company and the realm of its activities. * Also add information on how to contact you by electronic and paper mail. ******************************************************/ #include "extrapolate.h" void gctl::cosine_extrapolate_1d(double *in_arr, int in_num, double *out_arr, int out_num, double end_value) { if (in_arr == nullptr || out_arr == nullptr) throw domain_error("Invalid pointer. Thrown by gctl::cosine_extrapolate_1d(...)"); if (in_num >= out_num) throw invalid_argument("Incompatible output array size. Thrown by gctl::cosine_extrapolate_1d(...)"); if (std::isnan(end_value)) { end_value = 0.5*(in_arr[0] + in_arr[in_num-1]); } // 确定左右侧外插值数量 int l_num = floor((out_num - in_num)/2.0); int r_num = ceil((out_num - in_num)/2.0); for (int i = 0; i < l_num; i++) { out_arr[i] = end_value + (in_arr[0]-end_value)*cos(-0.5*GCTL_Pi*(i-l_num)/l_num); } for (int i = 0; i < in_num; i++) { out_arr[l_num+i] = in_arr[i]; } for (int i = 0; i < r_num; i++) { out_arr[l_num+in_num+i] = end_value + (in_arr[in_num-1]-end_value)*cos(0.5*GCTL_Pi*(i+1)/r_num); } return; } void gctl::cosine_extrapolate_1d(const array &in_arr, array &out_arr, double end_val, int out_len) { if (in_arr.empty()) throw domain_error("The input array is empty. From gctl::cosine_extrapolate_1d(...)"); int ex_len; if (out_len < 0) { ex_len = pow(2,ceil(log(in_arr.size())/log(2))+1); } else ex_len = out_len; out_arr.resize(ex_len, 0.0); cosine_extrapolate_1d(in_arr.get(), in_arr.size(), out_arr.get(), ex_len, end_val); return; } void gctl::cosine_extrapolate_2d(const _2d_matrix &in_arr, _2d_matrix &out_arr, int &ori_row, int &ori_col, double end_valx, double end_valy, int out_lenx, int out_leny) { if (in_arr.empty()) throw domain_error("The input array is empty. From gctl::cosine_extrapolate_2d(...)"); int M_ex, N_ex; int M = in_arr.row_size(); int N = in_arr.col_size(); if (out_lenx < 0) { N_ex = pow(2,ceil(log(N)/log(2))+1); } else N_ex = out_lenx; if (out_leny < 0) { M_ex = pow(2,ceil(log(M)/log(2))+1); } else M_ex = out_leny; if (M >= M_ex || N >= N_ex) throw logic_error("Invalid input or output array size. From gctl::cosine_extrapolate_2d(...)"); out_arr.resize(M_ex, N_ex, 0.0); gctl::array endmean(N); gctl::array endmean2(M_ex); int M_ex_up = floor((M_ex - M)/2.0); int M_ex_down = ceil((M_ex - M)/2.0); int N_ex_left = floor((N_ex - N)/2.0); int N_ex_right = ceil((N_ex - N)/2.0); ori_row = M_ex_down; ori_col = N_ex_left; if (std::isnan(end_valy)) { for (int i = 0; i < N; i++) endmean[i] = (in_arr[0][i]+in_arr[M-1][i])/2.0; } else { for (int i = 0; i < N; i++) endmean[i] = end_valy; } for (int i = 0; i < N; i++) { for (int k = 0; k < M_ex_down; k++) { out_arr[k][N_ex_left+i] = endmean[i] + (in_arr[0][i]-endmean[i])*cos(-0.5*GCTL_Pi*(k-M_ex_down)/M_ex_down); } for (int k = 0; k < M; k++) { out_arr[M_ex_down+k][N_ex_left+i] = in_arr[k][i]; } for (int k = 0; k < M_ex_up; k++) { out_arr[M_ex_down+M+k][N_ex_left+i] = endmean[i] + (in_arr[M-1][i]-endmean[i])*cos(0.5*GCTL_Pi*(k+1)/M_ex_up); } } if (std::isnan(end_valx)) { for (int i = 0; i < M_ex; i++) endmean2[i] = (out_arr[i][N_ex_left]+out_arr[i][N_ex_left+N-1])/2.0; } else { for (int i = 0; i < M_ex; i++) endmean2[i] = end_valx; } for (int i = 0; i < M_ex; i++) { for (int k = 0; k < N_ex_left; k++) { out_arr[i][k] = endmean2[i] + (out_arr[i][N_ex_left]-endmean2[i])*cos(-0.5*GCTL_Pi*(k-N_ex_left)/N_ex_left); } for (int k = 0; k < N_ex_right; k++) { out_arr[i][N_ex_left+N+k] = endmean2[i] + (out_arr[i][N_ex_left+N-1]-endmean2[i])*cos(0.5*GCTL_Pi*(k+1)/N_ex_right); } } return; } void gctl::cosine_extrapolate_2d(const array &in_arr, array &out_arr, int in_row, int in_col, int &out_row, int &out_col, int &ori_row, int &ori_col, double end_valx, double end_valy) { if (in_arr.empty()) throw runtime_error("The input array is empty. Thrown by gctl::cosine_extrapolate_2d(...)"); if (in_row*in_col != in_arr.size()) throw invalid_argument("Incompatible array sizes. Thrown by gctl::cosine_extrapolate_2d(...)"); int M = in_row; int N = in_col; if (out_col < 0) out_col = pow(2,ceil(log(N)/log(2))+1); else if (out_col <= in_col) { throw invalid_argument("Incompatible output column size. Thrown by gctl::cosine_extrapolate_2d(...)"); } int N_ex = out_col; if (out_row < 0) out_row = pow(2,ceil(log(M)/log(2))+1); else if (out_row <= in_row) { throw invalid_argument("Incompatible output row size. Thrown by gctl::cosine_extrapolate_2d(...)"); } int M_ex = out_row; out_arr.resize(M_ex*N_ex, 0.0); gctl::array endmean(N); gctl::array endmean2(M_ex); int M_ex_up = floor((M_ex - M)/2.0); int M_ex_down = ceil((M_ex - M)/2.0); int N_ex_left = floor((N_ex - N)/2.0); int N_ex_right = ceil((N_ex - N)/2.0); ori_row = M_ex_down; ori_col = N_ex_left; if (std::isnan(end_valy)) { for (int i = 0; i < N; i++) endmean[i] = (in_arr[i]+in_arr[(M-1)*N+i])/2.0; } else { for (int i = 0; i < N; i++) endmean[i] = end_valy; } for (int i = 0; i < N; i++) { for (int k = 0; k < M_ex_down; k++) { out_arr[k*N_ex+N_ex_left+i] = endmean[i] + (in_arr[i]-endmean[i])*cos(-0.5*GCTL_Pi*(k-M_ex_down)/M_ex_down); } for (int k = 0; k < M; k++) { out_arr[(M_ex_down+k)*N_ex+N_ex_left+i] = in_arr[k*N+i]; } for (int k = 0; k < M_ex_up; k++) { out_arr[(M_ex_down+M+k)*N_ex+N_ex_left+i] = endmean[i] + (in_arr[(M-1)*N+i]-endmean[i])*cos(0.5*GCTL_Pi*(k+1)/M_ex_up); } } if (std::isnan(end_valx)) { for (int i = 0; i < M_ex; i++) endmean2[i] = (out_arr[i*N_ex+N_ex_left]+out_arr[i*N_ex+N_ex_left+N-1])/2.0; } else { for (int i = 0; i < M_ex; i++) endmean2[i] = end_valx; } for (int i = 0; i < M_ex; i++) { for (int k = 0; k < N_ex_left; k++) { out_arr[i*N_ex+k] = endmean2[i] + (out_arr[i*N_ex+N_ex_left]-endmean2[i])*cos(-0.5*GCTL_Pi*(k-N_ex_left)/N_ex_left); } for (int k = 0; k < N_ex_right; k++) { out_arr[i*N_ex+N_ex_left+N+k] = endmean2[i] + (out_arr[i*N_ex+N_ex_left+N-1]-endmean2[i])*cos(0.5*GCTL_Pi*(k+1)/N_ex_right); } } return; } void gctl::zeros_extrapolate_2d(const array &in_arr, array &out_arr, int in_row, int in_col, int &out_row, int &out_col, int &ori_row, int &ori_col) { if (in_arr.empty()) throw runtime_error("The input array is empty. Thrown by gctl::zeros_extrapolate_2d(...)"); if (in_row*in_col != in_arr.size()) throw invalid_argument("Incompatible array sizes. Thrown by gctl::zeros_extrapolate_2d(...)"); int M = in_row; int N = in_col; if (out_col < 0) out_col = pow(2,ceil(log(N)/log(2))+1); else if (out_col <= in_col) { throw invalid_argument("Incompatible output column size. Thrown by gctl::zeros_extrapolate_2d(...)"); } int N_ex = out_col; if (out_row < 0) out_row = pow(2,ceil(log(M)/log(2))+1); else if (out_row <= in_row) { throw invalid_argument("Incompatible output row size. Thrown by gctl::zeros_extrapolate_2d(...)"); } int M_ex = out_row; out_arr.resize(M_ex*N_ex, 0.0); int M_ex_down = ceil((M_ex - M)/2.0); int N_ex_left = floor((N_ex - N)/2.0); ori_row = M_ex_down; ori_col = N_ex_left; for (int i = 0; i < M; i++) { for (int j = 0; j < N; j++) { out_arr[(M_ex_down+i)*N_ex+N_ex_left+j] = in_arr[i*N+j]; } } return; }