精华内容
下载资源
问答
  • 针对现有彩色图像插值算法的实时性及可靠性不能兼备的问题,提出了一种高效的实时彩色图像缩放算法,算法基于Lanczos核生成可应用于整个目标图像的核查找表,并且目标图像所有像素的插值仅涉及定点小数运算,计算...
  • Lanczos逼近 Lanczos逼近是Cornelius Lanczos于1964年发布的一种用于数值计算伽玛函数的方法。它是替代实用的以固定精度计算更流行的Stirling逼近的一种实用方法。 Lanczos近似由以下公式组成: 用于伽玛函数,...
  • 某公司实现lanczos插值的定点化实现代码,可以高效实现任意大小的resize,同时保证图像质量
  • 对求解大规模稀疏Hamilton矩阵特征问题的辛Lanczos算法给出了舍人误差分析。分析表明辛Lanczos算法在无中断时,保Hamilton结构的限制没有破坏非对称Lanczos算法 的本质特性。本文还讨论了辛Lanczos算法计算出的辛...
  • 带有 Lanczos 平滑的“单次旋转”元胞自动机演示。 另见3D版本! 您也可以观看。 是一个非常简单但丰富的可逆元胞自动机。 与著名不同,它保留了活细胞的总数。 这允许跟踪每个细胞在其进化过程中穿过场时的轨迹。 ...
  • 本代码使用多相位插值法实现图像缩放,实际上在4 x 4领域大小内进行多相位插值和三次插值几乎是一样的,只是对应插值函数值...多相位插值法是通过对输出点对应原图中的领域进行Lanczos2 函数移相插值来产生输出点的。
  • Lanczos双对角化方法中近似奇异三元组的改进
  • Lanczos 方法 Lanczos 方法是求解大规模稀疏矩阵特征值问题最常用的方法之一,其主要思想是将对称矩阵三对角化,通过对三对角矩阵进行特征分解,从而得到原对称矩阵的特征值与特征向量。 由于三对角矩阵是稀疏矩阵...

    1. 特征值分解(EVD)

    如果 A A A 是一个 m × m m\times m m×m实对称矩阵 A = A T A=A^T A=AT) ,如果存在 m m m 维列向量 q q q 和实数 λ \lambda λ 满足 A q = λ q Aq=\lambda q Aq=λq,则称 q q q A A A特征向量 λ \lambda λ A A A特征值

    如果我们求出 A A A m m m 个特征值与特征向量,并将所有特征向量标准化,则可以得到标准正交矩阵 Q = [ q 1 , q 2 , ⋯   , q m ] , Q ∈ R m × m Q=[q_1, q_2, \cdots, q_m], Q\in \mathbb{R}^{m\times m} Q=[q1,q2,,qm],QRm×m,与对角矩阵 Σ = d i a g [ λ 1 , λ 2 , ⋯   , λ m ] \Sigma = diag[\lambda_1, \lambda_2, \cdots, \lambda_m] Σ=diag[λ1,λ2,,λm],他们满足如下形式:
    A = Q Σ Q T = Q [ λ 1 ⋯ ⋯ ⋯ ⋯ λ 2 ⋯ ⋯ ⋯ ⋯ ⋱ ⋯ ⋯ ⋯ ⋯ λ m ] Q T A = Q\Sigma Q^T= Q\left[ \begin{matrix} \lambda_1 & \cdots & \cdots & \cdots\\ \cdots & \lambda_2 & \cdots & \cdots\\ \cdots & \cdots & \ddots & \cdots\\ \cdots & \cdots & \cdots & \lambda_m\\ \end{matrix} \right]Q^T A=QΣQT=Qλ1λ2λmQT

    numpy 中对一个矩阵进行特征值分解使用 np.linalg.eig 函数,从下例可以看出,该函数求解得到的特征值大小是随机排列的,如果需要前 n n n 个大的特征值与特征向量,需要进一步进行处理:

    #创建一个对称矩阵
    sym_matrix = np.array([[1,2,3,4],
                          [2,5,6,7],
                          [3,6,8,9],
                          [4,7,9,10]])
    #对矩阵求解特征值(se)与特征向量(q)
    se, q = np.linalg.eig(sym_matrix)    
    print('se:')
    print(se)
    print('q:')
    print(q)
    """
    se:
    [24.06253512 -0.80548492  0.5580362   0.18491359]
    q:
    [[-0.22593827 -0.72531367 -0.56047033 -0.32976505]
     [-0.44322186 -0.31846973  0.77633831 -0.3153256 ]
     [-0.57278788 -0.14246073 -0.05844028  0.805111  ]
     [-0.6514755   0.59346614 -0.28241204 -0.37897369]]
    """
    

    在特征值分解中,对要分解的矩阵有着较高的要求,它需要被分解的矩阵 A A A 为实对称矩阵,但是现实中我们遇到的矩阵并不都是实对称的,对于一个 m × n m\times n m×n 的矩阵,如果我们想分解为如上的形式,我们需要使用 SVD分解。

    2. 奇异值分解(SVD)

    对于一个 m × n m\times n m×n 的实数矩阵 A A A,我们定义其奇异值分解(SVD)的形式为:
    A = U Σ V T A=U\Sigma V^T A=UΣVT

    其中 U ∈ R m × m U\in \mathbb{R}^{m\times m} URm×m V ∈ R n × n V\in \mathbb{R}^{n\times n} VRn×n 为单位正交矩阵,即有 U U T = I UU^T=I UUT=I V V T = I VV^T=I VVT=I U U U 为左奇异矩阵, V V V 为右奇异矩阵, Σ ∈ R m × n \Sigma \in \mathbb{R}^{m\times n} ΣRm×n 仅在主对角线上有值,我们称它们为奇异值。分解示意图为:
    在这里插入图片描述

    奇异值分解的求解通过下式:
    A A T = U Σ V T V Σ T U T = U Σ Σ T U T AA^T =U\Sigma V^TV\Sigma^TU^T =U\Sigma \Sigma^TU^T AAT=UΣVTVΣTUT=UΣΣTUT

    A T A = V Σ T U T U Σ V T = V Σ T Σ V T A^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T ATA=VΣTUTUΣVT=VΣTΣVT

    因为 A A T AA^T AAT A A T AA^T AAT 是对称矩阵,则通过对上面两个等式求它们的特征值和特征向量,就可以求出 U , V , Σ U,V,\Sigma U,V,Σ

    我们可以使用 numpy 中的 np.linalg.svd 函数对矩阵求解 SVD :

    #创建一个5×3的矩阵
    T = np.array([[16,25,35],
                  [45,23,6],
                  [7,12,9],
                  [10,18,12],
                  [15,14,15]])
    #进行SVD分解      
    u,ss,vh = np.linalg.svd(T, full_matrices=False)
    
    print('u:')
    print(u)
    
    print('ss:')
    print(ss)
    
    print('vh')
    print(vh)
    
    #输出 u * ss * vh
    print('u*ss*vh:')
    print(u.dot(np.diag(ss)).dot(vh))
    """
    u:
    [[-0.57781685 -0.62450998  0.36719022]
     [-0.63112018  0.74837707  0.09682015]
     [-0.21931497 -0.12213601 -0.43502282]
     [-0.31447256 -0.15727956 -0.7906354 ]
     [-0.34759596 -0.10131625  0.20358785]]
    ss:
    [72.50580717 29.74834802  6.47639687]
    vh
    [[-0.65566262 -0.58091643 -0.48233041]
     [ 0.66347606 -0.13833268 -0.73529829]
     [ 0.3604248  -0.80212229  0.47612372]]
    u*ss*vh:
    [[16. 25. 35.]
     [45. 23.  6.]
     [ 7. 12.  9.]
     [10. 18. 12.]
     [15. 14. 15.]]
    """
    

    3. QR 分解

    一个矩阵 A ∈ R m × n ( m ≤ n ) A\in \mathbb{R}^{m\times n}(m \leq n) ARm×n(mn) 可以被分解为 A = Q R A=QR A=QR 的形式,其中:

    • Q ∈ R m × m Q\in \mathbb{R}^{m\times m} QRm×m 是一个正交矩阵,即 Q T Q = Q Q T = I Q^TQ=QQ^T=I QTQ=QQT=I
    • R = [ R ^ 0 ] ∈ R m × n R=\left[\begin{matrix} \hat{R} \\ 0 \end{matrix}\right] \in \mathbb{R}^{m\times n} R=[R^0]Rm×n R ^ ∈ R m × m \hat{R} \in \mathbb{R}^{m\times m} R^Rm×m是上三角矩阵

    numpy 中,进行 QR 分解使用 np.linalg.qr 函数:

    A = np.array([[16,25,35],
                  [45,23,6],
                  [7,12,9],
                  [10,18,12],
                  [15,14,15]])
    q,r = np.linalg.qr(A)
    print('q:')
    print(q)
    
    print('r:')
    print(r)
    
    print('q*r')
    print(q.dot(r)) 
    """
    q:
    [[-0.31051867  0.63947628  0.59444254]
     [-0.87333376 -0.44331977 -0.09225387]
     [-0.13585192  0.33011618 -0.35331103]
     [-0.19407417  0.51220924 -0.66985902]
     [-0.29111125  0.15232425  0.25414072]]
    r:
    [[-51.52669211 -37.04875904 -24.02638224]
     [  0.          21.10425203  31.12417123]
     [  0.           0.          12.84596908]]
    q*r
    [[16. 25. 35.]
     [45. 23.  6.]
     [ 7. 12.  9.]
     [10. 18. 12.]
     [15. 14. 15.]]
    """
    

    4. Lanczos 方法

    Lanczos 方法是求解大规模稀疏矩阵特征值问题最常用的方法之一,其主要思想是将对称矩阵三对角化,通过对三对角矩阵进行特征分解,从而得到原对称矩阵的特征值与特征向量。

    由于三对角矩阵是稀疏矩阵,因此对其进行特征值求解的时间复杂度相对是很低的。

    4.1 Lanczos 方法思路

    设矩阵 A ∈ R n × n A \in \mathbb{R}^{n\times n} ARn×n,如果 A = A T A=A^T A=AT,那么一定存在 正交矩阵 Q n Q_n Qn ,使得 T = Q n T A Q n \boldsymbol{T=Q_n^TAQ_n} T=QnTAQn,其中 T T T 是对称三对角矩阵。

    可以将以上分解中的 Q Q Q 写为:
    Q = [ q 1 , q 2 , ⋯   , q n ] Q=[q_1, q_2,\cdots,q_n] Q=[q1,q2,,qn]

    其中 q i q_i qi Q Q Q 的列向量。

    T T T 记为:
    T = [ α 1 β 1 0 0 0 0 ⋯ 0 0 β 1 α 2 β 2 0 0 0 ⋯ 0 0 0 β 2 α 3 β 3 0 0 ⋯ 0 0 0 0 β 3 α 4 β 4 0 ⋯ 0 0 0 0 0 β 4 α 5 β 5 ⋯ 0 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 0 0 0 0 ⋯ α n − 1 β n − 1 0 0 0 0 0 0 ⋯ β n − 1 α n ] T= \left[ \begin{matrix} \alpha_1 & \beta_1 & 0 & 0 & 0& 0& \cdots & 0 & 0\\ \beta_1 & \alpha_2 & \beta_2 & 0 & 0& 0& \cdots & 0 & 0\\ 0 & \beta_2 &\alpha_3 &\beta_3 & 0& 0 & \cdots & 0 & 0\\ 0 & 0 & \beta_3 &\alpha_4 &\beta_4 & 0& \cdots & 0 & 0\\ 0 & 0 & 0 & \beta_4 &\alpha_5 &\beta_5 & \cdots & 0 & 0\\ \vdots & \vdots &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots &\vdots\\ 0 & 0 & 0 & 0 & 0 & 0 & \cdots & \alpha_{n-1} & \beta_{n-1}\\ 0 & 0 & 0 & 0 & 0 & 0 & \cdots & \beta_{n-1} & \alpha_n \end{matrix} \right] T=α1β100000β1α2β200000β2α3β300000β3α4β400000β4α5000000β50000000αn1βn100000βn1αn

    从上面我们可以知道,对于一个矩阵 A A A ,我们如果想要求其三对角矩阵,要求的是 T T T 中的 α 、 β \alpha、\beta αβ 参数,以及 Q Q Q 矩阵,因此通过如下的变形,我们求其中的变量。

    由于 Q Q Q 是正交的,我们可以将求解三对角矩阵的式子写为:
    A Q = Q T AQ=QT AQ=QT

    对左边的式子,我们取 Q Q Q 的第一列与 A A A 作乘积,得到左边结果的第一列为 A q 1 Aq_1 Aq1。同时我们取右式 T T T 的第一列与 Q Q Q 作乘积,得到 α 1 q 1 + β 1 q 2 \alpha_1q_1 + \beta_1q_2 α1q1+β1q2 ,即
    A q 1 = α 1 q 1 + β 1 q 2 Aq_1=\alpha_1q_1 + \beta_1q_2 Aq1=α1q1+β1q2

    同理,有:
    A q 2 = β 1 q 1 + α 2 q 2 + β 2 q 3 Aq_2 = \beta_1q_1 + \alpha_2q_2 + \beta_2q_3 Aq2=β1q1+α2q2+β2q3

    A q 3 = β 2 q 2 + α 3 q 3 + β 3 q 4 Aq_3 = \beta_2q_2 + \alpha_3q_3 + \beta_3q_4 Aq3=β2q2+α3q3+β3q4

    由此可以得出:
    A q i = β i − 1 q i − 1 + α i q i + β i q i + 1 , i = 1 , 2 , 3 , ⋯   , n Aq_i = \beta_{i-1}q_{i-1} + \alpha_iq_i + \beta_iq_{i+1},i=1,2,3,\cdots,n Aqi=βi1qi1+αiqi+βiqi+1i=1,2,3,,n

    由于 β 0 , β n , q 0 , q n + 1 \beta_0, \beta_n, q_0, q_{n+1} β0,βn,q0,qn+1 都是不存在的,所以我们令 β 0 q 0 = β n q n + 1 = 0 \beta_0q_0 = \beta_nq_{n+1} = 0 β0q0=βnqn+1=0

    我们将上式两边同时左乘以 q i T q_i^T qiT,由于 Q Q Q 是正交的,所以其列向量互相正交,且 q i T q i = 1 q_i^Tq_i=1 qiTqi=1,我们可以得到:
    q i T A q i = α i q_i^T A q_i = \alpha_i qiTAqi=αi

    通过上式可以看出,如果我们知道了 q i q_i qi ,就能求出 α i \alpha_i αi,这里就是求解 α \boldsymbol{\alpha} α 的方法。

    我们再对等式两边移项可以得到:
    β i q i + 1 = A q i − β i − 1 q i − 1 − α i q i \beta_i q_{i+1} = Aq_i - \beta_{i-1}q_{i-1} - \alpha_iq_i βiqi+1=Aqiβi1qi1αiqi

    上式左侧是一个向量,由于 Q Q Q 是正交矩阵,所以 q i + 1 q_{i+1} qi+1 是一个归一化的单位向量,我们令
    r i = β i q i + 1 = A q i − β i − 1 q i − 1 − α i q i r_i=\beta_i q_{i+1} = Aq_i - \beta_{i-1}q_{i-1} - \alpha_iq_i ri=βiqi+1=Aqiβi1qi1αiqi

    β i = ∥ β i q i + 1 ∥ 2 \beta_i = \|\beta_i q_{i+1}\|_2 βi=βiqi+12 ,即

    β i = ∥ r i ∥ 2 = ∥ β i q i + 1 ∥ 2 = ∥ A q i − β i − 1 q i − 1 − α i q i ∥ 2 \beta_i = \|r_i\|_2=\|\beta_i q_{i+1}\|_2 = \|Aq_i - \beta_{i-1}q_{i-1} - \alpha_iq_i\|_2 βi=ri2=βiqi+12=Aqiβi1qi1αiqi2

    从上式可以看出,如果我们知道了 α i 、 q i 、 β i − 1 、 q i − 1 \alpha_i、q_i、\beta_{i-1}、q_{i-1} αiqiβi1qi1 就能求出 β i \beta_i βi,又因为 q q q 都是单位向量,所以我们还能通过 ∥ β i q i + 1 ∥ 2 \|\beta_i q_{i+1}\|_2 βiqi+12求出 q i + 1 q_{i+1} qi+1,这样就能循环求解了。

    我们从初始条件看,如果我们知道了 q 1 q_1 q1 的值,就能根据上面求 α \boldsymbol{\alpha} α 的方法求出 α 1 \alpha_1 α1,因为我们假设 β 0 q 0 = 0 \beta_0q_0=0 β0q0=0,所以我们就能通过 β 0 q 0 、 α 1 、 q 1 \beta_0q_0、\alpha_1、q_1 β0q0α1q1 求出 β 1 \beta_1 β1,再通过令 q i + 1 = r i ∥ r i ∥ 2 q_{i+1}=\frac{r_i}{\|r_i\|_2} qi+1=ri2ri,这时 q i + 1 q_{i+1} qi+1 就是一个单位向量,我们就得到了下一轮循环的 q q q。因此所有的重点就在 q 1 q_1 q1 ,在实际求解中,我们可以随机初始化一个单位向量作为 q 1 q_1 q1,由此我们可以得到计算 A A A 的三对角矩阵的流程:

    1. 随机初始化一个单位向量 q 1 q_1 q1
    2. 计算 α 1 = q 1 T A q 1 \alpha_1=q_1^TAq_1 α1=q1TAq1
    3. r 1 = A q 1 − α 1 q 1 − β 0 q 0 r_1 = Aq_1-\alpha_1q_1-\beta_0q_0 r1=Aq1α1q1β0q0
    4. β 1 = ∥ r 1 ∥ 2 \beta_1 = \|r_1\|_2 β1=r12
    5. q 2 = r 1 β 1 q_2=\frac{r_1}{\beta_1} q2=β1r1
    6. 重复以上步骤,依次求出 { q 1 , q 2 , ⋯   , q n } \{q_1,q_2,\cdots,q_n\} {q1,q2,,qn} α 1 , β 1 , α 2 , β 2 , ⋯   , α n \alpha_1, \beta_1, \alpha_2, \beta_2, \cdots, \alpha_n α1,β1,α2,β2,,αn

    此即 Lanczos 迭代,其中 q i q_i qi 为 Lanczos 向量,若其中某一步 β i = 0 \beta_i=0 βi=0 ,则此时得到的 T j T_j Tj 的特征值都将是 A A A 的特征值。

    到这里,我们就得到了 T = Q T A Q T=Q^TAQ T=QTAQ,如果此时对 T T T 进行特征分解得到 T = W S W T T=WSW^T T=WSWT,所以有
    A = ( Q W ) S ( W T Q T ) A=(QW)S(W^TQ^T) A=(QW)S(WTQT)

    我们得到了通过将矩阵 A A A 变换为三对角矩阵,从而求解特征值与特征向量的方法。将矩阵转换为三对角矩阵的 python 实现为:

    #生成一个对称矩阵
    sym_matrix = np.array([[1,2,3,4],
                          [2,5,6,7],
                          [3,6,8,9],
                          [4,7,9,10]])
    Q = []
    A = []
    B = []
    sum_dim = sym_matrix.shape[0]
        
    #随机初始化q1
    q1 = np.random.random((sum_dim,1))
    q1 = q1/np.linalg.norm(q1)
    Q.append(q1)
    
    #计算a0
    a0 = (q1.transpose()[0]).dot(sym_matrix).dot(q1)[0]
    A.append(a0)
    
    #计算r0
    r = sym_matrix.dot(q1) - a0*q1
    
    #计算b0
    b0 = np.linalg.norm(r)
    B.append(b0)
    
    #循环迭代,得到A、B和Q
    for i in range(1,sum_dim):
        qi = r/B[i-1]
        Q.append(qi)
    
        ai = qi.transpose()[0].dot(sym_matrix).dot(qi)[0]
        A.append(ai)
    
        r = sym_matrix.dot(qi) - ai*qi - B[i-1]*Q[i-1]
        if i < sum_dim - 1:
            b_next = np.linalg.norm(r)
            B.append(b_next)
    
    #组成三对角矩阵
    L = np.diag(A) + np.diag(B,1) + np.diag(B,-1)
    
    print(L)
    """
    [[10.99996429  6.7727395   0.          0.          0.        ]
     [ 6.7727395  11.4584147   5.03926679  0.          0.        ]
     [ 0.          5.03926679 -6.36314172  4.98605503  0.        ]
     [ 0.          0.          4.98605503  1.00512957  1.8902924 ]
     [ 0.          0.          0.          1.8902924   0.89963316]]
    """
    

    4.2 三对角矩阵的特征值分解

    三对角化后的矩阵 T T T 有三种常用的方法进行特征分解:

    1. QR 法,即对 T T T 进行 QR 分解,令 T = Q R T=QR T=QR,再令 T ′ = R Q T'=RQ T=RQ,如此反复,直至收敛。其中产生的 Q Q Q 矩阵的连乘得到的结果就是 T T T 的特征矩阵。
    2. 分治法。这个方法将 T T T 划均匀分成两块三对角阵和一个秩为 1 的 2 × 2 2\times 2 2×2 的块状方阵的和。分别对两个三对角阵分解,然后合并。合并时对特征向量作秩1的修正。
    3. MRRR 方法。这个方法理论上是最快的方式,复杂度为 O ( l 2 ) O(l^2) O(l2),主要采用逆迭代来加快收敛。

    4.3 利用 QR 算法求解特征值

    给定一个实对称矩阵 X X X,假设其通过特征值分解得到 X = P S P T X=PSP^T X=PSPT,其中 P P P 为正交阵, S S S 是对角阵,利用 QR 算法求解其特征值分解的算法为:

    X 1 = X X_1=X X1=X
    P = E P=E P=E
    f o r   k = 1 , 2 , ⋯ for \ k=1,2,\cdots for k=1,2,
    X k = Q k × R k \quad X_k=Q_k\times R_k Xk=Qk×Rk
    X k + 1 = R k × Q k \quad X_{k+1}=R_k\times Q_k Xk+1=Rk×Qk
    P = P × Q k \quad P=P\times Q_k P=P×Qk

    这个算法很简单,每次迭代只需要做一次 QR 分解和一个矩阵乘法即可,用 python 实现为:

    #创建对称矩阵
    sym_matrix = np.array([[1,2,3,4],
                          [2,5,6,7],
                          [3,6,8,9],
                          [4,7,9,10]])
    dim = sym_matrix.shape
    #初始化P为单位矩阵
    P = np.eye(dim[0])
    X_n = sym_matrix
    i = 1
    while i<20:   #迭代次数可以根据需要改变
        i += 1
        (Q, R) = np.linalg.qr(X_n)  #QR分解
        X_n = R.dot(Q)
        P = P.dot(Q)
    
    print("QR分解求解特征矩阵")
    print(P)
    
    se, q = np.linalg.eig(sym_matrix)
    print("特征值函数求解特征矩阵:")
    print(q)
    
    """
    QR分解求解特征矩阵
    [[-0.22593827 -0.72442125 -0.56162332 -0.32976505]
     [-0.44322186 -0.31970419  0.77583077 -0.3153256 ]
     [-0.57278788 -0.1423676  -0.05866681  0.805111  ]
     [-0.6514755   0.5939146  -0.28146771 -0.37897369]]
    特征值函数求解特征矩阵:
    [[-0.22593827 -0.72531367 -0.56047033 -0.32976505]
     [-0.44322186 -0.31846973  0.77633831 -0.3153256 ]
     [-0.57278788 -0.14246073 -0.05844028  0.805111  ]
     [-0.6514755   0.59346614 -0.28241204 -0.37897369]]
    """
    
    展开全文
  • 滤波程序Lanczos_Filter

    2013-07-18 18:54:33
    算出了Lanczos_Filter滤波的权重系数
  • Lanczos算法计算伽玛函数的自然对数值,数值计算,计算结果较为精确
  • 图像Lanczos3滤波C实现——优化

    千次阅读 2019-03-30 14:57:40
    上一篇博文给出了图像Lanczos3滤波的直观实现,但是整个算法实现过于低效,其原因在于对于每个插值点,都需要重新计算Lanczos系数,本文参考https://github.com/richgel999/imageresampler中的实现,将其C++代码转换...

    上一篇博文给出了图像Lanczos3滤波的直观实现,但是整个算法实现过于低效,其原因在于对于每个插值点,都需要重新计算Lanczos系数,本文参考https://github.com/richgel999/imageresampler中的实现,将其C++代码转换成C实现。
    优化思路:每一行其插值的x位置相同,每一列其插值的y位置相同,因此可以将x位置的插值系数做成table,将y位置的插值系数也做成table,在插值某个点 ( x , y ) (x,y) (x,y)时只需查表,即可得系数。
    代码实现
    数据结构准备:

    #define LANCZOS_SIZE 3.0f
    
    //Lanzcos系数
    typedef struct lanc_coef_t
    {
    	float weight;//lanzcos核函数值
    	int   pos; //位置,可以为x或y
    }lanc_coef;
    
    //Lanzcos系数向量,存储所有的系数
    typedef struct lanc_coef_vec_t
    {
    	lanc_coef* coefs;
    	int size;
    }lanc_coef_vec;
    
    //某个插值点所有Lanzcos系数
    typedef struct lanc_coef_node_t
    {
    	lanc_coef* p;//某个插值点Lanzcos系数的起始点地址
    	int n;//非0Lanzcos系数的个数
    }lanc_coef_node;
    
    //所有的插值点Lanzcos系数向量,所有插值点所有Lanzcos系数组成的向量
    typedef struct lanc_coef_node_vec_t
    {
    	lanc_coef_node*  nds;
    	int size;
    }lanc_coef_node_vec;
    //某个插值点的位置和领域的范围
    typedef struct lanc_bound_t
    {
    	float c;
    	int l;
    	int r;
    }lanc_bound;
    
    //所有插值点范围组成的向量,total为所有点系数和
    typedef struct lanc_bound_vec_t
    {
    	lanc_bound* lbs;
    	int size;
    	int total;
    	float scale;
    }lanc_bound_vec;
    
    //插值结果缓冲区
    typedef struct real_vec_t
    {
    	float* data;
    	int	  size;
    }real_vec;
    
    #define MAX_SCAN_BUF_SIZE 16384
    
    typedef struct lanc_sampler_t
    {
    	int intermediate_w;
    
    	int resample_src_w;
    	int resample_src_h;
    	int resample_dst_w;
    	int resample_dst_h;
    
    	float* dst_buf;
    	float* tmp_buf;
    
    	int delay_x_resample;
    
    	int* src_y_count;
    	uint8_t* src_y_flag;
    
    	int scan_buf_y[MAX_SCAN_BUF_SIZE];
    	float* scan_buf_l[MAX_SCAN_BUF_SIZE];
    
    	int cur_src_y;
    	int cur_dst_y;
    }lanc_sampler;
    
    

    相关数据结构创建与释放:

    real_vec*  create_real_vec(int size);
    void  destroy_real_vec(real_vec** _rv);
    
    real_vec*  create_real_vec(int size)
    {
    	if (size < 0)return NULL;
    	real_vec* rv = (real_vec*)calloc(1, sizeof(real_vec));
    	if (!rv) return rv;
    	rv->data = (float*)calloc(1, size*sizeof(float));
    	if (!rv->data)
    	{
    		destroy_real_vec(&rv);
    		return NULL;
    	}
    	rv->size = size;
    	return rv;
    }
    
    void  destroy_real_vec(real_vec** _rv)
    {
    	if (!_rv)return;
    	real_vec* rv = *_rv;
    	if (!rv) return;
    	if (rv->data)
    	{
    		free(rv->data);
    		rv->data = NULL;
    	}
    	free(rv);
    	*_rv = NULL;
    }
    
    void destroy_lanc_coef_vev(lanc_coef_vec** _lcv)
    {
    	if (!_lcv)return;
    	lanc_coef_vec* lcv = *_lcv;
    	if (lcv)
    	{
    		if (lcv->coefs)
    			free(lcv->coefs);
    		free(lcv);
    	}
    	*_lcv = NULL;
    }
    
    lanc_coef_vec* create_lanc_coef_vec(int size)
    {
    	lanc_coef_vec* lcv = NULL;
    	lcv = (lanc_coef_vec*)calloc(1, sizeof(lanc_coef_vec));
    	if (!lcv)return lcv;
    	lanc_coef* coefs = (lanc_coef*)calloc(size, sizeof(lanc_coef));
    	if (!coefs)
    	{
    		destroy_lanc_coef_vev(&lcv);
    		return NULL;
    	}
    	lcv->coefs = coefs;
    	lcv->size = size;
    	return lcv;
    }
    
    void destroy_lanc_bound_vec(lanc_bound_vec** _lbv)
    {
    	if (!_lbv)return;
    	lanc_bound_vec* lbv = *_lbv;
    	if (lbv)
    	{
    		if (lbv->lbs)
    			free(lbv->lbs);
    		free(lbv);
    	}
    	*_lbv = NULL;
    }
    
    lanc_bound_vec* create_bound_vec(int size)
    {
    	lanc_bound_vec* lbv = NULL;
    	lbv = (lanc_bound_vec*)calloc(1, sizeof(lanc_bound_vec));
    	if (!lbv)return lbv;
    	lanc_bound* bs = (lanc_bound*)calloc(size, sizeof(lanc_bound));
    	if (!bs)
    	{
    		destroy_lanc_bound_vec(&lbv);
    		return NULL;
    	}
    	lbv->lbs = bs;
    	lbv->size = size;
    	return lbv;
    }
    
    void destroy_lanc_coef_node_vev(lanc_coef_node_vec** _cnv)
    {
    	if (!_cnv)return;
    	lanc_coef_node_vec* cnv = *_cnv;
    	if (cnv)
    	{
    		if (cnv->nds)
    			free(cnv->nds);
    		free(cnv);
    	}
    	*_cnv = NULL;
    }
    
    lanc_coef_node_vec* create_lanc_coef_node_vec(int size)
    {
    	lanc_coef_node_vec* cnv = NULL;
    	cnv = (lanc_coef_node_vec*)calloc(1, sizeof(lanc_coef_node_vec));
    	if (!cnv)return cnv;
    	lanc_coef_node* nds = (lanc_coef_node*)calloc(size, sizeof(lanc_coef_node));
    	if (!nds)
    	{
    		destroy_lanc_coef_vev(&cnv);
    		return NULL;
    	}
    	cnv->nds = nds;
    	cnv->size = size;
    	return cnv;
    }
    

    Lanzcos3滤波核函数实现:

    static double sinc(double x)
    {
    	x = (x * M_PI);
    
    	if ((x < 0.01f) && (x > -0.01f))
    		return 1.0f + x*x*(-1.0f / 6.0f + x*x*1.0f / 120.0f);
    
    	return sin(x) / x;
    }
    
    static float clean(double t)
    {
    	const float EPSILON = .0000125f;
    	if (fabs(t) < EPSILON)
    		return 0.0f;
    	return (float)t;
    }
    
    static float lanczos_coef(float t)
    {
    	if (t < 0.0f)
    		t = -t;
    	if (t < 3.0f)
    		return clean(sinc(t) * sinc(t / 3.0f));
    	else
    		return (0.0f);
    }
    

    生成x或y方向上每个插值点的邻域范围,c表示当前插值点,l表示领域的左侧,对y方向来说是上侧,r表示领域右侧,对y方向是下侧.

    static void make_lan_bounds(lanc_bound_vec* bv, float scale)
    {
    	int total = 0;
    	int i;
    	float hw;
    	int size;
    	int l, r;
    	float c;
    	lanc_bound bound;
    	if (!bv)return;
    	if (scale > 1.0f)
    		hw = LANCZOS_SIZE;
    	else
    	    hw = LANCZOS_SIZE / scale;
    	size = bv->size;
    	for (i = 0; i < size; i++)
    	{
    		c = (float)i / scale;
    		l = (int)floor(c - hw);
    		r = (int)ceil(c + hw);
    
    		bound.c = c;
    		bound.l = l;
    		bound.r = r;
    
    		total += (r - l + 1);
    		bv->lbs[i] = bound;
    	}
    	bv->total = total;
    	if (scale >= 1.0f)
    		scale = 1.0f;
    	bv->scale = scale;
    	return;
    }
    

    计算每个插值点Lanzcos核函数值:

    static void make_lanc_coef_nodes(lanc_coef_vec* cv, lanc_bound_vec* bv, lanc_coef_node_vec* nv, int img_size)
    {
    	if (!cv || !bv || !nv)return;
    	float c;
    	int l, r;
    	lanc_bound bound;
    	lanc_coef_node node;
    	float s, norm;
    	float coef;
    	int i, j;
    	int k, max_k;
    	float max_coef;
    	int pos;
    	int size = bv->size;
    	lanc_coef* next = cv->coefs;
    	for (i = 0; i < size; i++)
    	{
    		bound = bv->lbs[i];
    		l = bound.l;
    		r = bound.r;
    		c = bound.c;
    		s = 0.0;
    		//累加系数值,需要归一化
    		for (j = l; j <= r; j++)
    		{
    			coef = lanczos_coef((c - (float)j)*bv->scale);
    			s += coef;
    		}
    		norm = 1.0f / s;
    
    		node = nv->nds[i];
    		node.n = 0;
    		node.p = next;
    		next += (r - l + 1);
    
    		s = 0.0f;
    		max_coef = FLT_MIN;
    		max_k = -1;
    		//计算归一化系数,s为系数和,找到系数最大值和位置,如果s不等于1,将1-s值加上最大值,作为修正,为了省去没必要的乘法,只是保留非0值
    		for (j = l; j <= r; j++)
    		{
    			coef = lanczos_coef((c - (float)j)*bv->scale);
    			coef *= norm;
    			if (coef == 0.0f)
    				continue;
    			pos = j;
    			if (j < 0)pos = 0;
    			if (j >= img_size) 
    				pos = img_size - 1;
    			k = node.n++;
    			node.p[k].pos = pos;
    			node.p[k].weight = coef;
    
    			s += coef;
    			if (max_coef < coef)
    			{
    				max_coef = coef;
    				max_k = k;
    			}
    		}
    		if (max_k == -1 || node.n == 0)
    		{
    			printf("make_lanc_coef_nodes fail\n");
    			return;
    		}
    		if (s != 1.0f)
    		{
    			node.p[max_k].weight += (1.0 - s);
    		}
    		nv->nds[i] = node;
    	}
    }
    

    整合上面的方法,生成x和y方向的插值系数:

    void make_lanc_coef(uint8_2d* src, uint8_2d* dst, lanc_coef_node_vec** _xcnv, lanc_coef_node_vec** _ycnv)
    {
    	lanc_coef_vec* xcv = NULL;
    	lanc_coef_vec* ycv = NULL;
    	lanc_bound_vec* xbv = NULL;
    	lanc_bound_vec* ybv = NULL;
    	lanc_coef_node_vec* xcnv = NULL;
    	lanc_coef_node_vec* ycnv = NULL;
    	if (!src || !dst)return;
    	float scale;
    	int dw, dh;
    	int sw, sh;
    	dw = dst->cols;
    	dh = dst->rows;
    	sw = src->cols;
    	sh = src->rows;
    	scale = (float)dw / (float)sw;
    	xbv = create_bound_vec(dw);
    	ybv = create_bound_vec(dh);
    	xcnv = create_lanc_coef_node_vec(dw);
    	ycnv = create_lanc_coef_node_vec(dh);
    	if (!xbv || !ybv || !xcnv || !ycnv)
    	{
    		destroy_lanc_bound_vec(&xbv);
    		destroy_lanc_bound_vec(&ybv);
    		destroy_lanc_coef_node_vev(&xcnv);
    		destroy_lanc_coef_node_vev(&ycnv);
    		return;
    	}
    	make_lan_bounds(xbv, scale);
    	scale = (float)dh / (float)sh;
    	make_lan_bounds(ybv, scale);
    	xcv = create_lanc_coef_vec(xbv->total);
    	ycv = create_lanc_coef_vec(ybv->total);
    	if (!xcv || !ycv)
    	{
    		destroy_lanc_bound_vec(&xbv);
    		destroy_lanc_bound_vec(&ybv);
    		destroy_lanc_coef_vev(&xcv);
    		destroy_lanc_coef_vev(&ycv);
    		return;
    	}
    	make_lanc_coef_nodes(xcv, xbv, xcnv,sw);
    	make_lanc_coef_nodes(ycv, ybv, ycnv,sh);
    
    	destroy_lanc_bound_vec(&xbv);
    	destroy_lanc_bound_vec(&ybv);
    	*_xcnv = xcnv;
    	*_ycnv = ycnv;
    }
    
    

    创建采样器:

    void destroy_lanc_sampler(lanc_sampler** _sampler)
    {
    	if (!_sampler)return;
    	lanc_sampler* sampler = *_sampler;
    	if (sampler)
    	{
    		if (sampler->src_y_count)
    			free(sampler->src_y_count);
    		if (sampler->src_y_flag)
    			free(sampler->src_y_flag);
    		if (sampler->dst_buf)
    			free(sampler->dst_buf);
    		if (sampler->tmp_buf)
    			free(sampler->tmp_buf);
    		free(sampler);
    	}
    	*_sampler = NULL;
    }
    //边界检查,可以是x或y
    static inline int resampler_range_check(int v, int h) { (void)h; assert((v >= 0) && (v < h)); return v; }
    //计算操作数量,程序规模估算用
    inline int count_ops(lanc_coef_node_vec* cnv, int k)
    {
    	int i, t = 0;
    	for (i = 0; i < k; i++)
    		t += cnv->nds[i].n;
    	return (t);
    }
    //创建采样器
    lanc_sampler* create_lanc_sampler(uint8_2d* src, uint8_2d* dst, lanc_coef_node_vec* xcnv, lanc_coef_node_vec* ycnv)
    {
    	if (!src || !dst)return;
    	lanc_sampler*  sampler = NULL;
    	sampler = (lanc_sampler*)calloc(1, sizeof(lanc_sampler));
    	if (!sampler)return sampler;
    	sampler->resample_src_w = src->cols;
    	sampler->resample_src_h = src->rows;
    	sampler->resample_dst_w = dst->cols;
    	sampler->resample_dst_h = dst->rows;
    	sampler->dst_buf = (float*)calloc(dst->cols, sizeof(float));
    	if (!sampler->dst_buf)
    	{
    		destroy_lanc_sampler(&sampler);
    		return NULL;
    	}
    	sampler->src_y_count = (int*)calloc(sampler->resample_src_h, sizeof(int));
    	if (!sampler->src_y_count)
    	{
    		destroy_lanc_sampler(&sampler);
    		return NULL;
    	}
    	sampler->src_y_flag = (uint8_t*)calloc(sampler->resample_src_h, sizeof(uint8_t));
    	if (!sampler->src_y_flag)
    	{
    		destroy_lanc_sampler(&sampler);
    		return NULL;
    	}
    	int i, j;
    	for (i = 0; i < sampler->resample_dst_h; i++)
    		for (j = 0; j < ycnv->nds[i].n; j++)
    			sampler->src_y_count[resampler_range_check(ycnv->nds[i].p[j].pos, sampler->resample_src_h)]++;
    	for (i = 0; i < MAX_SCAN_BUF_SIZE; i++)
    	{
    		sampler->scan_buf_y[i] = -1;
    		sampler->scan_buf_l[i] = NULL;
    	}
    	   {
    		   // Determine which axis to resample first by comparing the number of multiplies required
    		   // for each possibility.
    		   int x_ops = count_ops(xcnv, dst->cols);
    		   int y_ops = count_ops(ycnv, dst->rows);
    
    		   // Hack 10/2000: Weight Y axis ops a little more than X axis ops.
    		   // (Y axis ops use more cache resources.)
    		   int xy_ops = x_ops * sampler->resample_src_h +
    			   (4 * y_ops * sampler->resample_dst_w) / 3;
    
    		   int yx_ops = (4 * y_ops * sampler->resample_src_w) / 3 +
    			   x_ops * sampler->resample_dst_h;
    #if TEST
    		   printf("src: %i %i\n", sampler->resample_src_w, sampler->resample_src_h);
    		   printf("dst: %i %i\n", sampler->resample_dst_w, sampler->resample_dst_h);
    		   printf("x_ops: %i\n", x_ops);
    		   printf("y_ops: %i\n", y_ops);
    		   printf("xy_ops: %i\n", xy_ops);
    		   printf("yx_ops: %i\n", yx_ops);
    #endif
              //根据x、y方向的操作计算,决定x先插值或y方向先插值
    		   // Now check which resample order is better. In case of a tie, choose the order
    		   // which buffers the least amount of data.
    		   if ((xy_ops > yx_ops) ||
    			   ((xy_ops == yx_ops) && (sampler->resample_src_w < sampler->resample_dst_w))
    			   )
    		   {
    			   sampler->delay_x_resample = 1;
    			   sampler->intermediate_w = sampler->resample_src_w;
    		   }
    		   else
    		   {
    			   sampler->delay_x_resample = 0;
    			   sampler->intermediate_w = sampler->resample_dst_w;
    		   }
    #if TEST
    		   printf("delaying: %d\n", sampler->delay_x_resample);
    #endif
    	   }
    	   if (sampler->delay_x_resample)
    	   {
    		   sampler->tmp_buf = (float*)malloc(sampler->intermediate_w * sizeof(float));
    		   if (!sampler->tmp_buf)
    		   {
    			   destroy_lanc_sampler(&sampler);
    			   return NULL;
    		   }
    	   }
    	   return sampler;
    }
    

    插值x方向,读入一行,输出一行:

    static void lanc_sample_x(float* src, float* dst, lanc_sampler* sampler, lanc_coef_node_vec* xcnv)
    {
    	int i, j;
    	float total;
    	lanc_coef* p;
    	if (!xcnv)return;
    	for (i = 0; i < sampler->resample_dst_w; i++)
    	{ 
    		total = 0.0;
    		p = xcnv->nds[i].p;
    		for (j = 0; j < xcnv->nds[i].n; j++)
    		{ 
    			total += src[p[j].pos] * p[j].weight;
    		}
    		*dst++ = total;
    	}
    }
    

    生成一行最终结果:

    static void put_line(real_vec* buf, lanc_sampler* sampler, lanc_coef_node_vec* xcnv)
    {
    	if (!buf || !sampler || !xcnv)
    		return;
    	int i;
    
    	if (sampler->cur_src_y >= sampler->resample_src_h)
    		return;
    
    	/* Does this source line contribute
    	* to any destination line? if not,
    	* exit now.
    	*/
    
    	if (!sampler->src_y_count[resampler_range_check(sampler->cur_src_y, sampler->resample_src_h)])
    	{
    		sampler->cur_src_y++;
    		return ;
    	}
    
    	/* Find an empty slot in the scanline buffer. (FIXME: Perf. is terrible here with extreme scaling ratios.) */
    
    	for (i = 0; i < MAX_SCAN_BUF_SIZE; i++)
    		if (sampler->scan_buf_y[i] == -1)
    			break;
    
    	/* If the buffer is full, exit with an error. */
    
    	if (i >= MAX_SCAN_BUF_SIZE)
    	{
    		return ;
    	}
    
    	sampler->src_y_flag[resampler_range_check(sampler->cur_src_y, sampler->resample_src_h)] = 1;
    	sampler->scan_buf_y[i] = sampler->cur_src_y;
    
    	/* Does this slot have any memory allocated to it? */
    
    	if (!sampler->scan_buf_l[i])
    	{   
    		sampler->scan_buf_l[i] = (float*)calloc(sampler->intermediate_w,sizeof(float));
    		if (!sampler->scan_buf_l[i])
    		{
    			return;
    		}
    	}
    
    	// Resampling on the X axis first?
    	if (sampler->delay_x_resample)
    	{
    		assert(sampler->intermediate_w == sampler->resample_src_w);
    
    		// Y-X resampling order
    		memcpy(sampler->scan_buf_l[i], buf->data, sampler->intermediate_w * sizeof(float));
    	}
    	else
    	{
    		assert(sampler->intermediate_w == sampler->resample_dst_w);
    
    		// X-Y resampling order
    		lanc_sample_x(buf->data,sampler->scan_buf_l[i],sampler, xcnv);
    	}
    
    	sampler->cur_src_y++;
    }
    
    //乘以系数后移动到缓冲区
    static void scale_y_mov(float* Ptmp, const float* Psrc, float weight, int dst_x)
    {
    	int i;
    	for (i = dst_x; i > 0; i--)
    		*Ptmp++ = *Psrc++ * weight;
    }
    
    //一行累加到另一行
    static void scale_y_add(float* Ptmp, const float* Psrc, float weight, int dst_x)
    {
    	for (int i = dst_x; i > 0; i--)
    		(*Ptmp++) += *Psrc++ * weight;
    }
    
    

    截取范围:

    static inline float clamp_sample(float f)
    {
    	if (f < 0.0f)
    		f = 0.0f;
    	else if (f > 1.0f)
    		f = 1.0f;
    	return f;
    }
    
    static void clamp(float* Pdst, int n)
    {
    	while (n > 0)
    	{
    		*Pdst = clamp_sample(*Pdst);
    		++Pdst;
    		n--;
    	}
    }
    

    插值y方向:

    static void lanc_sample_y(lanc_sampler* sampler, lanc_coef_node_vec* xcnv,lanc_coef_node_vec* ycnv)
    {
    	int i, j;
    	float* Psrc;
    	//Contrib_List* Pclist = &m_Pclist_y[cur_dst_y];
    	if (!sampler || !ycnv || !xcnv)return;
    	lanc_coef_node node = ycnv->nds[sampler->cur_dst_y];
    
    	float* Ptmp = sampler->delay_x_resample ? sampler->tmp_buf : sampler->dst_buf;
    	assert(Ptmp);
    
    	/* Process each contributor. */
    
    	for (i = 0; i < node.n; i++)
    	{
    		/* locate the contributor's location in the scan
    		* buffer -- the contributor must always be found!
    		*/
    
    		for (j = 0; j < MAX_SCAN_BUF_SIZE; j++)
    			if (sampler->scan_buf_y[j] == node.p[i].pos)
    				break;
    
    		assert(j < MAX_SCAN_BUF_SIZE);
    
    		Psrc = sampler->scan_buf_l[j];
    
    		if (!i)
    			scale_y_mov(Ptmp, Psrc, node.p[i].weight, sampler->intermediate_w);
    		else
    			scale_y_add(Ptmp, Psrc, node.p[i].weight, sampler->intermediate_w);
    
    		/* If this source line doesn't contribute to any
    		* more destination lines then mark the scanline buffer slot
    		* which holds this source line as free.
    		* (The max. number of slots used depends on the Y
    		* axis sampling factor and the scaled filter width.)
    		*/
    
    		if (--sampler->src_y_count[resampler_range_check(node.p[i].pos, sampler->resample_src_h)] == 0)
    		{
    			sampler->src_y_flag[resampler_range_check(node.p[i].pos, sampler->resample_src_h)] = 0;
    			sampler->scan_buf_y[j] = -1;
    		}
    	}
    
    	/* Now generate the destination line */
    
    	if (sampler->delay_x_resample) // Was X resampling delayed until after Y resampling?
    	{
    		assert(sampler->dst_buf != Ptmp);
    		lanc_sample_x(Ptmp, sampler->dst_buf,sampler, xcnv);
    	}
    	else
    	{
    		assert(sampler->dst_buf == Ptmp);
    	}
    	clamp(sampler->dst_buf, sampler->resample_dst_w);
    }
    

    获取一行结果:

    static float* get_line(lanc_sampler* sampler, lanc_coef_node_vec* xcnv,lanc_coef_node_vec* ycnv)
    {
    	int i;
    
    	/* If all the destination lines have been
    	* generated, then always return NULL.
    	*/
    	if (!sampler || !ycnv)return NULL;
    
    	if (sampler->cur_dst_y == sampler->resample_dst_h)
    		return NULL;
    
    	/* Check to see if all the required
    	* contributors are present, if not,
    	* return NULL.
    	*/
    	lanc_coef* p;
    	int n = ycnv->nds[sampler->cur_dst_y].n;
    	for (i = 0; i < n; i++)
    	{
    		p = ycnv->nds[sampler->cur_dst_y].p;
    		if (!sampler->src_y_flag[resampler_range_check(p[i].pos, sampler->resample_src_h)])
    			return NULL;
    	}
    	lanc_sample_y(sampler, xcnv, ycnv);
    	sampler->cur_dst_y++;
    
    	return sampler->dst_buf;
    }
    
    

    整合所有方法:

    void resize_using_lanczos(uint8_2d* src, uint8_2d* dst)
    {
    	int src_rows, src_cols;
    	int dst_rows, dst_cols;
    	int i, j;
    	int a = 3;
    	if (!src || !dst)return;
    	lanc_coef_node_vec* xcnv = NULL;
    	lanc_coef_node_vec* ycnv = NULL;
    	lanc_sampler*   sampler = NULL;
    	real_2d* s = NULL;
    	real_2d* d = NULL;
    	real_vec* row_buf = NULL;
    	s = uint8_to_real(src->data, src->cols, src->rows);
    	d = create_real_2d(dst->cols, dst->rows);
    	row_buf = create_real_vec(src->cols);
    	if (!s || !d || !row_buf)
    	{   
    		destroy_real_vec(&row_buf);
    		destroy_real_2d(&s);
    		destroy_real_2d(&d);
    		return;
    	}
    	make_lanc_coef(src,dst, &xcnv, &ycnv);
    	normalize_real_2d(s, 255.0f);
    	sampler = create_lanc_sampler(src, dst, xcnv, ycnv);
    	int rows, cols;
    	int row = 0;
    	uint32_t row_bytes;
    	rows = src->rows;
    	cols = src->cols;
    	row_bytes = cols*sizeof(float);
    	for (i = 0; i < rows; i++)
    	{
    		memcpy(row_buf->data, s->arr[i], row_bytes);
    		put_line(row_buf,sampler, xcnv);
    		for (;;)
    		{
    			float* out = get_line(sampler, xcnv, ycnv);
    			if (!out)
    				break;
    			else
    			{
    				for (int x = 0; x < dst->cols; x++)
    				{
    					dst->arr[row][x] = (uint8_t)(out[x] * 255.0f);
    				}
    				row++;
    			}
    		}
    	}
    	destroy_lanc_coef_node_vev(&xcnv);
    	destroy_lanc_coef_node_vev(&ycnv);
    	destroy_lanc_sampler(&sampler);
    	destroy_real_vec(&row_buf);
    	destroy_real_2d(&s);
    	destroy_real_2d(&d);
    }
    

    测试:

    void test_lancos_resize(dai_image* img, float factor)
    {
    	uint8_2d* r = NULL;
    	uint8_2d* g = NULL;
    	uint8_2d* b = NULL;
    
    	uint8_2d* r1 = NULL;
    	uint8_2d* g1 = NULL;
    	uint8_2d* b1 = NULL;
    
    	dai_image* img1 = NULL;
    	if (!img)return;
    	split_img_data(img, &r, &g, &b);
    	int w, h;
    	int w1, h1;
    	w = img->width;
    	h = img->height;
    	w1 = factor*w;
    	h1 = factor*h;
    
    	r1 = create_uint8_2d(h1, w1);
    	g1 = create_uint8_2d(h1, w1);
    	b1 = create_uint8_2d(h1, w1);
    
    	resize_using_lanczos(b, b1);
    	resize_using_lanczos(g, g1);
    	resize_using_lanczos(r, r1);
    
    	merge_img_data(r1, g1, b1, &img1);
    	if (img1)
    	{
    		img1->type = EJPEG;
    		dai_save_image("resize_lancos3.jpg", img1);
    		dai_destroy_image(&img1);
    	}
    	destroy_uint8_2d(&r);
    	destroy_uint8_2d(&g);
    	destroy_uint8_2d(&b);
    	destroy_uint8_2d(&r1);
    	destroy_uint8_2d(&g1);
    	destroy_uint8_2d(&b1);
    }
    
    
    展开全文
  • 图像Lanczos3滤波——C实现

    千次阅读 2019-03-30 09:41:06
    Lanczos 滤波器,因发明者而得名,在信号处理领域,主要用在增加采样率和低通滤波,在图像处理领域用于图像缩放、旋转,其可以在减小锯齿、锐度、振铃效应( aliasing, sharpness, and minimal ringing)取得最好的...

    Lanczos 滤波器,因发明者而得名,在信号处理领域,主要用在增加采样率和低通滤波,在图像处理领域用于图像缩放、旋转,其可以在减小锯齿、锐度、振铃效应( aliasing, sharpness, and minimal ringing)取得最好的平衡。
    滤波函数如下:
    L ( x ) = { sinc ⁡ ( x ) &ThinSpace; sinc ⁡ ( x / a ) if  − a &lt; x &lt; a , 0 otherwise . {\displaystyle L(x)={\begin{cases}\operatorname {sinc} (x)\,\operatorname {sinc} (x/a)&amp;{\text{if}}\ -a&lt;x&lt;a,\\0&amp;{\text{otherwise}}.\end{cases}}} L(x)={sinc(x)sinc(x/a)0if a<x<a,otherwise.
    等价于:
    L ( x ) = { 1 if  x = 0 , a sin ⁡ ( π x ) sin ⁡ ( π x / a ) π 2 x 2 if  − a ≤ x &lt; a  and  x ≠ 0 , 0 otherwise . {\displaystyle L(x)={\begin{cases}1&amp;{\text{if}}\ x=0,\\{\dfrac {a\sin(\pi x)\sin(\pi x/a)}{\pi ^{2}x^{2}}}&amp;{\text{if}}\ -a\leq x&lt;a\ {\text{and}}\ x\neq 0,\\0&amp;{\text{otherwise}}.\end{cases}}} L(x)=1π2x2asin(πx)sin(πx/a)0if x=0,if ax<a and x̸=0,otherwise.
    当a取3时,就是Lanczos3滤波器。
    插值公式:
    S ( x ) = ∑ i = ⌊ x ⌋ − a + 1 ⌊ x ⌋ + a s i L ( x − i ) S(x) = \sum_{i=\lfloor x \rfloor - a + 1}^{\lfloor x \rfloor + a} s_{i} L(x - i) S(x)=i=xa+1x+asiL(xi)
    a是滤波尺寸, ⌊ x ⌋ \lfloor x \rfloor x是向下取整。
    性质:
    1、a>0,Lanczos 核函数对每个点都是连续,并且其倒数也是连续,因此 S ( x ) S(x) S(x)连续,其倒数连续;
    2、Lanczos 核函数在整数点值为零,除了 x = 0,核函数为1,相当于原始x处信号值不变,仅仅对小数位置进行插值;
    二维的情形:
    L ( x , y ) = L ( x ) L ( y ) {\displaystyle L(x,y)=L(x)L(y)} L(x,y)=L(x)L(y)
    算法实现:
    代码实现参考https://github.com/richgel999/imageresampler

    #define  M_PI       3.14159265358979323846
    //sinc函数值
    static double sinc(double x)
    {
    	x = (x * M_PI);
    	if ((x < 0.01f) && (x > -0.01f))
    		return 1.0f + x*x*(-1.0f / 6.0f + x*x*1.0f / 120.0f);
    	return sin(x) / x;
    }
    

    对于0附近的sinc函数值,采用其泰勒展开式计算:
    sinc ⁡ ( x ) = sin ⁡ ( x ) x = ∑ n = 0 ∞ ( − x 2 ) n ( 2 n + 1 ) ! \operatorname {sinc} (x)={\frac {\sin(x)}{x}}=\sum _{n=0}^{\infty }{\frac {\left(-x^{2}\right)^{n}}{(2n+1)!}} sinc(x)=xsin(x)=n=0(2n+1)!(x2)n

    //计算Lanczos3核函数
    static float clip(double t)
    {
    	const float eps = .0000125f;
    	if (fabs(t) < eps)
    		return 0.0f;
    	return (float)t;
    }
    
    static inline float lancos(float t)
    {
    	if (t < 0.0f)
    		t = -t;
    
    	if (t < 3.0f)
    		return clip(sinc(t) * sinc(t / 3.0f));
    	else
    		return (0.0f);
    }
    
    //x方向的插值:
    static inline float lancos3_resample_x(uint8_t** arr,int src_w, int src_h, int y, int x,float xscale)
    {
    	float s = 0;
    	float coef_sum = 0.0f;
    	float coef;
    	float pix;
    	int i;
    
    	int l, r;
    	float c;
    	float hw;
    	//对于缩小的情况hw相当于扩大了领域像素个数,这里如果不作此处理,最终缩小图片效果跟最近领域插值法没多少区别,其效果相当于先低通滤波,再插值
    	if (xscale > 1.0f)
    		hw = 3.0f;
    	else
    		hw = 3.0f / xscale;
    
    	c = (float)x / xscale;
    	l = (int)floor(c - hw);
    	r = (int)ceil(c + hw);
    
    	if (y < 0)y = 0;
    	if (y >= src_h)y = src_h - 1;
    	if (xscale > 1.0f)xscale = 1.0f;
    	for (i = l; i <= r; i++)
    	{   
    		x = i;
    		if (i < 0)x = 0;
    		if (i >= src_w)x = src_w - 1;
    		pix = arr[y][x];
    		coef = lancos((c-i)*xscale);
    		s += pix * coef;
    		coef_sum += coef;
    	}
    	s /= coef_sum; 
    	return s;
    }
    
    
    void img_resize_using_lancos3(uint8_2d* src, uint8_2d* dst)
    {
    	int src_rows, src_cols;
    	int dst_rows, dst_cols;
    	int i, j;
    	if (!src || !dst)return;
    	uint8_t** src_arr;
    	uint8_t** dst_arr;
    	float xratio;
    	float yratio;
    	float pix;
    	int val;
    	int k;
    	float hw;
    	src_arr = src->arr;
    	dst_arr = dst->arr;
    	src_rows = src->rows;
    	src_cols = src->cols;
    	dst_rows = dst->rows;
    	dst_cols = dst->cols;
    
    	xratio = (float)(dst_cols) / (float)src_cols;
    	yratio = (float)(dst_rows) / (float)src_rows;
    
    	float scale = 0.0f;
    	if (yratio > 1.0f)
    	{
    		hw = 3.0;
    		scale = 1.0f;
    	}
    	else
    	{
    		hw = 3.0 / yratio;
    		scale = yratio;
    	}
    	
    	for (i = 0; i < dst_rows; i++)
    	{
    		for (j = 0; j < dst_cols; j++)
    		{
    			int t, b;
    			float c;
    
    			float s = 0;
    			float coef_sum = 0.0f;
    			float coef;
    			float pix;
    
    			c = (float)i / yratio;
    			t = (int)floor(c - hw);
    			b = (int)ceil(c + hw);
                //先对x方向插值,再进行y方向插值。
    			for (k = t; k <= b; k++)
    			{
    				pix = lancos3_resample_x(src_arr, src_cols,src_rows, k, j,xratio);
    				coef = lancos((c - k)*scale);
    				coef_sum += coef;
    				pix *= coef;
    				s += pix;
    			}
    			val = (int)s / coef_sum;
    			if (val < 0)val = 0;
    			if (val > 255)val = 255;
    			dst_arr[i][j] = val;
    		}
    	}
    }
    
    

    测试代码:

    void test_lancos3_resize(dai_image* img, float factor)
    {
    	uint8_2d* r = NULL;
    	uint8_2d* g = NULL;
    	uint8_2d* b = NULL;
    
    	uint8_2d* r1 = NULL;
    	uint8_2d* g1 = NULL;
    	uint8_2d* b1 = NULL;
    
    	dai_image* img1 = NULL;
    	if (!img)return;
    	split_img_data(img, &r, &g, &b);
    	int w, h;
    	int w1, h1;
    	w = img->width;
    	h = img->height;
    	w1 = factor*w;
    	h1 = factor*h;
    
    	r1 = create_uint8_2d(h1, w1);
    	g1 = create_uint8_2d(h1, w1);
    	b1 = create_uint8_2d(h1, w1);
    
    	img_resize_using_lancos3(r, r1);
    	img_resize_using_lancos3(g, g1);
    	img_resize_using_lancos3(b, b1);
    
    	merge_img_data(r1, g1, b1, &img1);
    	if (img1)
    	{
    		img1->type = EJPEG;
    		dai_save_image("resize_lancos3.jpg", img1);
    		dai_destroy_image(&img1);
    	}
    	destroy_uint8_2d(&r);
    	destroy_uint8_2d(&g);
    	destroy_uint8_2d(&b);
    	destroy_uint8_2d(&r1);
    	destroy_uint8_2d(&g1);
    	destroy_uint8_2d(&b1);
    }
    

    效果图:
    原图:
    在这里插入图片描述
    缩小3倍的图:

    在这里插入图片描述
    另外,缩小之后,图像相对变暗,这时,可以采用gama校正改变亮度。

    参考文献:
    1、http://www.realitypixels.com/turk/computergraphics/ResamplingFilters.pdf
    2、https://en.wikipedia.org/wiki/Lanczos_resampling
    3、https://github.com/richgel999/imageresampler

    展开全文
  • lanczos算法用来计算大型稀疏矩阵的最大最小本征值及相应的本征矢量。
  • 因此需要有一种更高效的方法来求解,今天给大家带来的是兰索斯法(Lanczos法)振型求解思路方法。 关于子空间迭代法可以观看下列文章: 【JY】振型求解之子空间迭代 【JY|STR】求解器之三维结构振型分析 目前科学计算...

    一、写在文前

    【前言】子空间迭代法可同时求解几个极端特征值和相应的特征向量,但它有收敛较慢,运算量较大,积累误差的缺点;随后,人们对其作了进一步的研究,出现了预处理子空间迭代法,这种方法的运算量较之子空间迭代法有所减小,但还是比较大,另外,当要求的特征值与不要求的特征值之间的间隔较大时,预处理子空间迭代法收敛也比较慢。因此需要有一种更高效的方法来求解,今天给大家带来的是兰索斯法(Lanczos法)振型求解思路方法。

    关于子空间迭代法可以观看下列文章:

    【JY】振型求解之子空间迭代

    【JY|STR】求解器之三维结构振型分析

        目前科学计算、理论研究、科学实验并列为当今世界科学活动的三种主要方式。许多科学和工程领域如果没有科学计算不可能有一流的研究成果。而在工程领域,矩阵计算具有举足轻重的地位。

    矩阵计算的基本问题有三个:

    • 其一,求解线性方程组问题;

    • 其二,线性最小二乘问题;

    • 其三,矩阵特征值问题。

        在土木工程上,矩阵的特征值具有广泛的应用,如建筑结构中的固有频率分析问题以及钢结构的稳定性,建筑结构的振动问题,大型桥梁的颤振问题等等,都与矩阵的特征值密切相关。一个复杂结构的动力分析和稳定性分析可转化为大型稀疏矩阵的特征值/特征向量问题。

        求解结构振型和频率实质上是对以下式子的特征值和特征周期进行求解。但是,根据伽罗瓦理论,对于次数大于四的多项式求根不存在一种通用的方法。换句话说,对于大型矩阵并无法直接求解得到通法解析解,从而得到特征值和特征向量。于是,人们通过寻求其他的求解途径,提出了许多很好的思想方法和具体算法,促进了该学科的不断发展。

        本期介绍一种针对稀疏矩阵高效的特征值、特征向量计算方法——Lanczos迭代法。Lanczos方法存储量个大,计算速度较快,而且适用于求解矩阵的极端特征值。由于在实际问题中,矩阵的阶数虽然很大,但往往只是需要求其依模最大或最小的特征值,因此Lanczos方法得以广泛的应用,特别是适合求大型稀疏对称矩阵的部分特征值。

        对于Lanczos方法的本质:将实对称矩阵转换成为三对角矩阵(稀疏矩阵),从而便于计算机储存和后续的计算。(这样就把一个求实对称矩阵的特征值问题转化为求一个三对角矩阵的特征值问题。)

        有兴趣的小伙伴也可以了解下另外两种方法将一个实对称矩阵转化为三对角矩阵的方法Givens法及Householdes。

    、方程原理

    结构特征矩阵为:

        假定该结构特征矩阵是大型、稀疏的实对称矩阵,求其几个最大或最小的特征值的问题,可以用Lanczos算法解决,它的原理是先产生一个三对角阵T,于是问题转化为求一个三对角矩阵的特征值,变得相对简便。

        因此,我们的目的是将矩阵A转化为T,并且只要求得T的特征值和特征向量,则可以通过一定关系得到原结构特征矩阵的特征值和特征向量,也即可以得到结构振型和频率。

        首先对于一般结构来说,均需要进行初始向量的预设进行迭代,但是大部分振型都难以提前预估,我们可以听下威尔金森(Wilkinson)的建议预先不必猜测,而统一初始的向量为:

        然后进行预设需求的振型数量为 i≤n (结构特征矩阵维度)进行求解,以下流程图对于自己编写Lanczos方法中的T矩阵集成有较好的理解帮助。

        通过上述的迭代方法,并将求得的α和β带入矩阵T中,即可将结构特征矩阵是大型、稀疏的实对称矩阵化为一个三对角阵T。

        那么结构的频率和振型与该三对角矩阵T的关系是什么呢?

        接下来讨论下T矩阵和结构的频率和振型的关系,继上面的公式推导得到三对角矩阵T,由Gram-Schmidt正交化条件得到:

        继上述推导公式可得到下式:

        将Gram-Schmidt正交化条件带入上式中可得到:

        考虑列向量间的正交性,并将得到的上述公式往会带入可得到:

        得到这个标注了三星的重要级公式!其中标记红框部分是左乘部分。若将红框这么一画,就变成了:

        再通过上述这样一个化解,原结构和变换后的T的频率和振型的关系一目了然。若令原结构的特征向量(振型)为:

    则通过上述可知:

        也就是求解T矩阵的完全解为,则该完全解的特征向量与原结构的特征向量(振型)的关系式之间相差一个系数矩阵X,而特征值(频率)是原结构的倒数。因此取Lanczos向量的个数i等于系统的阶数n,则Lanczos变换法给出了原特征值问题的精确解。

        一般情况下我们取 i ≤ n 。由于组装结构特征矩阵时候,刚度求逆变换了一次,相当于进行了一次逆迭代,因此Lanczos变换法给出了系统低阶特征值的很好的近似解。值得注意的是,为了保证求解精度,在迭代过程中可以用Gram-Schmidt对迭代向量进行重新正交化,并采用移轴法可提高其效率。

        至此,该问题变为求一个i阶三对角阵T的特征值和特征向量。

        对于求该降阶的特征值、特征向量,有许多方法如:比较高效的QR分解法等等,这里就不再多赘述,小伙伴们可以自行编程尝试。

    三、分析对比

        建立一个无楼板的简易框架结构,体验一把Abaqus的兰索斯算法,并对比下自编的Lanczos。模型如下:

    求解前4阶模态

        利用Abaqus 直接提取该结构的质量矩阵、刚度矩阵、阻尼矩阵(非必要),一个简单的结构密密麻麻的数,感兴趣的小伙伴可以试一试:

        并将这三个矩阵放入自行编写的动力计算器(由于该动力计算器尚未完善,所以暂不公开使用,小伙伴们按这个思路也可以自己编写得到属于自己的求解器。)

    结果对比:

        可以看得到结果非常接近,误差来源可能来源于计算机计算的舍入误差。而且计算速度非常高效,因此针对稀疏矩阵高效的特征值、特征向量的计算方法可以采用可存储量个大,计算速度较快的Lanczos迭代法。

    往期精彩

    #性能分析

    【JY】基于性能的抗震设计浅析(一)

    【JY】基于性能的抗震设计浅析(二)

    【JY】浅析消能附加阻尼比

    【JY】近断层结构设计策略分析与讨论

    【JY】浅析各动力求解算法及其算法数值阻尼(人工阻尼)

    理念

    【JY|体系】结构概念设计之(结构体系概念)

    【JY|理念】结构概念设计之(设计理念进展)

    【JY】有限单元分析的常见问题及单元选择

    【JY】结构动力学之显隐式

    【JY】浅谈结构设计

    【JY】浅谈混凝土损伤模型及Abaqus中CDP的应用

    #概念机理

    【JY】基于Ramberg-Osgood本构模型的双线性计算分析

    【JY】结构动力学初步-单质点结构的瞬态动力学分析

    【JY】从一根悬臂梁说起

    【JY】反应谱的详解与介绍

    【JY】结构瑞利阻尼与经济订货模型

    【JY】主成分分析与振型分解

    【JY】浅谈结构多点激励之概念机理(上)

    【JY】浅谈结构多点激励之分析方法(下)

    【JY】板壳单元的分析详解

    【JY】橡胶支座的简述和其力学性能计算

    【JY】振型求解之子空间迭代

    【JY】橡胶支座精细化模拟与有限元分析注意要点

    #软件讨论

    【JY】复合材料分析利器—内聚力单元

    【JY】SDOF计算教学软件开发应用分享

    【JY】Abaqus案例—天然橡胶隔震支座竖(轴)向力学性能

    【JY】Abaqus6.14-4如何关联fortran?

    【JY】如何利用python来编写GUI?

    【JY】如何解决MATLAB GUI编程软件移植运行问题?

    【JY】浅谈结构分析与设计软件

    【JY|STR】求解器之三维结构振型分析

    【JY】SignalData软件开发应用分享

    【JY】基于Matlab的双线性滞回代码编写教程

    #其他

    【JY】位移角还是有害位移角?

    【JY】如何利用python来编写GUI?

    【JY】今日科普之BIM

    ~关注未来更精彩~

    展开全文
  • 360Lib:Lanczos插值

    千次阅读 2017-10-11 17:16:20
    在格式转换中使用的是Lanczos插值滤波器,色度分量使用的是Lanczos-2,亮度分量使用的是Lanczos-3。 Lanczos插值滤波器原理见: http://blog.csdn.net/lin453701006/article/details/78205170 360Lib Lanczos插值...
  • boost::gil::scale_lanczos用法的测试程序实现功能C++实现代码 实现功能 boost::gil::scale_lanczos用法的测试程序 C++实现代码 #include <boost/gil/image.hpp> #include <boost/gil/image_processing/...
  • a2,...an)TA=T^{-1}diag(a_1,a_2,...a_n)TA=T−1diag(a1​,a2​,...an​)T 对于无法进行特征分解的矩阵,以及非方阵,我们常用的处理手段是将矩阵转化为方阵,然后使用双对角化,或者三对角化进行简化,Lanczos对方...
  • Lanczos插值滤波器

    2019-10-18 00:03:13
    Lanczos插值滤波器 https://blog.csdn.net/lin453701006/article/details/78205170 参考: https://en.wikipedia.org/wiki/Lanczos_resampling http://blog.csdn.net/trent1985/article/details/45150677 Lanczos...
  • 多相位图像插值算法(Lanczos、sinc)

    千次阅读 2019-10-18 00:01:08
    多相位图像插值算法(Lanczos、sinc) https://blog.csdn.net/swimmingfish2004/article/details/7312361 Lanczos Algorithm Analyse 在公司时候研究过的Lanczos图像缩略算法,今天整理出来给大家分享,分析的是...
  • Lanczos 重采样

    2019-10-13 00:08:00
    Lanczos 重采样 https://blog.csdn.net/u010468553/article/details/79359488 版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。 本文链接:...
  • 针对光声成像的病态性和较大的系统矩阵会导致重建速度慢的问题,提出了一种基于Lanczos双对角化的快速指数滤波重建方法,并通过数值仿真证实了该方法的有效性。仿真结果表明,所提方法在保证重建图像高质量的同时极...
  • 第三十二篇 Lanczos转化到三对角形式 在之前的篇章里,有许多求解线性方程的迭代方法,如最陡下降法,可以通过向量乘法和各种简单的向量运算,简化为一个单个矩阵的循环。将矩阵化为三对角形式的Lanczos方法,保留其...
  • 图像处理之Lanczos采样放缩算法 一:什么是Lanczos采样 参见这里:http://en.wikipedia.org/wiki/Lanczos_resampling 二:大致算法流程 三:算法运行结果 1.向下采样, 生成缩略图, 左边为原图,右边为缩略图 ...
  • 提出利用多重多级子结构技术与Lanczos方法求解超大型复杂结构动力特性的子结构算法。该算法利用子结构周游树技术,分别对每个子结构进行Lanczos迭代,通过累加各个子结构的正交化系数组成全局三对角矩阵,最后求解...
  • 采用非对称Lanczos 算法研究线性分数阶系统的模型降阶问题, 提出一种保持系统传递函数一定数量的分数阶矩的模型降阶方法. 根据Caputo 导数的运算法则给出线性分数阶系统的分数阶矩的计算方法; 利用非对称Lanczos ...
  • 文章利用Lanczos迭代法求解复杂机械结构振动模态形式,并运用Lanczos迭代法建立了振动筛的有限元动力学模型,进而研究了固有特性,并求解了振动筛的振动模态,为振动筛结构动态特性求解方法提供了理论依据。...
  • 根据青藏高原大气热源季内振荡分析的实际需要,设计了适于提取准双周(10~20 d)、准一月(20~40 d)振荡的Lanczos滤波器(L f.)。通过与Butterworth滤波器(B.f.)滤波效果的定量分析,确定了准双周、准一月L f.窗宽...
  • Lanczos法的目的:将实对称矩阵转换成为三对角矩阵(稀疏矩阵),从而便于计算机储存和后续的计算。 在三对角矩阵矩阵上,采用QR分解,得到矩阵的特征值。 # -*- coding: utf-8 -*- """ Created ...
  • 把求解大规模对称矩阵端部特征对问题的基本方法-Lanczos算法应用于核主成分分析模型的求解,设计了大样本核主成分分析模型求解的实用算法。在clapack和nu-TRLan两个软件包的基础上,开发了大样本核主成分分析模型...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,036
精华内容 1,614
关键字:

lanczos