Three-Dimensional Krawtchouk Descriptors for Protein Local Surface Shape Comparison
Atilla Sit Daisuke Kihara



























label1]Atilla Sit
[label1]Department
of Mathematics and Statistics, Eastern Kentucky University, Richmond, KY, 40475 USA
atilla.sit@eku.edu


label2,label3,label4]Daisuke Kiharacor1
[label2]Department of Biological Sciences, Purdue University, West Lafayette, IN, 47907 USA
[label3]Department of Computer Science, Purdue University, West Lafayette, IN, 47907 USA
[label4]Department of Pediatrics,  Cincinnati Children's Hospital Medical Center, University of Cincinnati, Cincinnati, OH, 45229 USA
[cor1]Corresponding author.
dkihara@purdue.edu







Direct comparison of three-dimensional (3D) objects is computationally expensive due to the need for translation, rotation, and scaling of the objects to evaluate their similarity. In applications of 3D object comparison, often identifying specific local regions of objects is of particular interest. We have recently developed a set of 2D moment invariants based on discrete orthogonal Krawtchouk polynomials for comparison of local image patches. In this work, we extend them to 3D and construct 3D Krawtchouk descriptors (3DKD) that are invariant under translation, rotation, and scaling. The new descriptors have the ability to extract local features of a 3D surface from any region-of-interest. This property enables comparison of two arbitrary local surface regions from different 3D objects. We present the new formulation of 3DKD and apply it to the local shape comparison of protein surfaces in order to predict ligand molecules that bind to query proteins.





3D image retrieval local image comparison region of interest discrete orthogonal functions Krawtchouk polynomials weighted Krawtchouk polynomials 3D Krawtchouk moments protein surface ligand binding site pocket comparison structure-based function prediction















§ INTRODUCTION



Moment-based approaches have become very popular in 2D and 3D image processing due to their compact representation of images. A moment-based approach characterizes a 2D or 3D image by considering its shape as a mathematical function and computes integral of the function multiplied by specific base functions. The approach has been used in many problems including reconstruction, detection, pattern recognition, and compression of images <cit.>. The theory of moment invariants in 2D has been well established since the foundation of algebraic Hu invariants <cit.>. Sadjadi and Hall <cit.> extended the algebraic 2D invariants to 3D and explicitly derived the second order moment invariants, which were later reproduced by Guo <cit.>. Using a group theoretic approach, Lo and Don <cit.> constructed twelve complex moment invariants including both second and third order moments. Galvez and Canton <cit.> defined the 3D moments by evaluating them on the 3D object's surface and extracting global descriptors from normalized surface shapes. An extension of moment invariants to n-dimension can be found in Mamistvalov's work <cit.>, in which the zeroth and second order moment invariants of n-dimensional regular solids were established. Other examples of 3D moments that are invariant to rotation and blur were provided by Flusser et al. <cit.>.

We have recently developed a set of 2D local moment invariants based on the discrete Krawtchouk polynomials and successfully applied them to the comparison of local image patches <cit.>. Krawtchouk polynomials were introduced by Mikhail Krawtchouk <cit.> and used for the first time in image analysis by Yap et al. <cit.>. They have been employed in image processing and pattern recognition fields including image reconstruction <cit.>, image watermarking <cit.>, and face recognition <cit.>. In our previous 2D work <cit.>, while constructing a set of local descriptors that are rotation, position, and size independent, we have also preserved their ability to extract features from any local interest region in an image.

In this paper, we extend our previous 2D work to 3D for local comparison of 3D surface shapes. Our new method is based on 3D Krawtchouk polynomials. 3D Krawtchouk moments were earlier defined and used in content-based search and retrieval of 3D objects <cit.>. Despite the compact representation and discriminative powers of these moments, the theory of invariants based on 3D Krawtchouk polynomials has not been well studied. Also, the very critical local retrieval property of the 3D moments has been noticed in <cit.>, but much of the focus is given to their fast computation.
We propose a new approach on this long-standing issue of local image comparison by constructing 3D Krawtchouk descriptors (3DKD) for describing local 3D surfaces. The new formulation has many advantages over many similar moment-based approaches, such as TRS invariants <cit.> and Zernike descriptors <cit.>: 1) Krawtchouk polynomials are defined on a discrete space, so the moments derived from them do not carry any error due to discretization unlike many other moments related to continuous functions <cit.>. 2) These polynomials are orthogonal; each moment brings in a new feature of the image, where minimum redundancy is critical in their discriminative performance. Moreover, they are directly defined in the image coordinate space, and hence their orthogonality property is well retained in the computed moments. 3) They are complete with a finite number of functions (equal to the image size) while many other polynomial spaces have infinitely many members. 4) They have the ability to retrieve local image patches by only changing the resolution of reconstruction and using low order moments. 5) The location of the patch can also be controlled by changing three parameters and hence shifting the region-of-interest along each dimension. 6) We also prove that these moments can be transformed into local descriptors, which are invariant under translation, rotation, and scaling. Therefore, using only a small number of invariant descriptors per image will make it possible to develop an efficient method for quick local image retrieval.

Moment-based approaches, particularly Krawtchouk moments, are very useful for representing biological and medical images as they are pixelized or voxelized data. In medical imaging, such as computerized tomography (CT) scan and magnetic resonance imaging (MRI), objects are observed at different viewpoints and local images need to be extracted and examined. In digital pathology, for instance, pathologists are interested in information about specific structures rather than the whole image <cit.>. Thus, it is necessary to construct moment invariants that do not change by translation, rotation, and scaling and can retrieve local image patches or subimages.

Local shape search methods have many applications also in structural biology, which deals with 3D structures of biomolecules. An important application is the identification of ligand molecules (i.e., small chemical compounds including drug molecules) that bind to local protein surface regions, which is important for predicting biological function of proteins <cit.> and for computational drug design <cit.>. Ligand molecules that bind to a local surface region in a protein can be predicted by finding similar local regions (binding pockets) of known ligand-binding proteins in the protein structure database.
In this work, we applied the developed 3DKD for the protein ligand binding pocket comparison. A ligand binding pocket is represented as a combination of overlapping local surface patches, each of which is characterized by its geometric shape. The shapes of surface patches are compactly represented by 3DKD. The method is benchmarked on a dataset, which contains a total of 463 proteins that bind to one of 11 ligand molecules. The 3DKD-based method shows higher prediction accuracy than a previously developed method, Pocket-Surfer <cit.>, for seven of the 11 ligand types.

This paper is organized as follows. In Section <ref>, we give a brief background of one-dimensional Krawtchouk polynomials. Then, after introducing the 3D weighted Krawtchouk polynomials and their moments in Section <ref>, we present the theory and formulation of our new 3D Krawtchouk descriptors in Section <ref>. In Section <ref>, we provide a detailed scheme for efficient computation of these descriptors. In Section <ref>, we show numerical results from local surface recognition performances of 3DKD using protein structures placed in different orientations. Finally, we discuss the application of 3DKD on the comparison of ligand binding pockets on protein surfaces. We finish the paper with a conclusion and summary of this work in Section <ref>.







§ KRAWTCHOUK POLYNOMIALS



We start with introducing one-dimensional Krawtchouk polynomials, which can also be found in <cit.>. A more general and abstract form of these polynomials was provided as Hahn polynomials in <cit.>.

The nth order Krawtchouk polynomials are defined as

    K_n(x;p,N)     = ∑_i=0^n a_i,n,p,N  x^i  = _2 F_1 (-n,-x;-N;1/p)

where x,n=0,…,N, N>0, p∈ (0,1) and the function _2 F_1 is the hypergeometric function which is defined as:

    _2 F_1 (a,b;c;z) = ∑_i=0^∞(a)_i(b)_i/(c)_iz^i/i!.

The symbol (a)_i in (<ref>) is the Pochhammer symbol given by

    (a)_i = a(a+1)(a+2)… (a+i-1) = Γ(a+i)/Γ(a).


Note that the series in (<ref>) terminates if either a or b is a nonpositive integer. Hence, the polynomial coefficients a_i,n,p,N in (<ref>) can be obtained by simplifying the summation. The first three Krawtchouk polynomials are

    K_0(x;p,N)      = 1, 
    
    	K_1(x;p,N)      = 1 - (1/Np)x, 
    
    	K_2(x;p,N)      = 1 - (2/Np + 1/N(N-1)p^2)x 
           + (1/N(N-1)p^2)x^2.


It is shown in <cit.> that the range of Krawtchouk polynomials expands rapidly with the increase of the order. Besides, these polynomials are not numerically stable for large values of N. Hence, a more stable set of polynomials can be obtained from the classical Krawtchouk polynomials by normalizing with the norm and scaling by the square root of a weight function <cit.>. The weighted Krawtchouk polynomials is then defined by

    K̅_n(x;p,N) = K_n(x;p,N) √(w(x;p,N)/ρ(n;p,N))

where

    w(x;p,N)      = Nx p^x (1-p)^N-x, 
    ρ(n;p,N)      = (-1)^n (1-p/p)^n n!/(-N)_n.

The set of weighted Krawtchouk polynomials

    S̅ = {K̅_n (x;p,N) : n = 0,…,N }

becomes a complete orthonormal set of basis functions on the discrete space { 0,…,N} with the orthonormality condition

    ∑_x=0^NK̅_n(x;p,N) K̅_n'(x;p,N) = δ_nn'.


To compute the weighted Krawtchouk polynomials, the three-term recurrence relation given in <cit.> can be used. Such a recursive computation is shown to be more efficient than computing high order polynomials directly using (<ref>) and (<ref>). However, due to error propagation, computing polynomials recursively may still be numerically unstable for large N as noted by Zhang et al. <cit.>. To achieve numerical stability, we use symmetry and bi-recursive algorithm given in <cit.>.




§ THREE-DIMENSIONAL WEIGHTED KRAWTCHOUK MOMENTS



In this section, we give a brief formulation of 3D weighted Krawtchouk moments, which are also introduced in <cit.>. Note that the functions K̅_n defined by (<ref>) are orthonormal in the one-dimensional discrete set {0,…,N}, but they can be easily extended to three-dimension as follows:

Let

    A = {0,…,N}×{0,…,M}×{0,…,L}

be a discrete field in the 3D space. We define the set of 3D weighted Krawtchouk polynomials on A as

    S̅ =     {K̅_n (x;p_x,N) ·K̅_m (y;p_y,M) ·K̅_l (z;p_z,L): 
         n = 0,…,N,   m = 0,…,M,    l = 0,…,L }.


Note that S̅ is orthonormal on A with the orthonormality condition

    ∑_x=0^N∑_y=0^M   ∑_z=0^L K̅_n (x;p_x,N)  K̅_m (y;p_y,M)  K̅_l (z;p_z,L) 
    ·    K̅_n' (x;p_x,N)  K̅_m' (y;p_y,M)  K̅_l' (z;p_z,L) 
    
    	      =      δ_nn'δ_mm'δ_ll',

which follows immediately from the orthonormality of 1D functions given by (<ref>).

Let f(x,y,z) be a 3D function defined on the grid A in (<ref>). The 3D weighted Krawtchouk moments of order n+m+l of f(x,y,z) are defined by

    Q̅_nml =     ∑_x=0^N∑_y=0^M∑_z=0^L f(x,y,z) 
    ·    K̅_n (x;p_x,N)  K̅_m (y;p_y,M)  K̅_l (z;p_z,L).


Note that by using (<ref>) and solving (<ref>) for f(x,y,z), the 3D function f(x,y,z) can be written in terms of the 3D weighted Krawtchouk polynomials, i.e.,

    f(x,y,z) =    ∑_n=0^N∑_m=0^M∑_l=0^L Q̅_nml
    ·    K̅_n (x;p_x,N)  K̅_m (y;p_y,M)  K̅_l (z;p_z,L).


This means that the object f(x,y,z) can be reconstructed perfectly if all the moments Q̅_nm are used for n=0,…,N, m=0,…,M, l=0,…,L. An approximate reconstruction f̂(x,y,z) of f(x,y,z) can be written as

    f̂(x,y,z) =    ∑_n=0^N̂∑_m=0^M̂∑_l=0^L̂ Q̅_nml
    ·    K̅_n (x;p_x,N)  K̅_m (y;p_y,M)  K̅_l (z;p_z,L).

where 0≤N̂≤ N, 0≤M̂≤ M, 0≤L̂≤ L. 






































Fig. <ref> presents some reconstructions of 3D binary images using 3D weighted Krawtchouk polynomials for N̂, M̂, L̂ values of 5, 10, 25, and 50, and different (p_x,p_y,p_z) triplets. The 3D polygonal models for the horse and the mug image are downloaded from Princeton Shape Benchmark <cit.> and voxelized using the algorithm in <cit.>. The 3D weighted Krawtchouk moments from the original image are first computed using (<ref>), and then these moments are used in (<ref>) for reconstructing the image.The center of a local region corresponding to (p_x,p_y,p_z) is at (x_c,y_c,z_c) where x_c=N p_x, y_c=M p_y, and z_c=L p_z. These points are (97,126,167) near the horse's mouth and (169,145,100) near the mug's handle. Since N=M=L=200 in this example, the (p_x,p_y,p_z) triplets at these centers will correspond to (0.485,0.630,0.835) and (0.845,0.725,0.500), respectively.

As can be seen from left to right in Fig. <ref>, the reconstructions start at a local region corresponding to (Np_x, Mp_y, Lp_z) and expand as larger values of N̂, M̂, and L̂ are used.
Theoretically, using N̂=N=200, M̂=M=200, and L̂=L=200, the original image will be fully reconstructed regardless of the choice of (p_x,p_y,p_z). Using smaller numbers for N̂, M̂, and L̂, the reconstructed surfaces contain only local information which may actually be more useful for local comparison of 3D images. The parameters p_x, p_y, and p_z play a vital role here to determine the center of local region-of-interest. 

In the third and fourth row of Fig. <ref>, we show the voxelized surface of a protein, nucleosome recognition module of imitation SWI ATPase from fruit fly (Drosophila melanogaster). The atomic structure of this protein is downloaded from the Protein Data Bank (PDB) <cit.> (PDB ID: 1OFC) and then is voxelized using 3D-Surfer <cit.>. The mean radius[The mean radius is calculated as the average of all distances between the center of mass of the protein and each amino acid.] of this protein is about 27.8 Å. In the grid shown, 1 unit corresponds to 0.6 Å. For this example, two (p_x,p_y,p_z) triplets are selected: (0.66,0.67,0.82) and (0.21,0.44,0.31) in the third and the fourth row, respectively. So, the reconstruction centers will be (x_c,y_c,z_c)=(132,134,164) and (42,88,62), respectively, both chosen from salient parts on the surface of the protein. Again, using smaller numbers for N̂, M̂, and L̂, the reconstructed surfaces will contain only local information. Local retrieval of structures may reveal important information about the function of a protein, and this may be used to locally compare protein structures in a large database and quickly identify their ligand-binding sites. This may be very useful for identifying biological functions of proteins and further computational drug design for target proteins. In this paper, we will employ p_x, p_y, and p_z parameters for detecting the local region-of-interest by changing them between 0 to 1.











§ 3D KRAWTCHOUK DESCRIPTORS



In this section, we introduce a new set of invariants, called 3D Krawtchouk descriptors. We show that these invariants are not only rotation, size, and position independent, but also contain discriminative local features from any region-of-interest in a 3D image. Such invariants in 2D have been introduced in our previous work <cit.>. In this work, we extend them to 3D to locally compare 3D images.

Let f(x,y,z) be a function representing a 3D image defined on an orthogonal grid A given in (<ref>) and define the 3D weight function corresponding to triplets 𝐩=(p_x,p_y,p_z) and 𝒩=(N,M,L) by

    W(x,y,z;𝐩,𝒩) = w(x;p_x,N) w(y;p_y,M) w(z;p_z,L).

where x=0,…,N, y=0,…,M, and z=0,…,L. Similarly, one can define the 3D norms corresponding to triplets 𝐩 and 𝒩 by

    Ω (n,m,l;𝐩,𝒩)=ρ(n;p_x,N) ρ(m;p_y,M) ρ(l;p_z,L)

for n=0,…,N, m=0,…,M, and l=0,…,L. Using (<ref>), the 3D weighted Krawtchouk moments Q̅_nml in (<ref>) become

    Q̅_nml =      [Ω (n,m,l;𝐩,𝒩)]^-1/2·∑_x=0^N∑_y=0^M∑_z=0^L f̃(x,y,z)
       ·  K_n (x;p_x,N)   K_m (y;p_y,M)   K_l (z;p_z,L)

where

    f̃(x,y,z) = [W(x,y,z;𝐩,𝒩)]^1/2  f(x,y,z).
 
Now, substituting K_n, K_m, and K_l in (<ref>) by their definitions from (<ref>), reordering summations, and grouping terms, we obtain

    Q̅_nml =    [Ω (n,m,l;𝐩,𝒩)]^-1/2
    ·   ∑_i=0^n∑_j=0^m∑_k=0^l  a_i,n,p_x,N  a_j,m,p_y,M  a_k,l,p_z,L M̃_ijk,

where

    M̃_ijk = ∑_x=0^N∑_y=0^M∑_z=0^L x^i y^j z^k f̃(x,y,z)

are the geometric moments of the auxiliary function in (<ref>).

Notice that the geometric moments M̃_ijk and hence the weighted Krawtchouk moments Q̅_nml are not invariant under translation, rotation, and scaling. The translation invariant central moments of f̃(x,y,z) can be defined as

    μ̃_nml =  ∑_x=0^N∑_y=0^M∑_z=0^L    f̃(x,y,z) (x-x̃)^n (y-ỹ)^m (z-z̃)^l ,

where x̃=M̃_100/M̃_000, ỹ=M̃_010/M̃_000, and z̃=M̃_001/M̃_000 are the centroids of the auxiliary image f̃(x,y,z). 

If μ̃_nml are the central moments, then  we can define geometric moments of f̃(x,y,z), which are invariant under translation and scaling as follows:

    η̃_nml = μ̃_nml/(M̃_000)^n+m+l/3+1.


Obtaining rotation invariant geometric moments of f̃(x,y,z) is, however, not as straightforward. To achieve rotational invariance, we rotate the auxiliary image f̃(x,y,z) so that its principal axes lie in the x, y, z-directions, respectively. When the auxiliary image f̃(x,y,z) is centered at the origin, the principal axes of f̃(x,y,z) can be defined as the eigenvectors of the inertia matrix 


    Ĩ = [ 
    [  Ĩ_xx -Ĩ_xy -Ĩ_xz; -Ĩ_yx  Ĩ_yy -Ĩ_yz; -Ĩ_zx -Ĩ_zy  Ĩ_zz ]]

where


    Ĩ_xx    = μ̃_020+μ̃_002
    Ĩ_yy    = μ̃_200+μ̃_002
    Ĩ_zz    = μ̃_200+μ̃_020

and


    Ĩ_xy    = Ĩ_yx = μ̃_110
    Ĩ_xz    = Ĩ_zx = μ̃_101
    Ĩ_yz    = Ĩ_zy = μ̃_011.


Here, Ĩ is a symmetric matrix with real eigenvalues {λ̃_1, λ̃_2, λ̃_3} and orthogonal eigenvectors {ũ_1, ũ_2, ũ_3} such that

    Ĩ ũ_i     = λ̃_i ũ_i     i=1,2,3.

The eigenvectors {ũ_1, ũ_2, ũ_3} define the rows of the rotation matrix R̃ that aligns the principal axes with the standard xyz coordinate system. However, for each eigenvalue λ̃_i, both ũ_i and -ũ_i are eigenvectors, so they define eight different rotation matrices {±ũ_1, ±ũ_2, ±ũ_3} specifying the same principal axes <cit.>.

Although it is possible to reduce this ambiguity to four combinations by only keeping right-handed coordinate systems <cit.>, a heuristic approach is still needed to obtain a unique standard rotation matrix. In our work, we will use the invariants to locally compare 3D image surfaces. For this reason, we first locate points on the surface of the 3D image. Then, for each local patch around these points, a direction vector can be specified using the vertex normal at the surface point, pointing away from the surface. Among the eight rotations, the unique rotation matrix R̃ can be chosen, for example, as the one rotating the vertex normals to the octant in which x, y, and z coordinates are all nonpositive. Once R̃ is determined, we can define geometric moments of f̃(x,y,z), which are invariant under rotation, translation, and scaling by

    ν̃_ijk =      (M̃_000)^-i+j+k/3-1∑_x=0^N∑_y=0^M∑_z=0^L f̃(x,y,z) 
     · (ϕ̃_1   (x,y,z)) ^i  ·  ( ϕ̃_2(x,y,z) )^j ·  ( ϕ̃_3(x,y,z) )^k,

where ϕ̃_1(x,y,z) = R̃_11(x-x̃) + R̃_12(y-ỹ) + R̃_13(z-z̃), ϕ̃_2(x,y,z) = R̃_21(x-x̃) + R̃_22(y-ỹ) + R̃_23(z-z̃), and ϕ̃_3(x,y,z)=R̃_31(x-x̃) + R̃_32(y-ỹ) + R̃_33(z-z̃).



Fig. <ref> shows 2D views of the square root of the 3D weight function W(x,y,z;𝐩,𝒩) in (<ref>) for 𝒩=(200,200,200) and five different 𝐩=(p_x,p_y,p_z) triplets. To obtain the 2D views, we flatten the function by summing slices along each of the three dimensions.  Note that the coverage of each function differs as the parameters 𝐩=(p_x,p_y,p_z) change. The coverage is the largest for 𝐩=(0.5,0.5,0.5) and becomes smaller when 𝐩 approaches faces of the grid. Different weight functions result in loss of translation and scale invariance. One way to overcome this problem is to determine a unique suitable 𝐩 triplet, say 𝐩^∗=(p_x^∗,p_y^∗,p_z^∗), and use the corresponding weight W(x,y,z;𝐩^∗,𝒩) for every local region-of-interest by shifting the graph of W to that location, in other words, use the translated weight W^∗(x,y,z;𝐩^∗,𝒩)=W(x^∗,y^∗,z^∗;𝐩,𝒩) with x^∗=x-Np_x^∗+Np_x, y^∗=y-Mp_y^∗+Mp_y, and z^∗=z-Lp_z^∗+Lp_z. Whenever (x^∗, y^∗,z^∗) is situated outside the grid, we set W^∗(x,y,z;𝐩^∗,𝒩)=0. From now on in this paper, we will set 𝐩^∗=(0.5,0.5,0.5) due to the largest coverage and round shape of the corresponding weight function, which is also critical for rotational invariance. In order to preserve the round shape of the weight function, we will always use a cubical grid, i.e., N=M=L or set 𝒩=(N,N,N). Hence, the center of the grid will be at C=(N/2,N/2,N/2). In order to shift the centroid of the auxiliary image f̃(x,y,z) to the grid center C, ν̃_ijk in (<ref>) is modified to


    λ̃_ijk =      (M̃_000)^-1∑_x=0^N∑_y=0^N∑_z=0^N f̃(x,y,z) 
    ·    (  ϕ̃_1(x,y,z)  / (M̃_000)^1/3 +  N/2   ) ^i 
    ·    (  ϕ̃_2(x,y,z)  / (M̃_000)^1/3 +  N/2   ) ^j 
    ·    (  ϕ̃_3(x,y,z)   / (M̃_000)^1/3 +  N/2   ) ^k.

Using the binomial expansion, λ̃_ijk can be written as




    λ̃_ijk    =   (M̃_000)^-1∑_x=0^N∑_y=0^N∑_z=0^Nf̃(x,y,z) 
    ·   ∑_r=0^iir(ϕ̃_1(x,y,z)/(M̃_000)^1/3)^r (N/2)^i-r
    ·   ∑_s=0^jjs(ϕ̃_2(x,y,z)/(M̃_000)^1/3)^s (N/2)^j-s
    ·   ∑_t=0^kkt(ϕ̃_3(x,y,z)/(M̃_000)^1/3)^t (N/2)^k-t.

which can then be rearranged to

    λ̃_ijk    =  ∑_r=0^i∑_s=0^j∑_t=0^kirjskt  (N/2)^i+j+k-r-s-t
       ·   (M̃_000)^-r+s+t/3-1∑_x=0^N∑_y=0^N∑_z=0^Nf̃(x,y,z) 
       ·  (ϕ̃_1(x,y,z))^r · (ϕ̃_2(x,y,z))^s · (ϕ̃_3(x,y,z))^t.

Thus, 

    λ̃_ijk  =     ∑_r=0^i∑_s=0^j∑_t=0^kirjskt  (N/2)^i+j+k-r-s-t ν̃_rst

is a linear combination of invariants ν̃_rst where

    ν̃_rst =      (M̃_000)^-r+s+t/3-1∑_x=0^N∑_y=0^N∑_z=0^N f̃(x,y,z) 
     · (ϕ̃_1   (x,y,z)) ^r  ·  ( ϕ̃_2(x,y,z) )^r ·  ( ϕ̃_3(x,y,z) )^t,

for r=0,…,i, s=0,…,j and t=0,…,k. Therefore, these new geometric moments are rotation, translation, and scale invariant, and yet centered at the point (N/2,N/2,N/2). If we set p_x=0.5, p_y=0.5, and p_z=0.5 in (<ref>) and replace M̃_ijk by their invariant counterparts λ̃_ijk from (<ref>), we obtain a new set of moments which are invariant under rotation, translation, and scaling, i.e.,

    Q̃_nml =      [Ω (n,m,l;𝐩^∗,𝒩)]^-1/2
    ·∑_i=0^n   ∑_j=0^m∑_k=0^l  a_i,n,0.5,N  a_j,m,0.5,N  a_k,l,0.5,N λ̃_ijk.

where 𝐩^∗=(0.5, 0.5, 0.5) and 𝒩=(N,N,N).
This new set of moments will be called 3D Krawtchouk descriptors and referred as 3DKD in the rest of the paper. Note that 3DKD still depend on the number N, so it is important to use the same grid while performing the local comparison of 3D images.







§ COMPUTATION OF DESCRIPTORS



The descriptors defined in (<ref>) requires precomputation of λ̃_ijk in (<ref>) which is a linear combination of geometric moments ν̃_rst given by (<ref>). Notice that computation of ν̃_rst requires exponentiation of three 3D functions ϕ̃_1, ϕ̃_2, and ϕ̃_3, and then element-wise multiplication of four 3D functions, and finally summation over three variables. Thus, direct computation in (<ref>) of ν̃_rst can be quite time-consuming, especially when the grid size 𝒩=(N,N,N) is large. To save computational time, we will first separate the z variable from the x and y variables by rewriting (ϕ̃_1(x,y,z))^r, (ϕ̃_2(x,y,z))^s, and (ϕ̃_3(x,y,z))^t using the binomial expansion as follows:

    (ϕ̃_1(x,y,z))^r =     ∑_ε_1=0^r rε_1 (Ã_1(x,y))^r-ε_1 (D̃_1(z))^ε_1,


    (ϕ̃_2(x,y,z))^s =     ∑_ε_2=0^s sε_2 (Ã_2(x,y))^s-ε_2 (D̃_2(z))^ε_2,


    (ϕ̃_3(x,y,z))^t =     ∑_ε_3=0^t tε_3 (Ã_3(x,y))^t-ε_3 (D̃_3(z))^ε_3,

where Ã_τ(x,y) = R̃_τ 1(x-x̃) + R̃_τ 2(y-ỹ) and D̃_τ(z) = R̃_τ 3(z-z̃) with τ=1,2,3.

Hence, given a voxelized 3D surface function f(x,y,z), a triplet 𝐩=(p_x,p_y,p_z) corresponding to a point on the surface, and the surface normal at that point, an efficient computation of the geometric moments ν̃_rst in (<ref>) can be performed in the following steps.


	
  * Compute the auxiliary image f̃(x,y,z) using (<ref>).
    
    
  * Compute the function
    	
    T̃(x,y;ε_1,ε_2,ε_3) 
    
                    =     ∑_z=0^N f̃(x,y,z)   (D̃_1(z))^ε_1 (D̃_2(z))^ε_2 (D̃_3(z))^ε_3

       	for 0≤ε_1+ε_2+ε_3 ≤ 5. Note that the z variable is eliminated in this step and the rest of the computations will be carried out in the x and y variables only.
        
	
  * Compute 
    	
    T̃_3(x,y;ε_1,ε_2,t) 
    
                    =     ∑_ε_3=0^t tε_3  (Ã_3(x,y))^t-ε_3 T̃(x,y;ε_1,ε_2,ε_3)

       	for 0≤ε_1+ε_2+t ≤ 5.
        
	
  * Compute 
    	
    T̃_2(x,y;ε_1,s,t) 
    
                    =     ∑_ε_2=0^s sε_2  (Ã_2(x,y))^s-ε_2 T̃_3(x,y;ε_1,ε_2,t)

       	for 0≤ε_1+s+t ≤ 5.
        
	
  * Compute 
    	
    T̃_1(x,y;r,s,t) 
    
                    =     ∑_ε_1=0^r rε_1  (Ã_1(x,y))^r-ε_1 T̃_2(x,y;ε_1,s,t)

       	for 0≤ r+s+t ≤ 5.

	
  * Compute 
    	
    ν̃_rst =     (M̃_000)^-r+s+t/3-1 ∑_x=0^N∑_y=0^N T̃_1(x,y;r,s,t)

       	for 0≤ r+s+t ≤ 5.
        


Finally, we perform another couple of steps to compute 3DKD of the order up to 5.


	
	
  * Compute λ̃_ijk in (<ref>) for 0≤ i+j+k≤ 5 using ν̃_rst from step 6.
    
  * Compute Q̃_nml in (<ref>) for 0≤ n+m+l≤ 5 using λ̃_ijk from step 7.


Separating the z variable as described above makes the computations very efficient. The computational performance of the algorithm will be shown at the end of Section <ref>.







§ RESULTS AND DISCUSSION



In this section, we test the local discriminative performance of 3DKD. We use three feature vectors of descriptors
 
    K_3       =  {Q̃_nml: 0≤ n+m+l≤ 3}
    
         K_4       =  {Q̃_nml: 0≤ n+m+l≤ 4}
    
         K_5       =  {Q̃_nml: 0≤ n+m+l≤ 5}

namely, the descriptors of order up to 3, 4, and 5, that are computed using the algorithm summarized in Section <ref>. The number of elements in K_3, K_4, and K_5 is 20, 35, and 56, respectively. The seven descriptors (Q̃_000, Q̃_100, Q̃_010, Q̃_001, Q̃_011, Q̃_101, Q̃_110) involved in the normalization process were removed, because they take a constant value irrespective of the 3D patch we are working with. The numbers of descriptors used in local comparison are given in Table <ref>.



As the similarity measure, we use the (squared) Euclidean distance between feature vectors of the same size, namely

    d(^q,^t) =     ∑_i=1^T (_i^q - _i^t)^2

where ^q and ^t are the feature vectors for a query and a target object, respectively, to be compared, and T is dimension of the feature vector.








 §.§ Local Comparison Test I


We first test 3DKD for comparison of local patches on protein surfaces. We have downloaded the PDB file of a protein (PDB ID: 1GCO) from PDB <cit.>, and generated a voxelized surface of the protein using 3D-Surfer <cit.>. Using the GETPOINTS subroutine of LZerD docking suite <cit.>, we have specified 1608 vertex points on the protein surface as well as the normal vectors at these points pointing outside the surface. We have then reduced the number of points to 500 so that the minimal distance between neighboring points is more than 3 Å (See Fig. <ref>, left). Each of these points and the normal vectors are used to represent the center of a patch on the surface. These normal vectors are also used for determining the unique rotation matrix R as described in Section <ref>. The 3DKD corresponding to each patch is computed and stored as the query dataset. Then, the protein structure and the 500 surface points are rotated so that each patch moves to a different location and orientation (see Fig. <ref>). We have used three different rotation matrices S_1, S_2, and S_3 and obtained three target sets. The 3DKD corresponding to each patch in the target sets is computed and stored as Target 1, Target 2, and Target 3. The rotations S_1, S_2, and S_3 are given in (<ref>). S_1 rotates the query protein 90^∘ about the x-axis, S_2 rotates the query protein first 90^∘ about the x-axis, then 45^∘ about the y-axis, and then 30^∘ about the z-axis, and S_3 is a randomly generated rotation matrix.


    S_1   =     [
        	1      0     0 
     
        	0      0     1 
     
        	0     -1     0
        ] 
    
        S_2   =     [√(6)/4    √(6)/4     1/2 
     
        	-√(2)/4      -√(2)/4    √(3)/2 
    √(2)/2     -√(2)/2     0
        ] 
    
        S_3   =     [
        	-0.8256      0.4039     -0.3940 
     
        	-0.2008     -0.8629     -0.4637 
     
        	-0.5273     -0.3037      0.7936
        ]


Each of the 500 feature vectors in the first set is queried and compared with the 500 feature vectors in a target set. The results are ranked using the Euclidean distance. If the same patch as the query ranks top (Top 1), it is labeled as “correctly classified”. Otherwise, we look at the top five results in the ranking (Top 5), or the top ten results (Top 10). We have then calculated the recognition accuracies as 

    r  =

whose denominator is equal to 500.




We have obtained very high recognition accuracies ranging between 93.6% and 97% for the Top 1 case (see Table <ref>). The recognition accuracies stay between 97.6% and 99% for the Top 5 case. For the Top 10 case, it reaches up to between 98% and 100%. The results show that the 3DKD's are successful in local patch comparison for this small problem. For the Top 1 case, K_5 performs best with all three targets. K_4 gives the highest recognition accuracies for all targets when our classification criterion is relaxed to Top 5 or Top 10.







 §.§ Local Comparison Test II


We also test the local performance of 3DKD on a more difficult problem as follows. The surface grid, vertex points and normal vectors of the query protein, namely 1608 query patches, are generated as before. For the three target sets, the target proteins and their voxelized surface grids are also generated as before. However, the surface points and vertex normals of the target sets are redistributed using GETPOINTS so that we have a new set of points and normal vectors different than the ones obtained before. This occurs due to a randomized subroutine of the program. This time, we have obtained 1557, 1557, and 1546 such points in Target 1, Target 2, and Target 3, respectively. For each point and the corresponding surface patch, the 3DKD's are computed and stored. The number of query patches is then further reduced so that each query patch remains with at least five neighboring target patches when both query and target proteins are considered superimposed. Here, by a neighbor, we mean that the physical distance between centers of the query patch and the target patch is less than 3 Å. So, the number of query patches to be compared with those in the target sets now depends on each individual target protein, and we obtain 737, 736, and 701 such query patches corresponding to Target 1, Target 2, and Target 3, respectively.




Each of the feature vectors in the query set is compared with those from the corresponding target set. The results are ranked using the Euclidean distance and collected for three target sets. If one of the five neighbors of the query patch ranks within top k, it is labeled as “correctly classified”. We computed recognition accuracies for k=1,…,16 where k=16 approximately corresponds to the top 1-percentile.

The results are shown in Fig. <ref>. Compared to the results of the Test I (Table <ref>), the accuracies dropped about 25 % points for the k=1 results. When we only look at the top result (k=1) in the rankings, K_5 is the most successful set of descriptors. The performances of K_3 and K_5 are comparable and higher than K_4 for almost all k values in all three cases. In Test II, we have actually tested how each 3DKD vector is tolerant to slight changes in the location of patch centers. From Fig. <ref>, it is clear that the descriptors in K_4, in particular, the 4th order descriptors that are not in K_3, are quite sensitive to such patch center shifts. This deficiency appears to be corrected in K_5 with the addition of 5th order descriptors.








 §.§ Comparison with 3D Zernike Descriptors


Next, we compare the discriminative performances of 3DKD and 3D Zernike descriptors (3DZD) <cit.> computed for local protein surface patches. To compute 3DZD's, we used the same query and Target 1 protein from Test I. The 3D Zernike descriptors are not able to extract local features from an object directly as 3DKD's do; yet Zernike functions are defined as continuous functions whose domain is the unit ball. For this reason, we first mapped each patch into the unit ball by considering each surface point as the patch center and placing a sphere of 6 Å radius around each point. 6 Å here corresponds to the size of a typical ligand-binding pocket on a protein surface. The patch cut from the surface is then mapped into the unit ball so that the center of mass of the cropped patch is placed at the coordinate origin. The geometric moments and hence the 3D Zernike moments of the order up to 12 and 15 (3DZD_12 and 3DZD_15, respectively) are computed for each patch using the algorithm provided in <cit.>. Finally, the rotation invariant 3DZD's are computed from these moments using the formula provided in <cit.>. The above work is repeated using a larger sphere for each patch (9 Å radius) since there are some well-known ligands binding to proteins by forming larger pockets. In order to make the comparison fairer, we have recomputed 3DKDs (K_5) using the same surface function but taking on zero value at voxels that remain outside the spheres used above. By this occlusion, we ensure that 3DKDs use the same local information as 3DZDs do.















The recognition accuracies for 3DKDs (K_5) and 3DZDs using different order of descriptors and patch radii are shown in Table <ref>. The results show that our method is clearly better than 3DZDs in local feature extraction, even when 3DZDs are allowed to use more invariants. The Top 1 prediction in our method is 90%, whereas it is only 7.4% in 3DZDs. Increasing the order of invariants to 15 (i.e., the vector size to 72) does not significantly improve the performance of 3DZD. Both methods give better recognition when the patch size is increased from 6 to 9 Å, but 3DKD outperforms 3DZD in all cases shown in Table <ref>.















 §.§ Binding Ligand Prediction














Finally, we test 3DKD on binding ligand prediction for proteins, which is one of the important tasks in bioinformatics as it addresses a central question in molecular biology, protein function <cit.>, and has real-life application in computational drug design <cit.>. Ligand molecules that bind to a local surface region in a protein can be predicted by finding similar local pockets of known binding ligands in the structure database. An example of a ligand binding pocket is demonstrated in Fig. <ref>. In order to test 3DKD on binding ligand prediction, we have constructed a benchmark dataset of 463 proteins already known to bind to 11 different ligands. See Fig. <ref> for these ligand structures. For each protein in the dataset, surface vertices and normals are generated using the GETPOINTS subroutine as in the previous tests. In the dataset, we also have the PDB file of each protein, which not only contains the coordinates of all atoms in the protein, but also those that belong to the bound ligand. For each atom in the ligand structure, we have selected the nearest surface vertex on the protein and annotated it with the bound ligand type. The collection of all such points and the 3DKDs of the patches around these points are all stored in a `patch database' of 11,039 patches together with their annotations of binding ligands. Thus, for each query pocket, a database search is performed for each of the patches in the query pocket, and a patch score is assigned to each patch based on the database rankings. 



We query total 497 pockets for this task. See Table <ref> for the number of query pockets for each ligand type. In Table <ref>, we have also listed the number of patches that corresponds to a ligand type. The number of patches associated with each pocket depends on the size of the pocket. By the average size of a pocket, we mean the radius of a sphere that encapsulates all surface vertices annotated with that ligand type. As can be seen from Table <ref>, there is a high correlation between the size of a typical pocket and the average number of patches included in that pocket. For each patch in the query pocket, we compute a patch score based on the following formula:

    (p,F,k)  =     ∑_i=1^k(δ_l(i),Flog(n/i))·∑_i=1^kδ_l(i),F∑_i=1^nδ_l(i),F,

where l(i) denotes the ligand type (e.g. AMP, ATP, etc.) of the i-th closest patch to the query, n is the number of patches in the patch database, and the function δ_l(i),F is equal to 1 if i-th patch is of type F, and 0 otherwise. This formula is used before as a pocket score for binding-ligand prediction in <cit.>. The first term in (<ref>) is to only involve k closest patches in the patch database to each patch from the query pocket, assigning a higher score to a patch with a higher rank. The second term is to normalize the score by the number of patches of the same type F included in the patch database so that the results are not biased in favor of highly populated ligand types in the database.

For each of the patches that belong to the query pocket and k values from 1 to 300, patch scores are computed and then summed to obtain a unique pocket score for each ligand type F as follows: 

    (P,F,k)  =     ∑_j=1^N_P(p_j,F,k),

where p_j, j=1,…, N_P are the patches within the pocket P, and N_P is the number of such patches. Thus, a pocket has a certain number of patches (see the last row of Table <ref>) and the score for a query pocket P for a ligand type F is computed as the sum of the score of each patch for the ligand F.

For each query pocket, we compute the pocket score for each of the 11 ligand type in the database. The ligand with the highest Pocket_score is predicted to bind to the query pocket. We then compare these 11 scores and look at the largest one (Top 1) and the largest three (Top 3) to obtain the number of successful predictions (see Table <ref>). For each ligand type, the number of successful cases is divided by the number of query pockets of that type in the pocket database to obtain prediction accuracies. For each prediction, the results are shown for the k value which maximizes the average prediction accuracy. It turns out that the average prediction accuracy is maximized for small values of k as shown in Table <ref>. 



When we look at the average prediction accuracies, K_5 performs best (with 41.0% correct prediction) among the 3DKDs. For the individual ligand types, K_3 performs best only for FUC and HEM, while K_4 gives the highest accuracies for AMP, ATP, FMN, FUC, GLC, and MAN. Among the 3DKDs, K_5 is the one giving the best average prediction and highest prediction accuracies (including ties) for almost all ligand types except for AMP, GLC, and HEM. In Table <ref>, we take a closer look at the results for K_5 by forming the confusion matrix with the true positives being along the diagonal. According to Table <ref>, seven AMP queries are predicted as ATP, and similarly, five ATP queries are predicted as AMP. This may be excused due to similar shapes of AMP and ATP (only differing from each other by two phosphate groups.) Other similar shaped pairs can be observed among sugar molecules FUC, GAL, GLC, and MAN. FUC as a query gives one false positive (GLC), while GAL as a query is confused seven times with other sugar molecules (three times with FUC and four times with GLC). GLC as a query is confused five times with FUC, three times with GAL, and three times with MAN. Similarly, MAN as a query gives twelve false positives from the sugar group (FUC, nine times; GAL, two times; and GLC once). Thus, the confusion between ligands by 3DKD is quite reasonable, capturing similarity of binding pockets of similar ligand molecules. In the results for Top 3, K_5 is still the one showing the best average performance among the 3DKDs, and giving the highest accuracies (including ties) for almost all ligand types except FAD and GLC.

In Table <ref>, we also provide a comparison of our approach with a former binding prediction method named Pocket-Surfer <cit.>. In Pocket-Surfer, 3DZDs are used as shape descriptors of a pocket. In addition to shape information, the pocket size is also utilized in Pocket-Surfer as a classification measure by an optimal weighted average of scores from both shape and pocket size. In our work, we only use shape information from 3DKDs without employing the pocket size in the scoring functions. When Top 1 is the classification criterion, 3DKD predicts better for seven ligand types and the average prediction. With Top 3 classification, Pocket-Surfer performs better than 3DKD for nine (out of eleven) ligand types. We believe that this is partly due to the fact that 3DKD employs shape information only, whereas Pocket-Surfer also used pocket size, which is often critical to distinguish ligands of the different sizes. We also show results with random prediction, in which each query pocket is scored based on a randomly shuffled pocket database, averaged over 3000 randomizations. It is clear that 3DKD outperforms the random prediction in all cases.













Table <ref> shows the time took for computing 3DKD of order up to 5 as outlined in steps (1)-(8) in Section <ref>. The programs were tested on two different platforms: a Windows laptop with i3 CPU of 2.53 GHz and 4 GB memory using MATLAB 2009b 64-bit, version 7.9, and on a Linux workstation with Xeon CPU of 3.60 GHz and 94 GB memory and using MATLAB 2016b 64-bit, version 9.1. On each computer, MATLAB is limited to a single computational thread to perform its jobs. The first row of Table <ref> shows the average CPU times spent on computing 3DKD of order up to 5, where the average is taken among 1608 local protein surface patches from Test I. The computation of 3DKDs can be done in half a second in a laptop, whereas it finishes twice faster on the PC. In the second row of Table <ref>, we also show the computational time from the binding prediction test, assuming that the 3DKDs of all 11,039 patches in the patch database are precomputed and stored. Given a typical query pocket, the time it takes to compute the pocket scores in (<ref>) for 11 ligand types and k=1,2,…, 10 is about 0.4 second in a laptop, and 3.7 times faster in the second platform.










§ CONCLUSION



In this paper, we have developed a novel set of local descriptors, three-dimensional Krawtchouk descriptors (3DKD), for identification and comparison of local regions of 3D voxelized oriented surfaces. Our approach is based on 3D Krawtchouk moments. While obtaining the rotation, size, and position independent invariants, we have preserved the critical ability of Krawtchouk moments to extract local features of a 3D surface from any region-of-interest. The locality property is due to the 3D weight function given in the definition of Krawtchouk polynomials. The weight contains three parameters p_x, p_y, and p_z, shifting the center of the local surface region along the x, y, and z-axes, respectively. We have noticed that, for each triplet (p_x,p_y,p_z), the coverage of the weight function is different, which prevents Krawtchouk moments from being translation and rotation invariant. To overcome these problems, we have computed the 3D weight function corresponding to (0.5,0.5,0.5) (center of the 3D grid) and used it for other local regions by translating the graph of the weight function as needed. To achieve rotational invariance, we utilized the surface normal at the vertex of the local surface region and computed the eigenvectors of the local inertia matrix. Among the eight possible different orientations, we have chosen the one positioning the vertex normal in a fixed octant in 3D space. We have also provided a detailed scheme for efficient computation of 3D Krawtchouk descriptors. 

We have tested the discriminative performances of 3DKD on three test problems. For each test, we have used K_3, K_4, and K_5, namely 3DKD of the order up to 3, 4, and 5, respectively. In the first test, the results have been comparable, while K_5 has the best recognition accuracies for predicting the top match. The second test demonstrated that K_5, among the 3DKD, is the most robust set of descriptors to small changes in patch location. We have also compared 3DKD with 3DZD. 3DKD shows much better recognition accuracies than 3DZD in all cases reported. As the third test, we have employed 3DKD for prediction of ligand binding sites on protein surfaces. 3DKD showed better performance than Patch-Surfer when Top 1 prediction was considered. From the results of the second and the third tests suggest that 3DKD is more sensitive than 3DZD in detecting subtle shape similarity. The results on the binding ligand prediction were obtained by only considering geometric shape information of protein surface. Therefore further improvement is expected by integrating other features, such as the electrostatic potential or other physicochemical properties, to characterize protein surface regions, which is analogous to representing color images rather than black-and-white images of protein surfaces.







§ ACKNOWLEDGMENTS

This work was supported by the National Science Foundation (DMS1614777 and DMS1614661). DK also acknowledges support from National Institutes of Health, (R01GM123055) and the National Science Foundation (CMMI1825941) and the Purdue Institute for Drug Discovery. The authors would like to thank Xiaolei Zhu for the help in voxelizing protein surfaces and Juan Esquivel-Rodriguez for the help in running LZerD software.



















00




ShuEtAl2007
H. Shu, L. Luo, and J. L. Coatrieux, Moment-based approaches in imaging. Part 1, basic features, IEEE engineering in medicine and biology magazine: the quarterly magazine of the Engineering in Medicine & Biology Society, 26(5):70-74, 2007.

Hu1962
M. K. Hu, Visual pattern recognition by moment invariants, Information Theory, IRE Transactions on, 8(2), 179-187, 1962.

SadjadiHall1980
F. A. Sadjadi and E. L. Hall, Three-dimensional moment invariants, IEEE Transactions on Pattern Analysis and Machine Intelligence, 2:127-136, 1980.

Guo1993
X. Guo, Three dimensional moment invariants under rigid transformation, In Computer Analysis of Images and Patterns, pp. 518-522. Springer Berlin/Heidelberg, 1993.

LoDon1989
C-H. Lo and H-S. Don. 3-D moment forms: their construction and application to object identification and positioning, IEEE Transactions on Pattern Analysis and Machine Intelligence 11, no. 10:1053-1064, 1989.

GalvezCanton1992
J. M. Galvez and M. Canton, Normalization and shape recognition of three-dimensional objects by 3D moments, Pattern Recognition, 26(5), 667-681, 1993.

Mamistvalov1998
A. G. Mamistvalov, n-dimensional moment invariants and conceptual mathematical theory of recognition n-dimensional solids, IEEE Transactions on pattern analysis and machine intelligence 20, no. 8: 819-831, 1998.

FlusserEtAl2003
J, Flusser, J. Boldys, and B. Zitova, Moment forms invariant to rotation and blur in arbitrary number of dimensions, IEEE Transactions on Pattern Analysis and Machine Intelligence 25, no. 2:234-246, 2003.

MehtaEtAl2009
N. Mehta, R. S. Alomari, V. Chaudhary, Content based sub-image retrieval system for high resolution pathology images using salient interest points, Conference proceedings: Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 3719-3722, 2009.

SitKihara2014
A. Sit and D. Kihara, Comparison of Image Patches Using Local Moment Invariants, IEEE Transactions on Image Processing, 23(5), 2369-2379, 2014.

Pryzva1992
G. Y. Pryzva, Kravchuk orthogonal polynomials, Ukrainian Mathematical Journal, 44(7), 792-800, 1992.

YapRaveendranOng2002
P. T. Yap, P. Raveendran, S. H. Ong, Krawtchouk moments as a new set of discrete orthogonal moments for image reconstruction, In Neural Networks, 2002. IJCNN'02. Proceedings of the 2002 International Joint Conference on (Vol. 1, pp. 908-912). IEEE, 2002.

YapParamesranOng2003
P. T. Yap, R. Paramesran, S. H. Ong, Image analysis by Krawtchouk moments, IEEE Transactions on Image Processing, 12(11), 1367-1377, 2003.

VenkataramanaRaj2007
A. Venkataramana and P. A. Raj, Image watermarking using Krawtchouk moments In International Conference on Computing: Theory and Applications, 2007, ICCTA'07, (pp. 676-680), IEEE, 2007.

JassimRaveendran2012
W. A. Jassim and P. Raveendran, Face Recognition Using Discrete Tchebichef-Krawtchouk Transform, In Multimedia (ISM), 2012 IEEE International Symposium on (pp. 120-127), IEEE, 2012.

RaniDevaraj2012
J. S. Rani and D. Devaraj, Face recognition using Krawtchouk moment, Sadhana, 37(4), 441-460, 2012.

MademlisEtAl2006
A. Mademlis, A. Axenopoulos, P. Daras, D. Tzovaras, and M. G. Strintzis, 3D Content-Based Search Based on 3D Krawtchouk Moments, In 3DPVT (pp. 743-749), 2006.

DarasEtAl2006
P. Daras, D. Tzovaras, A. Mademlis, A. Axenopoulos, D. Zarpalas, and M. G. Strintzis, 3D search and retrieval using Krawtchouk moments, SHREC2006, 17-21, 2006.

XiangEtAl2006
P. Xiang, C. Qihua, and L. Zhi, Content-Based 3d retrieval by krawtchouk moments, In Image Analysis and Recognition (pp. 217-224), Springer Berlin Heidelberg, 2006.

MesbahEtAl2016
A. Mesbah, M. El Mallahi, Z. Lakhili, H. Qjidaa, and A. Berrahou, Fast and accurate algorithm for 3D local object reconstruction using Krawtchouk moments, 2016 5th International Conference on Multimedia Computing and Systems (ICMCS), Marrakech, pp. 1-6, 2016.

FlusserEtAl2009
J. Flusser, B. Zitova, and T. Suk, Moments and moment invariants in pattern recognition, John Wiley & Sons, 2009.

NovotniKlein2003
M. Novotni and R. Klein, 3D Zernike descriptors for content based shape retrieval, In Proceedings of the eighth ACM symposium on solid modeling and applications, pp. 216-225. ACM, 2003.

FlusserEtAl2016
J. Flusser, T. Suk, and B. Zitova, 2D and 3D Image Analysis by Moments, John Wiley & Sons, 2016.

ChikhiEtAl2011
R. Chikhi, L. Sael, and D. Kihara, Protein binding ligand prediction using moments-based methods, In Protein Function Prediction for Omics Era, pp. 145-163. Springer Netherlands, 2011.

SealKihara2012
L. Sael and D. Kihara, Detecting local ligand-binding site similarity in nonhomologous proteins by surface patch comparison, Proteins, 80(4):1177-1195, 2012.

ZhuEtAl2015
X. Zhu, Y. Xiong, and D. Kihara, Large-scale binding ligand prediction by improved patch-based method Patch-Surfer2.0, Bioinformatics, 31(5):707-713, 2015.

RosenbergGoldblum2006
M. Rosenberg and A. Goldblum, Computational protein design: a novel path to future protein drugs, Current pharmaceutical design, 12(31), 3973-3997, 2006.

ChikhiSaelKihara2010
R. Chikhi, L. Sael, and D. Kihara, Real-time ligand binding pocket database search using local surface descriptors, Proteins: Structure, Function, and Bioinformatics, 78(9), pp. 2007-2028, 2010.

ZhangEtAl2010
G. Zhang, Z. Luo, B. Fu, B. Li, J. Liao, X. Fan, and Z. Xi, A symmetry and bi-recursive algorithm of accurately computing Krawtchouk moments, Pattern Recognition Letters, 31(7), 548-554, 2010.

YapParamesranOng2007
P. T. Yap, R. Paramesran, and S. H. Ong, Image analysis using Hahn moments, IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(11), 2057-2062, 2007.

PrincetonShapeBenchmark2004
P. Shilane, P. Min, M. Kazhdan, and T. Funkhouser, The Princeton shape benchmark, In Shape modeling applications, 2004. Proceedings, pp. 167-178. IEEE, 2004.

PatilRavi2005
S. Patil and B. Ravi, Voxel-based representation, display and thickness analysis of intricate shapes, In Computer Aided Design and Computer Graphics, 2005. Ninth International Conference on, pp. 6-pp. IEEE, 2005.

PDB2000
H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov, and P. E. Bourne, The protein data bank, Nucleic acids research 28, no. 1:235-242, 2000.

3DSURFER
D. La, J. Esquivel-Rodriguez, V. Venkatraman, B. Li, L. Sael, S. Ueng, S. Ahrendt, D. Kihara, 3D-SURFER: Software for high throughput protein surface comparison and analysis, Bioinformatics 25: 2843-2844, 2009.

LZERD
V. Venkatraman, Y. D. Yang, L. Sael, and D. Kihara, Protein-protein docking using region-based 3D Zernike descriptors, BMC Bioinformatics 10:407, 2009.

SitEtAl2013
A. Sit, J. C. Mitchell, G. N. Phillips, S. J. Wright, An extension of 3D Zernike moments for shape description and retrieval of maps defined in rectangular solids, Computational and Mathematical Biophysics (formerly Molecular Based Mathematical Biology) 1:75-89, 2013.

UCSFChimera
E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng, T. E. Ferrin, UCSF Chimera – a visualization system for exploratory research and analysis, J Comput Chem. Oct;25(13):1605-12, 2004.

PLPS2
W. H. Shin, C. W. Christoffer, J. Wang, D. Kihara, PL-PatchSurfer2: Improved local surface matching-based virtual screening method that is tolerant to target and ligand structure variation, J. Chem. Inf. Model. 56(9):1676-91, 2016.

Patch
W. H. Shin, M. G. Bures, D. Kihara, PatchSurfers: Two methods for local molecular property-based binding ligand prediction, Methods 93:41-50, 2016.

LZERDSuite
J. Esquivel-Rodoriguez, V. Filos-Gonzalez, B. Li, D. Kihara, Pairwise and multimeric protein-protein docking using the LZerD program suite, Methods Mol. Biol. 1137:209-34, 2014.





Atilla Sit is currently an Assistant Professor with the Department of Mathematics and Statistics, Eastern Kentucky University, Richmond, KY, USA. He received the B.S. degree from Middle East Technical University, Ankara, Turkey, the M.S. degree from Bogazici University, Istanbul, Turkey, and the Ph.D. degree from Iowa State University, Ames, IA, USA, in 2001, 2005, and 2010, respectively. His current research interests include image/structure analysis, orthogonal systems, special functions, protein structure determination and comparison, and protein classification.



Daisuke Kihara is currently a Professor with the Department of Biological Sciences and the Department of Computer Science, Purdue University, West Lafayette, IN, USA. He received the B.S. degree from the University of Tokyo, Japan, in 1994, and M.S. and Ph.D. degrees from Kyoto University, Japan, in 1996 and 1999, respectively. His research projects include biomolecular shape comparison, computational drug design, protein tertiary structure prediction, protein-protein docking, and other geometrical problems in computational biology. He has authored over 150 research papers and book chapters. He is named Showalter University Faculty Scholar from Purdue University in 2013.


