Tensor template matching

Abstract

A method of template matching a tensor template of a target particle to identify an instance of the target within a tomographic image includes obtaining a tensor template of the target particle, wherein the tensor template describes a shape of the target particle under all rotations using a tensor field; and identifying an instance of the target particle within the tomographic image using the tensor template.

Description

Technical Field

[0001] This invention relates to an improved method for template matching in image analysis. Specifically, it relates to template matching of target molecules in cellular cryo-electron tomography, which can be performed faster than before. Background Technology

[0002] This invention relates to template matching, in which an image is scanned to determine a match with a target template. The invention can be applied to both 2D and 3D images (volume or tomographic images), although it is particularly advantageous in the case of 3D images. The invention can be used for template matching of targets that are present multiple times in an image, such as particles within the image, like ribosomes or proteasomes. More generally, template matching is widely used in industries such as image registration, alignment, object detection, counting, tracking, defect detection, volume reconstruction, etc. It is commonly used in such applications due to its ease of use and interpretation.

[0003] One specific application of this invention is cryo-electron tomography (Cryo-ET), a label-free imaging technique that provides 3D images (or tomographic images) of organelles and protein complexes at nanometer resolution. This is accomplished by opening a window within the cell using a focused ion beam milling technique applied to cryogenically frozen cells. A series of 2D images of this thinned cell sample are taken from different angles. The 2D images are then reconstructed into 3D tomographic images. This technique provides a 3D snapshot of proteins functioning within their functional environment and offers a window into how they work with other molecules to carry out key processes within the cell. Cryo-ET is used for two main purposes: protein structure discovery and functional cell behavior. In both cases, template matching is typically used to automatically pick up particles of interest to obtain their location, orientation, and class.

[0004] Template matching is a brute-force algorithm that computes the correlation between the template of the target particle and the template at each voxel in the tomographic image. The voxel showing the highest correlation value corresponds to the instance of the template position in the tomographic image that best matches the particle recorded in the tomographic image. Furthermore, because each instance of the template position can have any rotation of the target particle, the entire space of the target particle's rotations must be sampled, and the correlation must be computed for all rotations to find the optimal rotation at each voxel in the tomographic image. This method is computationally intensive and requires a significant amount of computation time even after algorithm optimization, the use of Fourier transforms, parallelization, and data downsampling. This constitutes a serious bottleneck in the Cryo-ET analysis workflow. Summary of the Invention

[0005] This paper describes a method for generating tensor templates for template matching to determine instances of targets within an image; the method includes: generating tensor templates of targets, wherein the tensor templates utilize tensor fields to describe the shape of the targets under all rotations.

[0006] Optionally, the tensor template comprises a tensor field of multiple tensors. Depending on the trade-off between the rotational accuracy and computational complexity required by the target application, the tensors can be of any order and any dimension. For example, 4th-order and 4-dimensional tensors can be used.

[0007] The target can be described by a model comprising an array of elements, such as pixels in a 2D image or voxels in a 3D image. Tensor templates can include per-element tensors of the model. Templates can be generated for different rotations of the target, each template comprising an array of elements. Optionally, tensor powers are generated for each rotation of the target, and the tensor powers can include components describing the shape of the target under that rotation. Each tensor power can be computed as a tensor power of a quaternion describing the shape of the target under that rotation. Each tensor power can be computed as a fourth-power multiplication of the four components of the quaternion.

[0008] Preferably, all tensors and all tensor powers have the same number of components, and each component of each tensor representing an element of the target comprises the integral over all rotations of the product of the corresponding component of the tensor power for one rotation and the value of an element of the template of the target for that one rotation. The value of each element may represent the intensity of that element. For example, the image may include a surface model of the target that uses intensity values ​​for each element to describe the shape of the target's surface.

[0009] Optionally, the method includes masking the model using an operator that defines the model of the target before generating the components of the tensor, thereby setting the size of the element array. This operator can be a symmetric positive semidefinite linear operator.

[0010] According to a first aspect of the present invention, a method for tensor template matching for identifying instances of a target within an image is provided, the method comprising: generating a tensor template of the target according to any of the methods described above; and using the tensor template to identify instances of the target within the image.

[0011] According to a second aspect, a method is provided for template matching of a target template to identify instances of the target within an image, the method comprising: obtaining a tensor template of the target, wherein the tensor template utilizes a tensor field to describe the shape of the target under all rotations; and using the tensor template to identify instances of the target within the image. Obtaining the tensor template of the target may include generating the tensor template or retrieving a tensor template from a stored tensor template library.

[0012] Further optional features of the first and second aspects of the invention will now be described.

[0013] Optionally, instances of using a tensor template to identify targets within an image include: convolving the tensor template with the image to obtain a convolutional tensor field that provides information about the correlation between the tensor template and the image at all locations within the image and under all rotations. Instances of using a tensor template to identify targets within an image may include determining the optimal correlation value for each location within the image by computing a scalar map of the convolutional tensor field at each location. The values ​​filling the scalar map can be a measure of the optimal correlation value at each location. Each location may correspond to an element of the image, such as a pixel in a two-dimensional image or a voxel in a three-dimensional image. Compiling the scalar map of the convolutional tensor field at each location may include computing a scalar map of the total energy of the convolutional tensor field at each location. The method may include computing a scalar map of the total energy of the convolutional tensor field at each location using the norm of the convolutional tensor field at each location. The Frobenius norm may be used. Compiling the Frobenius norm of the convolutional tensor field at each location may include computing the square root of the absolute square of the components of the tensor of the tensor template at each location.

[0014] The method may also include identifying the location of instances of a target within an image by determining the local maximum of the optimal correlation value. Identifying the rotation of the target at each identified location may include using global optimization techniques.

[0015] Identifying the rotation of the target at each of the identified locations may include identifying the eigenvectors and eigenvalues ​​of the convolutional tensor field at each of the identified locations. Identifying the rotation of the target at each identified location by identifying the eigenvectors and eigenvalues ​​of the tensor in the convolutional tensor field at each identified location may include unfolding the tensor in the convolutional tensor field into a matrix and performing operations on that matrix to determine the eigenvectors and eigenvalues. Performing operations on the matrix to determine the eigenvectors and eigenvalues ​​includes using several calls to a shift-symmetric higher-order power method that restarts at different random rotations.

[0016] Identifying the rotations of a target at each of the identified locations using global optimization techniques may involve sampling the rotation space of the target, represented as a quaternion, and computing the inner product between a selected tensor of the convolutional tensor field and the quaternion. This provides a list of rotations sorted by the value of this product. A list of peak locations with their most promising rotations (highest inner products) can be returned for each peak.

[0017] Optionally, the method may also include using global optimization techniques to identify the target's position and rotation. Position and rotation can be identified in a single step, rather than sequentially as may be done with the methods described above. The target's position and rotation can be identified by identifying the eigenvectors and eigenvalues ​​of the convolutional tensor field. Identifying the target's position and rotation by identifying the eigenvectors and eigenvalues ​​of the tensors in the convolutional tensor field may include unfolding the tensors in the convolutional tensor field into a matrix and performing operations on that matrix to determine the eigenvectors and eigenvalues. Performing operations on the matrix to determine the eigenvectors and eigenvalues ​​may include using several calls to a shift-symmetric higher-order power method that restarts at different random rotations. Identifying the rotation of the target at each of the identified positions may include sampling the rotation space of the target, represented as a quaternion, and computing the inner product between a selected tensor of the convolutional tensor field and the quaternion. This provides a list of rotations that can be sorted by the values ​​of the product. A list of peak positions with their most promising rotations (highest inner products) for each peak can be returned.

[0018] Preferably, using tensor templates to identify instances of a target within an image includes: using the identified position and rotation for each instance of the target within the image as input, and calculating a normalized correlation value for each identified instance of the target within the image by means of a non-tensor template of the target and a locally normalized cross-correlation of the image. Advantageously, this step produces normalized cross-correlation values ​​that allow for better comparison of instances identified for position and rotation. While calculating locally normalized cross-correlation values ​​is typically very time-consuming, this is not the case in this case because only values ​​for the identified position and rotation need to be calculated, which significantly reduces the computational space (compared to the need to calculate for all positions and all rotation directions in conventional locally normalized cross-correlation). The method may include generating each non-tensor template from a model of the target as a template describing the shape of the target at the corresponding identified rotation.

[0019] Optionally, the method further includes providing a list of identified instances of a target within the image. For each instance of a target, the list may include a calculated locally normalized cross-correlation value, the identified location, and the identified rotation. The list may be limited according to user-selected preferences.

[0020] The method may further include segmenting the image into blocks. For each block, the method may include convolving a tensor template with that block of the image to obtain a convolutional tensor field that provides information about the correlation between the tensor template and the image at all locations and under all rotations within that block of the image. Then, an optimal correlation value can be determined for each location within that block of the image, thereby computing a scalar map from the convolutional tensor field. The location of an instance of a target within that block of the image can then be identified by determining the local maxima of the optimal correlation value. Finally, a global optimization technique can be used to identify the rotation of the target at each of the identified locations.

[0021] Therefore, the previously described steps of convolving the tensor template, determining the relevant values ​​at each location, and identifying the position and rotation of each instance of the target in the image are broken down as a whole into a series of routines executed sequentially for each block. This can be advantageous when computer memory is limited.

[0022] Segmenting an image into blocks may include dividing the image into contiguous units, and defining each block as comprising one unit plus an overlapping portion extending into adjacent units. The overlap may be a uniform distance, and the size of the overlapping portion extending into adjacent units may be less than the maximum size of the target and may be greater than half the test size of the target. The method may include, for each block, retaining locations identified within the units of that block and discarding any locations identified within the overlapping portion of that block.

[0023] The method may also include incorporating rotations and positions identified in each block of the image before computing the normalized cross-correlation value calculated for each identified instance of a target within the image.

[0024] Optionally, the target is a particle, and the image is a tomographic image. Tomographic images can be obtained using cellular cryo-electron tomography. Alternatively, the target can be a defect, and the image can be an image of a manufactured part. The defect can be a defect in a semiconductor material, and the manufactured part can be a semiconductor component such as a chip or wafer.

[0025] A method for counting objects within an image is also described, comprising any of the methods described above, and providing the number of identified instances of the target in the image as output. A method for tracking objects within a series of images is also described, comprising performing any of the methods described above to identify instances of the target in each image of the series.

[0026] The present invention also relates to a computer program including computer program instructions that, when executed by a computer processor, cause the computer processor to perform any of the methods described above; a computer-readable medium on which such a computer program is stored; and a computer system including a computer processor and a computer memory in which such a computer program is stored.

[0027] The following clauses describe further features.

[0028] 1. A method for generating a tensor template for template matching to determine an instance of a target within an image; the method comprising: generating a tensor template of the target, wherein the tensor template utilizes a tensor field to describe the shape of the target under all rotations.

[0029] 2. The method according to Clause 1, wherein the tensor template comprises a tensor field of a plurality of tensors.

[0030] 3. The method according to Clause 2, wherein the objective is described by a model comprising an array of elements, and the tensor template comprises per-element tensors of the model.

[0031] 4. The method according to Clause 3, wherein templates are generated for different rotations of the target, each template comprising the array of elements.

[0032] 5. The method according to Clause 4, wherein a tensor power is generated for each rotation of the target, and wherein the tensor power includes a component describing the shape of the target under that rotation.

[0033] 6. The method according to Clause 5, wherein each tensor power is calculated as a tensor power of a quaternion describing the shape of the target under that rotation.

[0034] 7. The method according to Clause 6, wherein each tensor power is calculated as the fourth power multiplication of the four components of the quaternion.

[0035] 8. The method according to any one of clauses 5 to 7, wherein all tensors and all tensor powers have the same number of components, and each component of each tensor representing an element of the target comprises the integral over all rotations of the product of the corresponding component of the tensor power for one rotation and the value of the element of the template of the target for the one rotation.

[0036] 9. The method according to Clause 8, wherein the value of each element represents the intensity of that element.

[0037] 10. The method according to any one of clauses 3 to 9, the method comprising masking the model using an operator defining the model of the target before generating the components of the tensor, thereby setting the size of the array of elements.

[0038] 11. The method according to Clause 10, wherein the operator is a symmetric positive semidefinite linear operator.

[0039] 12. A method for tensor template matching for identifying instances of a target within an image, the method comprising: generating a tensor template of the target according to any one of claims 1 to 11; and using the tensor template to identify instances of the target within the image.

[0040] 13. A method for template matching of a target to identify instances of the target within an image, the method comprising: obtaining a tensor template of the target, wherein the tensor template describes the shape of the target under all rotations by means of a tensor field; and using the tensor template to identify instances of the target within an image.

[0041] 14. The method according to Clause 13, wherein obtaining the tensor template of the target includes generating the tensor template or retrieving the tensor template from a stored tensor template library.

[0042] 15. The method according to any one of clauses 12 to 14, wherein an example of using the tensor template to identify the target within an image includes: The tensor template is convolved with the image to obtain a convolutional tensor field, which provides information about the correlation between the tensor template and the image at all positions and under all rotations within the image.

[0043] 16. The method according to Clause 15, wherein an instance of using the tensor template to identify the target within an image comprises: determining the optimal correlation value for each location within the image by computing a scalar map of the convolutional tensor field at each location.

[0044] 17. The method according to Clause 16, wherein calculating a scalar graph of the convolution tensor field at each location includes calculating a scalar graph of the total energy of the convolution tensor field at each location.

[0045] 18. The method according to Clause 17, the method comprising using the norm of the convolutional tensor field at each location, optionally the Frobenius norm, to compute a scalar graph of the total energy of the convolutional tensor field at each location, and optionally, wherein computing the Frobenius norm of the convolutional tensor field at each location comprises computing the square root of the absolute square of the components of the tensor of the tensor template at each location.

[0046] 19. The method according to any one of clauses 16 to 18, the method further comprising: identifying the location of an instance of the target within the image by determining a local maximum value of the optimal correlation value.

[0047] 20. The method according to Clause 19, further comprising: using a global optimization technique to identify rotations of the target at each of the identified locations.

[0048] 21. The method according to clause 19 or 20, the method further comprising: identifying rotation of the target at each of the identified locations by identifying eigenvectors and eigenvalues ​​of the convolution tensor field at each of the identified locations.

[0049] 22. The method according to Clause 21, wherein identifying the rotation of the target at each identified position by identifying the eigenvectors and eigenvalues ​​of the tensors in the convolutional tensor field at each identified position comprises unfolding the tensors in the convolutional tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0050] 23. The method according to Clause 22, wherein performing operations on the matrix to determine the eigenvectors and eigenvalues ​​includes using several calls to a shift-symmetric higher-order power method that restarts under different random rotations.

[0051] 24. The method according to Clause 20, further comprising: identifying rotations of the target at each of the identified locations by sampling the rotation space of the target, represented as a quaternion, and calculating the inner product between a selected tensor of the convolution tensor field and the quaternion.

[0052] 25. The method according to Clause 15, further comprising: using a global optimization technique to identify the position and rotation of the target.

[0053] 26. The method according to Clause 27, further comprising: identifying the position and rotation of the target by identifying eigenvectors and eigenvalues ​​of the convolution tensor field.

[0054] 27. The method according to Clause 26, wherein identifying the position and rotation of the target by recognizing the eigenvectors and eigenvalues ​​of the tensors in the convolutional tensor field comprises unfolding the tensors in the convolutional tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0055] 28. The method according to Clause 27, wherein performing operations on the matrix to determine the eigenvectors and eigenvalues ​​includes using several calls to a shift-symmetric higher-order power method that restarts under different random rotations.

[0056] 29. The method according to Clause 19 or 20, the method further comprising: identifying rotations of the target at each of the identified locations by sampling the rotation space of the target, represented as a quaternion, and calculating the inner product between a selected tensor of the convolution tensor field and the quaternion.

[0057] 30. The method according to any one of clauses 20 to 29, wherein using the tensor template to identify instances of the target within an image comprises: calculating a normalized correlation value for each identified instance of the target within the image by using the position and rotation identified for each instance of the target within the image as input to calculate the local normalized cross-correlation between the non-tensor template of the target and the image.

[0058] 31. The method according to Clause 30, wherein calculating a normalized correlation value for each identified instance of the target within the image comprises: calculating the non-tensor template of the target and the local normalized cross-correlation of the image with respect to the position and rotation identified only for the instance of the target within the image.

[0059] 32. The method according to clause 30 or 31, the method further comprising: generating each non-tensor template based on a model of the target as a template describing the shape of the target under the corresponding identified rotation.

[0060] 33. The method according to any one of clauses 20 to 32, the method further comprising providing a list of identified instances of the target within the image; and optionally, wherein: for each instance of the target, the list includes a calculated locally normalized cross-correlation value, an identified position, and an identified rotation; and / or the list is limited according to user-selected preferences.

[0061] 34. The method according to any one of claims 12 to 33, the method further comprising segmenting the image into blocks; and, for each block, sequentially: (i) convolving the tensor template with the block of the image to obtain a convolutional tensor field, the convolutional tensor field providing information about the correlation between the tensor template and the image at all positions and under all rotations within the block of the image; (ii) determining the optimal correlation value at each of the positions within the block of the image, thereby computing a scalar map from the convolutional tensor field; (iii) identifying the position of an instance of the target within the block of the image by determining the local maximum of the optimal correlation value; and (iv) identifying the rotation of the target at each of the identified positions using a global optimization technique.

[0062] 35. The method according to clause 34, wherein segmenting the image into blocks comprises segmenting the image into consecutive units, and defining each block as comprising a unit plus an overlapping portion extending into the adjacent unit.

[0063] 36. The method according to Clause 35, wherein the overlap is a uniform distance, and optionally, the size of the overlapping portion extending into adjacent cells is less than the maximum size of the target and greater than half the test size of the target.

[0064] 37. The method according to Clause 36, the method comprising: for each block, retaining the locations identified within the cells of the block, and discarding any locations identified within the overlapping portions of the block.

[0065] 38. The method according to any one of clauses 34 to 37, the method further comprising, prior to calculating the normalized correlation value calculated for each identified instance of the target within the image, incorporating the rotation and position identified in each block of the image.

[0066] 39. The method according to any of the foregoing clauses, wherein the target is a particle and the image is a tomographic image, optionally a tomographic image obtained using cellular cryo-electron tomography.

[0067] 40. The method according to any one of clauses 1 to 38, wherein the target is a defect and the image is an image of a manufactured part.

[0068] 41. The method according to Clause 40, wherein the defect is a defect in a semiconductor material and the manufactured part is a semiconductor part, such as a chip or wafer.

[0069] 42. A method for counting objects within an image, the method comprising the method according to any one of clauses 12 to 41, and providing as output the number of identified instances of the object in the image.

[0070] 43. A method for tracking objects within a series of images, the method comprising performing a method according to any one of claims 12 to 41 to identify instances of the object in each of the series of images.

[0071] 44. A computer program comprising computer program instructions that, when executed by a computer processor, cause the computer processor to perform the method described in accordance with any of the foregoing provisions; a computer-readable medium on which such computer program is stored; or a computer system comprising a computer processor and a computer memory in which such computer program is stored.

[0072] 101. A method for template matching of a tensor template of a target particle to identify instances of the target within a tomographic image, the method comprising: obtaining a tensor template of the target particle, wherein the tensor template utilizes a tensor field to describe the shape of the target particle under all rotations; and using the tensor template to identify instances of the target particle within the tomographic image.

[0073] 102. The method according to Clause 101, wherein obtaining the tensor template of the target particle comprises retrieving the tensor template from a stored tensor template library.

[0074] 103. The method according to Clause 101, wherein obtaining a tensor template of the target particle includes generating the tensor template, wherein the tensor template utilizes a tensor field to describe the shape of the target particle over all rotations.

[0075] 104. The method according to Clause 103, wherein the tensor template comprises a tensor field of a plurality of tensors.

[0076] 105. The method according to Clause 104, wherein the target particle is described by a model comprising an array of voxels, and the tensor template comprises a per-voxel tensor of the model.

[0077] 106. The method according to Clause 105, wherein templates are generated for different rotations of the target particle, each template comprising the voxel array.

[0078] 107. The method according to Clause 106, wherein a tensor power is generated for each rotation of the target particle, and wherein the tensor power includes a component describing the shape of the target particle under that rotation.

[0079] 108. The method according to Clause 107, wherein each tensor power is calculated as a tensor power of a quaternion describing the shape of the target particle under that rotation.

[0080] 109. The method according to Clause 108, wherein each tensor power is calculated as a fourth power multiplication of the four components of the quaternion.

[0081] 110. The method according to any one of clauses 107 to 109, wherein all tensors and all tensor powers have the same number of components, and each component of each tensor representing a voxel of the target particle comprises the integral over all rotations of the product of the corresponding component of the tensor power for one rotation and the value of the voxel of the template of the target particle for the one rotation.

[0082] 111. The method according to Clause 110, wherein the value of each voxel represents the intensity of that voxel.

[0083] 112. The method according to any one of clauses 105 to 111, the method comprising masking the model using an operator defining the model of the target particle, thereby setting the size of the element array, prior to generating the components of the tensor.

[0084] 113. The method according to Clause 112, wherein the operator is a symmetric positive semidefinite linear operator.

[0085] 114. The method according to any one of clauses 101 to 113, wherein an example of using the tensor template to identify the target particle within the tomographic image comprises: convolving the tensor template with the tomographic image to obtain a convolution tensor field, the convolution tensor field providing information about the correlation between the tensor template and the tomographic image at all locations and under all rotations within the tomographic image.

[0086] 115. The method according to Clause 114, wherein an instance of using the tensor template to identify the target particle within the tomographic image comprises: determining the optimal correlation value for each of the locations within the tomographic image by computing a scalar map of the convolution tensor field at each location.

[0087] 116. The method according to Clause 115, wherein calculating a scalar graph of the convolution tensor field at each location includes calculating a scalar graph of the total energy of the convolution tensor field at each location.

[0088] 117. The method according to Clause 116, the method comprising using the norm of the convolutional tensor field at each location, optionally the Frobenius norm, to compute a scalar graph of the total energy of the convolutional tensor field at each location.

[0089] 118. The method according to Clause 117, the method comprising calculating a scalar graph using the Frobenius norm, Calculating the Frobenius norm of the convolutional tensor field at each location includes calculating the square root of the absolute square of the components of the tensor template at each location.

[0090] 119. The method according to any one of clauses 115 to 118, the method further comprising: The location of an instance of the target particle within the tomographic image is identified by determining the local maximum of the optimal correlation value.

[0091] 120. The method according to Clause 119, the method further comprising: using a global optimization technique to identify the rotation of the target particle at each of the identified positions.

[0092] 121. The method according to Clause 1119 or 20, the method further comprising: identifying rotation of the target particle at each identified position by identifying eigenvectors and eigenvalues ​​of the convolution tensor field at each identified position.

[0093] 122. The method according to Clause 121, wherein identifying the rotation of the target particle at each identified position by identifying the eigenvectors and eigenvalues ​​of the tensors in the convolutional tensor field at each identified position comprises unfolding the tensors in the convolutional tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0094] 123. The method according to Clause 122, wherein performing operations on the matrix to determine the eigenvectors and the eigenvalues ​​includes using several calls to a shift-symmetric higher-order power method that restarts under different random rotations.

[0095] 124. The method according to Clause 120, the method further comprising: identifying rotations of the target at each of the identified locations by sampling the rotation space of the target particle, represented as a quaternion, and calculating the inner product between a selected tensor of the convolution tensor field and the quaternion.

[0096] 125. The method according to Clause 115, further comprising: using a global optimization technique to identify the position and rotation of the target.

[0097] 126. The method according to Clause 127, further comprising: identifying the position and rotation of the target particle by identifying the eigenvectors and eigenvalues ​​of the convolution tensor field.

[0098] 127. The method according to Clause 126, wherein identifying the position and rotation of the target particle by identifying the eigenvectors and eigenvalues ​​of the tensors in the convolutional tensor field comprises unfolding the tensors in the convolutional tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0099] 128. The method according to Clause 127, wherein performing operations on the matrix to determine the eigenvectors and the eigenvalues ​​includes using several calls to a shift-symmetric higher-order power method that restarts under different random rotations.

[0100] 129. The method according to clause 119 or 120, the method further comprising: identifying the rotation of the target particle at each of the identified locations by sampling the rotation space of the target, represented as a quaternion, and calculating the inner product between a selected tensor of the convolution tensor field and the quaternion.

[0101] 130. The method according to any one of claims 120 to 129, wherein using the tensor template to identify instances of the target particle within the tomographic image comprises: calculating a normalized correlation value for each identified instance of the target particle within the tomographic image by using the position and rotation identified for each instance of the target particle within the tomographic image as input to calculate a local normalized cross-correlation between the non-tensor template of the target particle and the tomographic image.

[0102] 131. The method according to Clause 130, wherein calculating the normalized correlation value of each identified instance of the target particle within the tomographic image comprises: calculating the non-tensor template of the target particle and the local normalized cross-correlation of the tomographic image with respect to the position and rotation identified only for the instance of the target particle within the tomographic image.

[0103] 132. The method according to clause 130 or 131, the method further comprising generating each non-tensor template based on a model of the target particle as a template describing the shape of the target particle under the corresponding identified rotation.

[0104] 133. The method according to any one of clauses 120 to 132, the method further comprising providing a list of identified instances of the target particle within the tomographic image; and optionally, wherein: for each instance of the target particle, the list includes a calculated local normalized cross-correlation value, an identified location, and an identified rotation; and / or the list is limited according to user-selected preferences.

[0105] 134. The method according to any one of claims 101 to 133, the method further comprising segmenting the tomographic image into blocks; for each block, sequentially: (i) convolving the tensor template with the block of the tomographic image to obtain a convolutional tensor field, the convolutional tensor field providing information about the correlation between the tensor template and the tomographic image at all locations and under all rotations within the block of the tomographic image; (ii) determining the optimal correlation value at each of the locations within the block of the tomographic image, thereby calculating a scalar map from the convolutional tensor field; (iii) identifying the location of an instance of the target particle within the block of the tomographic image by determining the local maximum of the optimal correlation value; and (iv) identifying the rotation of the target particle at each identified location among the identified locations using a global optimization technique.

[0106] 135. The method according to clause 134, wherein segmenting the tomographic image into blocks comprises segmenting the tomographic image into continuous units, and defining each block as comprising a unit plus an overlapping portion extending into the adjacent unit.

[0107] 136. The method according to Clause 135, wherein the overlap is a uniform distance, and optionally, the size of the overlapping portion extending into adjacent cells is less than the maximum size of the target particle and greater than half the test size of the target particle.

[0108] 137. The method according to Clause 136, the method comprising: for each block, retaining the locations identified within the cells of the block, and discarding any locations identified within the overlapping portions of the block.

[0109] 138. The method according to any one of clauses 134 to 137, the method further comprising, prior to calculating the normalized correlation value calculated for each identified instance of the target particle within the tomographic image, incorporating the rotation and position identified in each block of the tomographic image.

[0110] 139. A computer program comprising computer program instructions that, when executed by a computer processor, cause the computer processor to perform the method described in any one of claims 101 to 138; a computer-readable medium on which such computer program is stored; or a computer system comprising a computer processor and a computer memory in which such computer program is stored. Attached Figure Description

[0111] To facilitate a clearer understanding of the invention, reference will now be made to the accompanying drawings by way of example only, in which: Figure 1A The image shows a cryotomymogram of Chlamydomonas cells, and Figure 1B Demonstrates the use of detection and identification Figure 1A The surface model of the protein template shown; Figure 2 This is a block diagram of a method for identifying molecules in tomographic images; Figure 3 It is identification Figure 2 A schematic diagram of the method for analyzing molecules in tomographic images; Figure 4 This is a flowchart of the method for generating tensor templates; and Figure 5 This is a schematic diagram of an alternative method for identifying molecules in tomographic images. Detailed Implementation

[0112] Figure 1A A cryo-tomographic image 10 of Chlamydomonas is shown, which includes structures 15 such as the Golgi apparatus 151, endoplasmic reticulum 152, and nuclear membrane 153. The tomographic image 10 also contains numerous macromolecules 15, most of which are ribosomes 154 and 155, as well as proteasomes and CDC 48-ATPase 156. These molecules 15 have been detected by template matching. Figure 1B It shows the target such as Figure 1A The surface model of template 20 for a specific rotation of macromolecule 15 is shown. These templates 20 are derived from a computer-generated model of macromolecule 15 and can be generated as a representation of a specific molecule 15 with a specific rotation. These templates 20 can be used for detection and identification using the conventional template matching described above. Figure 1A An example of macromolecule 15 in the tomographic image 10.

[0113] Figures 2 to 4 A method 200 for identifying target molecules within a tomographic image 10 is shown, such as... Figure 1A The tomographic image within 10 Figure 1B At least some of the macromolecules in macromolecule 15. Figure 3 For simplicity and clarity, the tomographic image 10 is shown as a two-dimensional representation.

[0114] Method 200 begins by obtaining a tensor template 222 of the target molecule. The tensor template 222 utilizes a tensor field to describe the target molecule for all rotations. In conventional template matching, a template 20 of the target molecule is generated, in which the model of the target molecule is reduced to a finite but considerable number of snapshots of the target molecule at different rotations (e.g., tens of thousands to hundreds of thousands, depending on the target accuracy). The tensor template 222 described herein replaces the conventional template 20 with a tensor field of a finite tensor order embedded in all rotations of the target molecule.

[0115] In step 202, a model of the target molecule of interest is obtained, and in step 204, a tensor template 222 is calculated for the target molecule. The tensor template 222 uses a tensor field generated as follows to describe the shape of the target molecule over all rotations.

[0116] Many different rotations of the target molecule are sampled, generating a mathematical description for each rotation. These descriptions are then mathematically combined using integration to derive each tensor of tensor template 222, which describes the shape of the target molecule over all rotations. Tensors are computed for each voxel of the model of the target molecule, thus forming a tensor field.

[0117] Each tensor template 222 is a tensor field that uses tensors to embed all rotations of the target molecule. This tensor has as many independent components as an n-order and d-dimensional tensor, for example, 35 for a 4-order and 4-dimensional tensor (the number of components of a tensor depends on the order and dimension of the tensor used). In this embodiment, quaternions are used to describe each rotation of the target molecule. To describe the target molecule in each sampled selection during a sampled rotation, a tensor power is represented as q raised to the nth power of the tensor, where n is the order of the tensor and q is a vector representing the dimension d of the rotation. In this embodiment, quaternions are used as vectors. Each quaternion is transformed into a tensor power describing the corresponding rotation by multiplying the four components of each quaternion to the fourth power, resulting in a tensor power of dimension 4. This dimension is preserved in the calculations performed to determine each tensor of tensor template 222. This leaves room for the order of the tensor power, and thus for the selection of tensors. This is a trade-off: the higher the order of the tensor used, the higher the achievable accuracy. However, the computational cost also increases with the order. Using a 4th-order 4-dimensional tensor provides a good trade-off, and this results in a tensor template 222 with 35 linear independent components (i.e., 35 different combinations of the four components of a quaternion).

[0118] To combine the tensor powers derived from each individual quaternion (each quaternion describing a single rotation of the target molecule), an integral is performed over each voxel of the target molecule's model across all sampled rotations, thereby deriving a voxel-by-voxel representation of the target molecule across all rotations. For this purpose, for each voxel, each of the 35 components of the tensor power is combined with the corresponding rotational version of the target molecule's model to provide the corresponding component in the tensor of the tensor field. Thus, a tensor is created for each voxel, and each tensor has as many components as the nth-order tensor describing the rotation (35 in this example). Specifically, each component of each tensor is computed as the integral over all rotations of the product between the corresponding component of the tensor power generated by the nth power of the quaternion defining the rotation and the corresponding intensity value of the voxel of the target molecule's corresponding rotational model.

[0119] Figure 4 An example of a method 400 for calculating tensor template 222 is shown. In this example, the method begins at step 252, in which tensor fields are created. The tensor field provides a tensor for each voxel, and each tensor has 35 components, which are populated according to the following method steps comprising three nested loops.

[0120] When the first sampling rotation of the target molecule is selected, the first loop 402 is entered at step 254. This sampling rotation is represented by a quaternion, and at step 256, the quaternion used for the current rotation is converted into a corresponding tensor power. As mentioned above, this tensor power will have 35 components. At step 258, a template 20 of the target molecule with the current rotation is generated from the model of the target molecule.

[0121] Next, at step 260, a second loop 404 is entered, in which the tensor of the first voxel is selected. A third loop 404 sequentially fills the components of this tensor, and initially, at step 262, the first of the 35 components of the tensor is selected. Then, at step 264, the product of the corresponding component of the tensor power of the first sampling rotation and the intensity of the first voxel of the template 20 of the target molecule of the first sampling rotation is calculated. The value of this product is a scalar, and at step 266, this value is stored as the first component of the tensor of the first voxel selected at step 260.

[0122] With the first component of the tensor filled, the process is repeated for each of the other 34 components of the tensor by iterating through loop 406 multiple times via loop paths 268 and 269 and steps 262 to 266. That is, a new component of the tensor is selected at step 262, and at step 264 the product of the corresponding component of the tensor power for the first rotation and the first voxel of the model of the target molecule is determined. This product is then added to the current component of the tensor at step 266. This process is repeated for all 35 components of the tensor, such that each component of the tensor is filled with an initial scalar value. By doing so, the test at step 268 determines that all components have been filled, and thus the next tensor can be selected.

[0123] At step 260, the next tensor is selected, and the third loop 406 is repeated for the new tensor so that all components of the tensor are filled. As the method iterates through loop 406, each component of the new tensor is computed. Loop 404 ensures that this operation is repeated for each tensor until all voxel tensors have been filled with initial scalar values ​​corresponding to the product of the corresponding component of the tensor power and the intensity of the corresponding voxel of the first rotation of the target molecule. When the test at step 270 determines that all tensors have been processed once, we obtain a set of tensors where all tensors have values ​​added for all 35 components of each tensor. However, these values ​​correspond to the first sample rotation of the target molecule. This process is repeated iteratively for each sample rotation of the target molecule through loop 401.

[0124] Each time method 400 returns to step 254 via loop paths 272 and 273, a new rotation is selected. Then, at steps 256 and 258, a tensor power sum representation of the target molecule under the new rotation is generated. The method then continues iterating through each component of the tensor via loops 406 and 404, one tensor at a time. At each step 264, the product of the current component of the current tensor power and the intensity of the current voxel for the current rotation is calculated. Each product includes a scalar value that is added to the scalar value already stored for that component of the current tensor in subsequent step 266. As the method continues and iterates through successive rotations via loop 402, the steps performed for that rotation cause the scalar values ​​stored in the tensors to be updated to add the determined scalar value. Thus, when method 400 completes all iterations over the sampled rotations, each component of each tensor becomes the integral of all products of the corresponding component of the tensor power and the intensity of the corresponding voxel over all rotations. It is this set of tensors that constitutes the tensor field.

[0125] In summary: Tensor template 222 is a tensor field comprising a set of voxels of a model for the target molecule; each voxel has an associated tensor; each tensor comprises the same number of components as the tensor power, 35 components in this example; each component of each tensor is a scalar value; and each scalar value is equal to the integral of each product of the corresponding component of the tensor power with the intensity of the corresponding voxel in each rotational model of the target.

[0126] Before calculating the integral, the model of the target molecule can be processed by a symmetric positive semidefinite linear operator, which includes a spherical mask defining the local neighborhood of the model. This limits the number of voxels and provides additional advantages in subsequent steps of method 200 described below.

[0127] Although its computation is not instantaneous (ranging from minutes to hours, depending on the number of rotations during sampling and the number of voxels in the target molecule model), it only needs to be performed once per target molecule (although it may be updated periodically as the target molecule model becomes more refined). Once computed, the tensor template 222 can be stored in a library for future use. Therefore, as... Figure 2 As shown in the alternative path, method 200 can begin at step 206, where only the tensor template 222 of the target molecule 20 is retrieved from the library. Although in Figure 2 and Figure 3 Not depicted in the text, but before calculating tensor template 222, a model of the target molecule can be masked using a sphere that tightly but completely surrounds the target molecule. Masking the target molecule model helps generate the functional tensor template 222.

[0128] At step 208, the tensor template 222, obtained in any way, is convolved with the tomographic image 10 to generate a convolutional tensor field 226 containing information about the correlation between the tensor template 222 and each voxel of the tomographic image 10 throughout the entire rotation space used to generate the tensor template 222. Although in Figure 2 and Figure 3 Not depicted, but when using a spherical template mask to generate tensor template 222, this also focuses the correlation to the neighborhood of the target molecule size at each voxel in the approximate tomographic image 10. This convolutional tensor field 226 can then be used to determine the best image match for instances of the target molecule in the tomographic image 10 in terms of voxel position and rotation. To identify instances of the target molecule, the following multi-step process is performed.

[0129] In the first step 210, scalar values ​​for each voxel of the tomographic image 10 are calculated based on the convolutional tensor field 226. These correlation values ​​are calculated using a scalar value function that has critical values ​​indicating the optimal correlation location. The result of this calculation is a scalar map 228, which provides an estimate of the optimal correlation value for the tensor template 222 at each voxel. Figure 3 In this paper, for simplicity and clarity, scalar plot 228 is shown as a two-dimensional representation; however, it should be understood that scalar plot 228 is three-dimensional because it provides the best-case correlation value for the estimate at each voxel. The scalar values ​​can be computed using norms, for example... Figure 2 The Frobenius norm is shown in step 210. The Frobenius norm of a tensor is calculated as the square root of the sum of the absolute squares of its components, making its calculation very efficient.

[0130] Step 210 generates a scalar map with scalar values ​​for each voxel of the tomographic image. To determine the best candidate for an instance of the target molecule in tomographic image 10, the peak positions of the peaks within the scalar map 228 are determined. This is found by using any conventional peak-finding technique to determine the location of the local maximum of the relevant values ​​in the scalar map 228.

[0131] Next, at step 214, for each potential peak position found in step 212, a rotation for generating the optimal correlation value of the target molecule is determined. This can be performed using eigenvalue decomposition or by sampling the rotation space defined by quaternions, expanding the corresponding quaternions into n-order tensors, and computing their inner product with respect to the convolutional tensor field.

[0132] Figure 2 Step 214 illustrates an example of using eigenvalue decomposition. As described above, the convolutional tensor field 226 retains information about all rotations, and the optimal matching positions and corresponding rotations—that is, those that produce the highest correlation values—can be determined using eigenvectors and eigenvalues. This is performed by unfolding the tensor of the convolutional tensor field 226 into a matrix and determining the eigenvalues ​​and eigenvectors from the matrix, where each eigenvalue corresponds to the optimal correlation value. The eigenvectors and eigenvalues ​​for each identified position can be determined using a shift-symmetric higher-order power method. Since the symmetric higher-order power method is an iterative eigenvalue decomposition that only allows local searches, its accuracy is highly correlated with the starting point, which can be randomly chosen to ensure coverage of the entire rotation space. The actual optimal value is always found because the process is repeated multiple times, for example, using 1000 uniformly distributed quaternions as the starting point.

[0133] This technique can be used to identify optimal rotation and optimal peak location, which suggests that step 212 is unnecessary. However, if applied to all voxels in a typical tomographic image 10, determining eigenvalues ​​and eigenvectors is computationally expensive and therefore very time-consuming. To avoid this problem, a much faster technique using the Frobenius norm (which cannot identify rotation) is first used at steps 210 and 212 to determine the peak location. The Frobenius norm is correlated with its largest eigenvalue and can therefore be used as an approximation for finding the peak location. Thus, the Frobenius norm approximates the correlation value at each voxel, and the peak location corresponding to the target molecule location in tomographic image 10 is determined based on this correlation value. Therefore, steps 210 and 212 use the Frobenius norm and peak lookup to identify the peak location, which is then fed into a more refined technique using eigenvectors to determine the optimal rotation of the target molecule at step 214. This results in a much faster execution time.

[0134] Alternatively, to ensure that optimal rotations are not missed during eigenvector computation, the rotation space of the target molecule, represented as a quaternion, can be sampled. The inner product between the selected tensor and the quaternion of the convolutional tensor field is then computed, allowing rotations to be sorted by this value. Finally, a list of peak positions for each peak with its most promising rotation (highest inner product) is returned. Although sampling the rotation space is required, this process is faster than computing standard cross-correlation because computing the inner product between the tensor and the quaternion is less costly than computing regular locally normalized cross-correlation. Furthermore, this computation is suitable for further acceleration using high-performance computing techniques.

[0135] As described above, when the tensor template 222 is computed in step 204, a symmetric positive semidefinite linear operator can be used. This ensures that the inner product between the tensor and the unit quaternion extended to the nth order tensor reaches a global maximum for optimal template rotation. In step 214, this global maximum is found either by sampling the rotation space defined by the quaternion or by unfolding the tensors of the convolution tensor field 226 into matrices and finding their eigenvalues ​​and eigenvectors. Therefore, the use of a symmetric positive semidefinite linear operator ensures that the necessary information is provided to offer relevant values ​​reflecting the proximity of the tensor template 222 to any specific voxel in the tomographic image 10.

[0136] Step 214 can generate a list of potential peaks as peak locations and optimal rotations relative to relevant values ​​(or provide several candidates for optimal rotations, as described in the preceding paragraph), such as Figure 3As shown at point 215 in the figure. However, the result of using tensors is that the correlation values ​​are not normalized across the entire tomographic image 10. For example, due to variations in illumination at various locations within tomographic image 10, variations across tomographic image 10 are likely to exist. The lack of normalized correlation values ​​limits any attempt to rank potential peaks. Furthermore, when the list contains more than one rotation for each peak, normalized correlation values ​​are needed to allow for a meaningful ranking of the best correlation values, which in turn allows for the identification and selection of the optimal rotation for each peak.

[0137] Step 216 provides a normalized cross-correlation value for each potential peak by utilizing conventional techniques of localized cross-correlation (albeit in a very limited manner). Advantageously, in this method 200, calculations are only required for each potential peak location and optimal rotation value (or for each peak location and a limited list of optimal rotations, if provided). That is, the localized cross-correlation calculation is performed by comparing the conventional template 20 of the target molecule with the tomographic image 10: however, it is performed by comparing the template 20 of the target molecule at the location found in step 212 with only a single rotation (or several selected optimal rotations identified) of the target molecule corresponding to the optimal rotation found in step 214. This contrasts with prior art, where localized cross-correlation is calculated for each single rotation of the target molecule at each individual voxel (which is typically around 30,000 rotations). Therefore, computation time is significantly reduced and / or a greater number of optimal rotations can be used to improve accuracy.

[0138] Therefore, at step 216, a regular template 20 is generated from the model of the target molecule: template 20 provides a snapshot of the target molecule at the identified optimal rotation (or, if more than one rotation is provided, template 20 is provided for each optimal rotation). Then, a locally normalized cross-correlation value is calculated for the template at the identified location. In the case of more than one rotation provided in step 214, the cross-correlation in step 216 provides a normalized cross-correlation score for each rotation. This allows for a meaningful comparison of the cross-correlation scores for each optimal rotation: the rotation that achieves the highest cross-correlation score can then be selected and provided to the identified location.

[0139] At step 218, a list 219 of potential peak locations is provided. For each instance of the target molecule, list 219 contains the identified locations, optimal rotation, and normalized cross-correlation values ​​(although in some embodiments, optimal rotation and / or normalized cross-correlation values ​​may be omitted). List 219 may be provided in different ways and with or without user input. For example, list 219 may be sorted according to normalized cross-correlation values.

[0140] Additionally, as indicated at step 217, the user can set the number of potential peaks returned in list 219. Other parameters set by the user at step 217 may include an absolute peak threshold, a relative peak threshold, and a minimum distance between peaks.

[0141] The number of potential peaks returned in list 219, set by the user, varies depending on the intended use. As a first example, all instances of a specific target molecule may be desired, in which case all potential peaks can be returned. Alternatively, to limit false positives, a threshold-normalized cross-correlation value can be set such that only potential peaks meeting or exceeding this cross-correlation value are included in list 219 of potential peaks. As a second example, the user may want to refine the structure of the target molecule by averaging images of aligned instances of the target molecule. This can be done by aligning each instance of the target molecule with the optimal rotation value and then averaging the aligned images of the target molecule. In this example, using fewer potential peaks—those with the optimal normalized cross-correlation scores—may be better, as these instances will optimally correspond to the already determined target molecule structure (thus allowing further refinement of that determined structure).

[0142] Figure 5 An alternative method 500 using block processing to reduce memory requirements is illustrated. The tomographic image 10 is segmented into individually processed three-dimensional blocks 502 before the results are merged. Figure 5 In the illustrated embodiment, the tomographic image 10 is cubic (although shown as a square for simplicity) and is divided into eight consecutive cubic blocks. As will be understood, other arrangements may be used depending on the size and shape of the tomographic image 10 and the available memory. Figure 5 A block from block 502 is shown, again as a two-dimensional square within a two-dimensional visualization of the diagram. Method 500 corresponds to method 200, which has already been described, although it is implemented on a block-by-block basis. Therefore, only the differences will be explained now.

[0143] Step 208 involves convolving tensor template 222 with a portion of tomographic image 10 to identify peaks within a block of block 502, and thus identify molecules. To account for the fact that one or more target molecules can span the boundary between adjacent blocks 502, tensor template 222 is convolved with a portion 505 of tomographic image 10, which includes block 502 and an overlapping region (halo) 504 surrounding the boundary of block 502 to overlap with adjacent blocks 502. Thus, the overlapping region 504 extends into the adjacent blocks 502. The width of the overlapping region 504 is chosen to be the maximum size (i.e., height, width, depth) around the target molecule, for example, between half and the full maximum size, as this ensures that any instance of the target molecule is enclosed within a block 502 and its overlapping region 504.

[0144] Section 505 is then analyzed in the same manner as described for method 200. The convolution at step 208 produces a convolutional tensor field 226, which is processed at step 210, for example, using the Frobenius norm, to produce a scalar graph 228. Step 212 sees the peaks identified based on local maxima among the relevant values ​​present in the scalar graph 228. Only peaks within block 502 are identified, or if peaks are identified in both block 502 and overlapping region 504, only peaks in block 502 are retained, and peaks in overlapping region 504 are discarded (because these peaks will be identified when processing block 502 containing overlapping region 504). At step 214, optimal rotations are determined, for example, using the eigenvectors and eigenvalues ​​found for each peak position in the peak list to determine the optimal rotation at each peak position, or alternatively using the inner product as described above to find several rotations for each peak.

[0145] Steps 208, 210, 212, and 214 are repeated for all blocks 502. At step 506, the lists of peaks and rotations found for each block 501 are merged into a single list. Then, at step 216, normalized cross-correlation values ​​are determined using regular localized cross-correlation, and these values ​​are calculated only for each peak position and optimal rotation pair (or for each rotation at each position if more than a single rotation is identified in step 214, in which case the optimal rotation is also identified from the normalized cross-correlation values). A refined list of peak positions, rotations, and normalized cross-correlation scores is provided in step 218, which can be guided by user input of the desired number of peaks at step 215.

[0146] Those skilled in the art will understand that changes can be made to the above embodiments in many different ways without departing from the scope of the invention as defined by the appended claims.

[0147] While the above embodiments relate to the use of the invention in Cryo-ET, the invention has many other applications. For example, the invention can be applied to both 2D and 3D images. The invention can be used for object detection in image processing, object counting in metrology, benchmark tracking for drift correction, benchmark detection for tomographic image reconstruction, defect detection in semiconductors, etc. Even within Cryo-ET, the invention can be used for both protein structure discovery and functional cell behavior studies.

[0148] The following is a mathematical explanation of the concepts and proofs used to demonstrate the effectiveness of the above methods.

[0149] 1 Introduction In image processing, particularly in pattern recognition, a classic problem is determining whether a large image contains copies of a smaller image known as a “template”—and the number, location, and orientation of these copies. The resulting algorithms are generally referred to as template matching algorithms [1, 2, 3]. The most classic solution is based on cross-correlation, although other methods exist based on, for example, metaheuristic algorithms [4] or deep learning [1, 5, 6]. In this paper, we present the mathematical foundation of a cross-correlation-based template matching algorithm (TM in all algorithms below), and we introduce a novel, fast algorithm that solves the problem using tensors.

[0150] Compared to machine learning-based algorithms, the main advantage of TM is that TM is a white-box model, which is directly applicable when you only have a template and a large image (no kind of training required, which can be a very difficult task in some applications), and it can locate rotations with arbitrary precision (current deep learning-based algorithms for template matching in 3D images cannot accurately estimate rotations [5]).

[0151] On the other hand, the main drawback of TM is its computational cost. The basic idea of ​​TM is to compute the inner product between the (rotated) template and the (translated) image and normalize the result. For each rotation, these computations are performed in the Fourier domain to efficiently solve the translation [7, 8, 9]. However, this process must be repeated for each rotation to be investigated, so the resulting complexity is dependent on the rotation being processed. The computational cost of this process can become limiting for 3D images because... The rotation space SO(3) is a (compact) manifold with dimension 3. In applications such as cryo-electron microscopy, more than ten thousand rotations are required to achieve angular accuracy of a few degrees.

[0152] We propose an algorithm called Tensor Template Matching (TTM), which integrates information about all rotations relative to the template into a unique symmetric tensor. In other words, the tensor template contains information about all rotations of the template in a unique object, allowing us to find the position and rotation of the template instance in any tomographic image with only a few correlations to the linearly independent components of the tensor. Each tensor template is computed only once, and once generated, it is capable of processing any image.

[0153] 2 Classic Template Matching Let's introduce some symbols; a d-dimensional image is simply an L-dimensional image. 2 ( The elements of a Hilbert space are those having an inner product. Naturally, the inner product is used to compare two images ƒ and g of the same size. Specifically, we can use... , where θ is the angle formed by f and g. Specifically, for some positive constant α, if Then ƒ = αg.

[0154] Template matching is often used to investigate whether instances of a "small" image t (the template) exist in a larger image ƒ, for example, to find instances of a specific macromolecule in a cryogenic electron tomography image (a 3D volumetric image). The size of the image is related to the set of points in the image that do not vanish (the image's support set). That is, when the set It is small (e.g., a small ball). (a subset of ƒ), where t represents "small". Let's assume ƒ and t have completely different sizes, so our interest is to compare t (the template, the small image) with only a subset of ƒ. In this case, we need to introduce some special operator S:L 2 ( ) → L 2 ( These operators fix our attention on a subset of the domain of ƒ. An interesting example of such an operator is... Where U: Let represent the unit step function of Heaviside, and r > 0. If the support set of template t is With 0 ∈ For a sphere with radius r centered at its center, the normalized inner product is... This reflects the relationship between t and f. The similarity between the restricted parts. Furthermore, if we introduce the translation operator τ x : L 2 ( ) → L 2 ( ), τ x (f)(z) = f(z + x), and calculate The result reflects the relationship between t and f. The similarity between the restricted parts. Of course, it is possible that f contains a rotated version of t, therefore, rotation is also necessary for a complete discussion of the problem. Therefore, given R∈SO(d), we define operator O R : L 2 ( ) → L 2 ( ), O R (t)(z) = t(R z ), and for t∈L 2 ( We define a rotational version of t.

[0155] (2) Normalized inner product Reflects t R With f in The similarity of the restricted parts. It's important to note that ||t||2 = ||t R ||2.

[0156] The operator S defined by (1) r It possesses some special properties. Specifically, it is symmetric, semidefinite positive, and commutative with rotations. To recap, given... A (real) Hilbert space (which we will also denote below as x·y if this simplifies the calculation) The operator S : X → X is called: • Symmetry (also known as self-adjoint), if For all f, g∈X • Semi-definite positive, if For all f∈X, • Positive definite, if it is semi-positive definite and This means f = 0.

[0157] If S : X → X is a symmetric semidefinite positive operator (SSP, in all below), then X becomes a pre-Hilbert space with the following inner product and seminorm: (3) and seminorm For example, it is observed that if S is given by (1), then ||ƒ||s = 0 means They are almost everywhere.

[0158] Theorem 2.1 Set X = L 2 ( Let S : L 2 ( ) → L 2 ( ) is an SSP operator, and consider the inner product given by (3). Then (a) For all f, g∈L 2 ( ), .

[0159] Furthermore, if f, g∈L 2 ( ), ||g|| S If ≠ 0, then the following statements are equivalent: (b) .

[0160] (c) .

[0161] Proof. Since S is an SSP, therefore for all α∈ , (5) Therefore, ||g|| S When ≠ 0, the only condition for the above (with respect to α) quadratic polynomial to be nonnegative everywhere is: This expression is equivalent to (6) On the other hand, if ||g|| S = 0, then the only way to satisfy formula (5) is In this case, (6) also holds true. This proves (α).

[0162] Now let us explain As long as ||g|| S ≠0. In fact, (c) is equivalent to If and only if , therefore .

[0163] Note that if ƒ, t∈L 2 ( ) are two images, α > 0, and we take S = S given by (1). r ,but This means that ƒ is within a unit sphere centered at x and t R Matching. In fact, there are many ways to define an operator S with the following properties: This means that f = g is in the neighborhood of 0, therefore This means that f matches the rotated version of t in the neighborhood of 0. Although arbitrary SSP operators may not enjoy this property, they allow for the creation of general ways to process such operators.

[0164] Therefore, in all the following, we assume S : L 2 ( ) → L 2 ( ) is an SSP operator and , where 1(x) = 1 is a constant image.

[0165] Note that 1 is not L. 2 ( However, this can be managed in several ways. In fact, although we use L... 2 ( Let f be used to represent the space of a d-dimensional image, but in practice we only consider images f with a compact support set K. Then, when we compute... We mean Furthermore, since our interest lies in the operator S that takes zero value on a function that is zero outside a certain neighborhood D of 0, we use... express .

[0166] Rotation and combination of operators play a crucial role in the embodiments of the invention described herein. Therefore, it is natural to ask how the composition of rotations affects the image. This is actually a simple calculation: therefore as well as Given an image ƒ, we consider its projection onto an image space that is S-orthogonal to the constant image 1. Note 2.2 These projections are important for studying the properties of an image under translation and rotation relative to constant brightness. Note that in an image f, the projections are of the form f + α1, α∈ There is no "real" difference between the images. When we change the constant α, what we observe is a uniform change in density or brightness in image f, rather than the emergence of new structures or forms. Therefore, f and its projection P S (ƒ) essentially represents the exact same image, because for some α∈ f = P S (ƒ) + α1.

[0167] Given two images ƒ, t, we have For some constants α and β, f = P S (f) + α1 and t = P S (t) + β1.

[0168] therefore Because P S (f), P S (t) ⊥ S 1. Therefore, if x∈ And if R∈SO(d), then there are two constants ρ = ρ(x) and δ = δ(R) such that .

[0169] Suppose S is interchanged with the rotation, and x∈ Fixed. Then, for each R∈SO(d), we have because .also, (Only take) And use det R = 1) (Because O) R ◦S = S◦O R ) (Because O) R (1) = 1) .

[0170] therefore And together with (10) means .

[0171] therefore t R = P S (t) R +β1 It is t R The S-orthogonal decomposition means that in t R In the S-orthogonal decomposition, the constant β multiplied by 1 does not depend on R, and .

[0172] In particular, for each x∈ question: • Maximize on rotation R .

[0173] • Maximize on rotation R .

[0174] • Maximize on rotation R .

[0175] They are equivalent.

[0176] Let us define it as follows: (11) Lemma 2.3 If S is an SSP operator with rotational commutation, then the parameter δ appearing in the orthogonal decomposition of S is... It does not depend on R. Therefore, given x∈ ,question • Maximize on rotation R .

[0177] • Maximize on rotation R .

[0178] They are equivalent.

[0179] Proof. We know... Therefore, we only need to prove It does not depend on R. In fact, (Make the change in variable z = R) y ) (Because S and) exchange) (because ) We can now state and prove the following: Theorem 2.4 (Classical Template Matching) Let S be the SSP operator, which is commuted with rotations and let x∈ It is fixed. Then, the following is the equivalence problem: (a) Maximize on rotation R .

[0180] (b) Maximize on rotation R .

[0181] (c) Maximize on rotation R .

[0182] (d) Maximize on rotation R .

[0183] Furthermore, if ||t|| S > 0 and S also has the following property ||f|| S = 0 means that in 0∈ In a certain neighborhood D, f |D = 0, this neighborhood contains all rotated templates t Q If f is the support set of (Q∈So(d)), then f and t are equal if any of the following statements are true. R A match exists at x: Finally, when we replace f with αf + β and t with δf + γ (where α, β, δ, γ ∈ ... When α,δ≠ 0, the normalized correlation values ​​described in (a*), (b*), (c*) and (d*) remain unchanged.

[0184] Proof. Equation and It has already been shown.

[0185] The following identity proof : (Only take Rz=y and use det R=1) (Because S and) exchange) The other claims are direct consequences of Theorem 2.1.

[0186] In all the following content, we assume that S is an SSP operator with rotational commutation, and that t is in t⊥S 1 and ||t|| S It is normalized in the sense that = 1. Then, for all rotations R, P S (t) R )=t R And ||t R || S =1. Therefore, as well as If and only if t in f and x R Its maximum value (= 1) is reached only when a perfect match exists. Furthermore, if we define... And consider the function ƒ, g∈L 2 ( ), which is defined as but (13) In general, a perfect match is never achieved. This is because the desired image represented by the template t is typically supported on a strict subset Ω of the domain D, where the operator S is capable of distinguishing functions. Therefore, the image f can well contain the image represented by t. R A copy of the image represented, but in t R In the supported neighborhood, f will contain t R Some information is missing in it. Furthermore, f is often corrupted by noise and distortion. This means that in Theorem 2.4... The normalized correlation described in the item will never equal 1. Therefore, a threshold should be introduced to determine whether a match has been generated (or not).

[0187] To find the rotation that maximizes c(x, R), cross-correlation is needed. The computation should be performed for a large number of rotations R, which makes classical matching an inefficient method for template matching. In fact, for d=3, the size of the rotation set R used to sample SO(3) to guarantee reliable results is 10. 4 Rotation with 10 7 Changes between rotations.

[0188] Due to numerical reasons, high frequencies can be altered during rotation transformation. Therefore, in practice, we do not apply the operator S to the original images f and t, but rather to their filtered versions that eliminate these high frequencies. Specifically, an isotropic (i.e., rotation-invariant) low-pass filter h is applied to both images, and subsequently, a template matching algorithm is applied to the resulting image. The idea behind this is that if a match exists between ƒ and t, then a template matching algorithm will be applied. and Similarly, the operator S is generated by applying a rotationally symmetric mask m(x) = ρ(||x||) to a given image. Therefore, we use... Replace f and use Replace t. Then, we use the SSP operator. Apply classical (or tensor) matching algorithms to this pair of images. Typically, the mask m is equal to 1 within a certain radius around 0, and equal to 0 outside a larger radius. Between these radii, the mask takes values ​​between 0 and 1. Under these constraints, it is clear that the operator S is an SSP and is interchangeable with rotations. Furthermore, if We have f |D = 0, where D is a sphere of positive radius centered at 0. Let's calculate the inner product. (Because each filter is translation invariant, and h is isotropic) (use (This is derived from the isotropy of h) in Furthermore, we use · to denote the standard product of real functions. This means we will have only considered the operator. The same effect of the associated template matching algorithm is applied to images f and t. Furthermore, the following holds true: Lemma 2.5 Let SS : L 2 ( ) → L 2 ( The following formula is given. (14) Where h defines an isotropic filter, and m defines a rotationally symmetric mask as described above. Then S is an SSP.

[0189] Proof. For the proof, we use the following (well-known) formula: for functions a, b, c, ∈ L 2 ( We have Make as well as .

[0190] Now let's consider the product. : (because (and m≥0). This proves that S is semidefinite positive. Let's show the symmetry: (because And · is the standard product of functions) Throughout the following content, we assume that the SSP operator S is of the form (14), where h and m satisfy the assumptions of Lemma 2.5. Therefore, the template matching algorithm is applied in conjunction with this operator, and fast computation of c(x, R) is required.

[0191] Direct calculation results in: (because ) (Because h is isotropic and m is rotationally symmetric). Furthermore, Now, using the definition of S (and applying h * 1 = 1), we can simplify the calculation as follows: It should be noted that the FFT algorithm can be used to calculate the inner product appearing at the end of the above formula, which helps to speed up the algorithm. In fact, if f and g are two images... , making (15) In addition, the following identity also holds: as well as .

[0192] therefore, An algorithm specifically designed for classical template matching was developed using the above formula.

[0193] In the embodiments of the invention described herein, an important tool we will use is the set of quaternions. In particular, we will use a unit quaternion (which can be equivalent to a unit 3-sphere). The rotation is parameterized by the following equation (see, for example, [10, 11]): •If x∈ If it has a norm of 1, then .

[0194] • Given x∈ x = a + bi + cj + dk, we equate x with a pair (a, v), where a ∈ and v = (b, c, d) ∈ And a is called the real part of x, denoted as a = Re(x).

[0195] Then, if x = (a, v), y = (b, w) ∈ Then we have Therefore, if x, y∈ If it has a norm of 1, then We conclude this section with a well-known result concerning the combination of SSP operators, which will be used to prove the main theorem described here: Lemma 2.6 If T, S : L 2 ( ) → L 2 ( If ) is semidefinite positive symmetric and commutative, then TS is also semidefinite positive symmetric.

[0196] Proof. If S and T are SSPs, then they have (unique) positive square roots. They are also symmetric operators. Furthermore, composite operators... It is also symmetrical. Since T and S are interchanged, their positive square roots are interchanged. Therefore... Furthermore, the square of any symmetric operator is positive. Algorithm 1 with rotational classical template matching

[0197] Image loading and Fourier transform Load the data into f. Load the mask and preprocess the image. Load the mask into m Load the template and preprocess it. Load the template into t. Calculation output 3-Tensor Template Matching This section introduces the Tensor Template Matching (TTM) algorithm. The goal is to efficiently handle both translation and rotation simultaneously. First, we introduce some background necessary to understand further mathematical developments. Second, we present the main theorem of TTM, which allows us to determine the optimal rotation of template t on each match in image f by computing some tensors without sampling SO(3). Finally, we explain how the matching position (template translation) can be determined directly from the computed tensors.

[0198] 3.1 Tensor Background A tensor A∈T of order n and dimension d n ( ) is a form of An array of elements All are real numbers. If for each permutation... (The permutation set of {1, ..., n}) has If , then the tensor A is said to be symmetric. We use S n ( Let f(n) represent the set of symmetric tensors of order n and dimension d. An important example of an n-order symmetric tensor is the vector υ = (υ1, ..., υ2). d )∈ The so-called nth power of the tensor is defined as follows: (17) As we all know, T n ( ) and S n ( ) is a real vector space with natural operations (pointwise summation and multiplication by a scalar), and for all n, d ≥ 1, dim T n ( ) = d n dim For example, T 4 ( ) = 4 4 =256, dim Furthermore, each symmetric tensor is a finite sum of tensor powers, which allows us to introduce the concept of the (symmetric) rank of a symmetric tensor as the minimum number of tensor powers used to represent a tensor in terms of its sum

[12] .

[0199] Mapping : T n ( ) × T n ( ) → It is given by the following formula Define the inner product. It is also often represented as... Furthermore, using this notation, if x, y∈ If it is a d-dimensional vector, then the direct application of the polynomial theorem shows... (19) Furthermore, if A∈S n ( ), and x = (x1, ⋯, x d )∈ Then we can also consider the inner product. It can be viewed as a homogeneous polynomial of degree n with d variables, which proves the use of notation. The rationality of this. Furthermore, if k < n, Ax k ∈S n–k ( If ) represents a symmetric tensor, then its components are: .(twenty one) Specifically, Ax n-1 ∈S 1 ( A vector whose i-th component is .(twenty two) In fact, if ,but in Representation function The gradient.

[0200] Note that the vector x can be selected from the definition above. This proves the following definition (see

[13] ): given A∈S n ( ), B∈S m ( If Au n–1 = λBu m–1 And Bu m = 1, then λ∈ Let u be an eigenvector of A, and u ∈ It is its corresponding eigenvector (equivalently, (λ, u) is called a B-feature pair of A).

[0201] Using the gradient, we can transform the equation Au n–1 = λBu m–1 Rewritten as .

[0202] Therefore, u is an eigenvector of A in B if and only if it is the critical point of the following optimization problem: Two particularly important examples are H-eigenvectors and Z-eigenvectors The optimization problem associated with finding the Z-eigenvectors of a given symmetric tensor is particularly important to us because our proposed tensor matching algorithm is reduced to one of these problems at each location, and fortunately, there are zero-order iterative algorithms that can approximate the solution (25) (see, for example, [14, 15]). These algorithms have linear convergence speeds. In Section 3.4, we show the heuristics that can be used to select possible matching locations, thus making the solution (25) necessary.

[0203] 3.2 Definition of Tensor Template In all of the following, S : L 2 ( ) → L 2 ( Let ) denote the SSP operator, which is related to rotation, and assume that the template t∈L 2 ( ) via t ⊥ S 1 and ||t|| S = 1 for normalization. In Section 2, we prove that for every x∈ , It reaches its maximum value (and that value equals 1) on rotation R if and only if there is a match between f and t at (x,R) (i.e., at τ). x (f) and t R (Matching between them). Let's define the symmetric tensor C using the following formula. n (x)∈S n ( ), where d' is the number of parameters used to describe the rotation SO(d) (in particular, for d = 3, we have d' = 4).

[0204] This means For all 1 ≤ i1, …, i n ≤ d′, therefore in It is a tensor template (or tensor pointer).

[0205] It is important to observe that T(z) is calculated only once and includes a component with a reduced number of components, because dim (In particular, dim) In fact, this is one of the main reasons why our introduced tensor template matching algorithm is fast. Another reason is that we define f and t... R The rotation R matched at x is a symmetric tensor C. n The Z-eigenvectors of (x) are very significant because the power method used to solve the corresponding optimization problem in [15, 14] does not require thousands or even millions of rotations like the classic matching algorithm, but only a small number of rotations: one for each iteration.

[0206] 3.3 Find the correct rotation Let us state the main results of this paper: Theorem 3.1 Let f,t∈L 2 ( ), x∈ and Given. If f and t R If there is a match at x, then it is defined in On rotation, a function parameterized by a unit quaternion It reaches its global maximum value (and that value is equal to 1) at Q = R.

[0207] Proof. We have proven the result as a corollary of Theorem 2.4. Therefore, our main objective is to express φ(Q) as a scalar product. (where S′ is some SSP operator), this will guarantee that if f and t R If there is a match at x, then φ(Q) reaches its global maximum at Q = R (and that value is equal to 1).

[0208] Let's calculate φ(Q): (By the definition of c(x, R)). Therefore, dividing by w(x), we get: (in ) Where K(R) = (Re(R)) n ,and Represents functions a, b∈L2 Convolution of (SO(3)).

[0209] in other words, Here, It must be interpreted as defined in The value is above. The function (actually, it is L) 2 ( (elements of) Where S(t)(z) : SO(3) → It is given by the following formula: .

[0210] In fact, generally speaking, each element t∈L 2 ( For each z∈ Both can be achieved by setting t(z)(R) = t R (z) is interpreted as L 2 An element of (SO(3)). Therefore, It is L 2 ( () elements.

[0211] (30) Therefore, we can conclude that... (In the last equation, let R = QP and use |Q| = 1) Where I d It indicates an identical rotation.

[0212] Now let's use S2 : L 2 ( ) → L 2 ( () represents the operator given by the following formula: And let S′ = S2◦S, then (Because S,S2 are interchanged with rotation) Therefore, the proof ends once we prove that S' is an SSP operator, and for this we need to use Lemma 2.6. In fact, S′ = S2◦S is a composition of operators, by the assumption that S is positive semi-definite, and S commutative with S2 because S commutates with rotations, and S2 is defined by convolution over SO(3). Therefore, Lemma 2.6 implies that S' is an SSP as long as S2 is an SSP.

[0213] To prove that S2 is symmetric positive semi-definite, we use the fact that convolution on SO(3) is interpreted as unit sphere S 3 Properties of convolution on a hypersphere. Recall that if S... d–1 = {x∈ : x·x t = 1} means If the (unit) sphere is , then SO(d) can act transitively on S. d–1 Above (This means that, given z1, z2∈S) d–1 There exists a rotation R∈SO(d) such that R(z1) = z2), which makes S d–1 It becomes a homogeneous space and allows the introduction of definitions in S d–1 The convolution of the function on is as follows: Where η∈S d–1 It is the north pole of the sphere, and f, g∈L 2 (S d–1 Therefore, if we use the elements of SO(3) consisting of quaternions with a norm of 1 (which are compatible with spheres in four-dimensional space), If the elements are equivalent, then assume S is parameterized. 3 The North Pole is exactly caused by an identity rotation I d Given that f, g∈L 2 The convolution of (SO(3)) can be interpreted as S 3 Hyperspherical convolution on: Now, as is widely known, L 2 (S 3 ) is Hilbert space, and the so-called hyperspherical harmonics This forms an orthogonal basis for the space. Therefore, for each function f∈L 2 (S 3 Allow Fourier expansion Furthermore, it has been proven in

[16] that if f, g∈L 2 (S 3 )and ,but Therefore, given a template t, for each z∈ The mapping t(z)(R) = t R (z) belongs to L 2 (S 3 (Here, the rotation R is parameterized as a unit quaternion, therefore R∈S) 3 )and as well as We need the following lemma, the proof of which is included in Section 4: Lemma 3.2 For all .

[0214] So This means the convolution with K, i.e., the operator C. K :L 2 (SO(3)) → L 2 (SO(3)), It is semidefinite positive. Furthermore, it is well known that the operator is symmetric (and we will use both of these things in the calculations below).

[0215] To prove that S2 is an SSP, we introduce the operator , defined as L(t)(R) = t R and operators Defined as .

[0216] So Where L(f)(w)(R)≔L(f)(R)(w)=f R (w)=f(R -1 w) and a(w)(R) ≔ a(R)(w). On the other hand, Therefore, if V = ʃ SO(3) Let dR be the volume of SO(3), then Where we use ,and Make (set up = Q –1 R –1 , making Q –1 = R) Therefore, we can conclude that... Therefore, S2(t) is semidefinite positive.

[0217] Furthermore, the same type of calculation shows This proves that S2 is symmetric. This concludes the proof of the theorem.

[0218] Note that Theorem 3.1 connects the following two problems: finding at a given position x such that f is equal to t. R The rotation R matched at x; by solving (25) (where A = C) n (x)∈S n ( (where n is an even number) to find the dominant Z-eigenvalue-eigenvector pairs.

[0219] 3.4 Find the correct location In the previous subsection, we showed how the correlation function based on tensors can be viewed using a slightly different degenerate inner product (based on S' instead of S from the proof of Theorem 3.1). Unfortunately, this is not considered in the normalization: not only for the template t, but more importantly, for the image f. So what are the implications? First, it is observed that the operation missing from S in the normalization is actually a convolution, and it should at most scale the constant components of the image. Therefore, if the image is orthogonal to 1 by S, it will also be orthogonal to 1 by S'. However, the norm will be affected.

[0220] For template t, this means the normalization deviates from a certain factor, but this factor is the same everywhere, so it doesn't affect the ratio between responses at different locations. For image f, the effect is less pronounced because... and Differences will arise in a non-uniform manner.

[0221] So, when will this shortcoming lead to misidentification? For this to happen, the normalization factor used at the mismatched locations would have to be much higher than the "correct" normalization factor, and / or the normalization factor would have to be (far) too low at the matched locations. Since the difference between S and S' is essentially a smoothing operation, and the normalization factor is the inverse of the norm of the projected image, the image must be (very) smooth at the mismatched locations while exhibiting a large amount of high-frequency energy around the matched locations. This is not impossible, but it is at least unusual in the context of typical applications such as the analysis of electron microscope images.

[0222] While we could find the spatial location of the peak by running an algorithm to find the dominant Z-eigenvalue-eigenvector pairs for each voxel, doing so using current high-order tensor decomposition algorithms is quite expensive. However, the Frobenius norm of a tensor is related to its spectral norm, and in practice, it has proven to be an excellent surrogate for finding the spatial location of the peak. In fact, we know that C... n (x)∈S n ( ) = S n ( Now, if we use ||T|| σ Let ||T|| denote the spectral norm of tensor T. F Let T represent its Frobenius norm. It is well known that the largest singular value of T is equal to its spectral norm, and... (See, for example, [17, 18]).

[0223] In fact, ||T|| σ With ||T|| F The connection between them is stronger than just this inequality. It is well known that each tensor is a finite sum of rank-1 tensors (in fact, if the tensors are symmetric, then rank-1 tensors can also be chosen to be symmetric)

[12] . Furthermore, if W1 is a rank-1 tensor, then it satisfies Then (see, for example

[19] ) as well as therefore, .

[0224] Therefore, if E1(T) is retained, then ||T|| σ (Correspondingly, ||T|| FThe increase in size is transformed into ||T|| F (Correspondingly, ||T|| σ ).

[0225] Furthermore, in 1938, it was proved in

[20] that for any symmetric tensor T, Therefore, large ||T|| F This implies a large tensor spectral norm, while C n The spectral norm of (x) is closely related to the optimization problem solved in Theorem 3.1, which proves the use of C n It is possible to use the Frobenius norm of (x) as a parameter to select the position x where a match may exist.

[0226] For each location x identified as a potential peak, the SS-HOPM algorithm (for the precise definition and implementation of the algorithm, see [[14, 15, 19]) is used to find the exact dominant Z-eigenvalue and its associated Z-eigenvector, which gives a matching rotational R candidate at x.

[0227] 3.5 Template Points Correctly integrating the “pointer” over all rotations is a key part of the process. In this work, this is achieved by summing over a large set of rotation samples. In conventional methods, the exact distribution of rotations is not very critical, as we only retain the maximum response. However, when integrating with samples, they must be distributed as uniformly as possible. Unfortunately, it is well known that it is impossible to find an arbitrarily large set of points that is perfectly uniform on a (hyper)sphere (except a circle). Here, a two-step approach is used: first, points representing the axes of rotation are generated on the sphere using the method proposed in

[21] (with some fine-tuning), and then these axes are combined with (1D) rotations that are distributed in a manner that is still uniformly distributed on the 3D sphere.

[0228] 4. Lemma Proof 3.2 Let's begin by reviewing the relationship with the sphere S. 3 The Fourier expansion formulas related to the hyperspherical harmonics are as follows. The parameterization of the sphere we consider is as follows: Where (a, b, c, d) ∈ S 3 Equivalent to the unit quaternion Q = a + bi + cj + dk, this quaternion represents three-dimensional Euclidean space. Rotation. Volume element (for S) 3 The integral over, and therefore also for SO(3), is given by the following equation. dV = sin2 θ sin ϕdθdϕdφ.

[0229] Then each function f ∈ L 2 (S3) can be decomposed into in L represents the formation of hyperspherical harmonics 2 (S 3 An orthogonal basis of ) and Let f be the Fourier coefficient of f in the basis. We want to prove that for all K( ,(0, 0) ≥ 0. Now, and K(Q) = (Rc(Q)) n = a n = (cos θ) n ,and in It is a positive constant, and express The second-order Geiger-Bauer polynomial, as its expansion: The first in The Taylor coefficient appears. It is well known that... ( (a second-order Chebyshev polynomial), and .therefore To estimate the above integral, we need to use several trigonometric formulas, as well as the assumption that n is even. Specifically, n being even means that n / 2 is an integer, and (cos(θ)) n = (cos(π – θ)) n Furthermore, for odd numbers, we have This makes for ∈ is +1, the integral equals 0.

[0230] Assumption ∈ is .So (For the last line, simply assume k = n / 2 – s). Furthermore, it is well known that... Therefore, by directly substituting the definition From the formula, we obtain We can now utilize to assert that The parity of and n implies that all factors multiplying the variable θ that appear inside the cosine function are even. This makes the corresponding integral (over [0, π]) equal to 0, except in the case where the factor itself is 0. In this case, cos(0) = 1 implies that only the cosine functions that appear with a negative sign in front of them in the formula can contribute negatively to the integral. Now clearly, +2 > 0 always holds (since ≥ 0), while + 2 + n – 2k = 0 implies 2k = n + + 2 > n, and thus, k > n / 2, which is impossible since the summation ranges from k = 1 to k = n / 2. This means that the term –cos(θ( + 2 +n – 2k)) never contributes negatively to the sum. On the other hand, if the cosine function with factor + 2 + 2k – n contributes, which means that + 2 + 2k – n = 0, then k = (n – – 2) / 2 < n / 2. In particular, taking k* = k +1, we have 1 ≤ k* ≤ n / 2, so cos(θ( + 2k* – n)) = cos(θ( + 2k + 2 – n)) =cos(0) = 1, and the corresponding term effectively appears in the summation. In particular, adding these two terms, we get since n is even and k < n / 2. This concludes the proof of Lemma 3.2 5. Conclusion We have revealed the mathematical principles of classical template matching with rotation. In addition, an alternative to the classical algorithm, called tensor template matching (or TTM), has been shown. TTM integrates the information related to all rotated versions of the template t into a unique symmetric tensor template T, which is computed only once per template. The main theorem (Theorem 3.1) shows finding the rotated version t of the template t with the image f at a given position x RAn exact match between them is equivalent to finding a match for a given tensor C. n The best rank-1 approximation of (x) (in the Frobenius norm). The resulting algorithm has reduced computational complexity compared to the classical algorithm. TTM only utilizes some correlation with the linear independent components of T to find the position and rotation of the template instance in any tomographic image. In particular, cryo-electron tomography (3D images) for macromolecule detection requires 7112, 45123 and 553680 rotations to achieve 13°, 7° and 3° accuracy, respectively

[22] . Therefore, considering the 4th order tensor (35 linear independent components), our method has a potential speedup of 203x, 1239x and 184560x relative to TM in these cases.

[0231] References [1] R. Brunelli. Template Matching Techniques in Computer Vision: Theory and Practice. John Wiley Sons, 2009.

[0232] [2] David Forsyth and Jean Ponce. Computer Vision: A Modern Approach. Prentice Hall, 2002.

[0233] [3] R. Gonzalez and R. Woods. Digital image processing, 4th Global Edition. Pearson Education, 2017.

[0234] [4] G. Corona, O. Makiel-Castillo, J. Morales-Castaned, A. Gonzalez and E. Cuevas. A new method to solve rotated template matching using meta-heuristicalgorithms and the structural similarity index. Mathematics and Computers in Simulation (MATCOM), Vol. 206 (C): 130-146, 2023.

[0235] [5] L. Lamm, RD. Righetto, Wietrzynski, Pöge, A. Martinez-Sanchez, T. Peng, and B. Engel. Membrain: A deep learning-aided pipeline for detection of membrane proteins in cryo-electron tomograms. Computer Methods and Programs in Biomedicine, Vol. 224: 106990, 2022.

[0236] [6] Moebel E., A. Martinz-Sanchez, L. Lamb, R. D. Gighetto, Wictrzynski, S. Albert, D. Lariviere, E. Fourmentin, S. Pfeffer, J. Ortiz, Baumcister, T. Peng, B. Engel, and C. Kervrann. Deep learning improves macro-molecule identification in 3d cellular cryo-electron tomograms. Nature Methods, Vol. 18 (No. 11): 1386-1394, 2021.

[0237] [7] JP Lewis. Fast template matching. In Visioninterface, Vol. 95, pp. 15-19. Quebec City, QC, Canada, 1995.

[0238] [8] J. Böhm, A.Ss Frangakis, R. Hegerl, S. Nickell, D. Typke, and W. Baumeister. Toward detecting and identifying macromolecules in a cellular context: template matching applied to electron tomograms. Proceedings of the National Academy of Sciences, Vol. 97 (No. 26): 14245–14250, 2000.

[0239] [9] AM Roseman. Particle finding in electron micrographs using a fastlocal correlation algorithm. Ultramicroscopy, Vol. 94 (Nos. 3-4): 225–236, 2003.

[0240]

[10] HD Ebbinghaus, H. Hermes, F. Hirzebruch, M. Koecher, K. Mainzer, J. Ncukirch, A. Prestel and R. Remmert. Numbers. Springer, 1991.

[0241]

[11] Lev Pontryagin. Generalization of numbers. CreateSpace, 2010.

[0242]

[12] P. Comon, G. Golub, L.-H. Lim and B. Mourrain. Symmetric tensors and symmetric tensor rank. SIAM Journal on Matrix Analysis and Applications, Vol. 30 (No. 3), 2008.

[0243]

[13] CFCui, Y.-H. Dai and J. Nie. All real eigenvalues ​​of symmetric tensors. SIAM Journal on Matrix Analysis and Applications, Vol. 35 (4): 1582-1601, 2014.

[0244]

[14] E. Kofidis and PARegalia. On the best rank-1 approximation of higher-order supersymmetric tensors. SIAM J. Matrix Anal. Appl., Vol. 23: 863-884, 2001.

[0245]

[15] TG Kolda and J. Mayo. Shifted power method for computing tensor eigenpairs. SIAM J. Matrix Anal. Appl., Vol. 32: 1095-1124, 2010.

[0246]

[16] I. Dokmanic and D. Petrinovic. Convolution on the n-sphere with application to pdf modeling. IEEE transactions on signal processing, 58(3): 1157–1170, 2009.

[0247]

[17] S. Cao, S. He, Z. Li and Z. Wang. Extreme ratio between spectral and frobenius norms of nonnegative tensors. SIAM Journal on Matrix Analysis and Applications, 44(2): 919–944, 2023.

[0248]

[18] K. Kozhasov and J. Tonelli-Cueto. Probabilistic bounds on best rank-one approximation ratio. Arxiv (arXiv preprint), 2022.

[0249]

[19] PARegalia and E. Kofidis. The higher-order power method revisited: convergence proofs and effective initialization. Included in the 2000 IEEE International Conference on Acoustics, Speech, and Signal Processing Proceedings (Cat. No. 00CH37100), Volume 5, pp. 2709-2712.

[0250]

[20] S. Banach. Über homogene polynome in (L 2 (in (L) 2 (homogeneous polynomials in ). Studio, Mathematica (Mathematical Research), Vol. 7 (Issue 1): 36-44, 1938.

[0251]

[21] CGKoay. A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere. J Comput Sci, Vol. 2 (4): 377–381, 2011.

[0252]

[22] ML. Chaillet, G. van der Nett, I. Gubins, S. Roet, R. C. Veltkamp, ​​and F. Foster. Extensive angular sampling enables the sensitive localization of macromolecules in electron tomograms. International Journal of Molecular Sciences, Vol. 24 (Issue 17): 13375, 2023.

Claims

1. A method for template matching of a tensor template of a target particle to identify instances of the target within a tomographic image obtained using electron tomography, such as cellular electron tomography, the method comprising: Obtain a tensor template of the target particle, wherein the tensor template uses a tensor field to describe the shape of the target particle under all rotations; as well as The tensor template is used to identify instances of the target particle within the tomographic image.

2. The method of claim 1, wherein obtaining the tensor template of the target particle comprises retrieving the tensor template from a stored tensor template library.

3. The method of claim 1, wherein obtaining a tensor template of the target particle comprises generating the tensor template, wherein the tensor template comprises a tensor field of a plurality of tensors, the tensor field utilizing the tensor field to describe the shape of the target particle over all rotations.

4. The method according to claim 3, wherein: The target particle is described by a model including a voxel array, and the tensor template includes a per-voxel tensor of the model; Templates are generated for different rotations of the target particle, each template including the voxel array; and A tensor power is generated for each rotation of the target particle, wherein the tensor power includes a component describing the shape of the target particle under that rotation.

5. The method of claim 4, wherein each tensor power is calculated as a tensor power of a quaternion describing the shape of the target particle under the rotation, optionally as a fourth-power multiplication of the four components of the quaternion.

6. The method of claim 4 or 5, wherein all tensors and all tensor powers have the same number of components, and each component of each tensor representing a voxel of the target particle comprises the integral over all rotations of the product of the corresponding component of the tensor power for one rotation and the value of the voxel of the template of the target particle for the one rotation.

7. The method according to any one of claims 1 to 6, the method comprising masking the model using an operator defining the model of the target particle before generating the components of the tensor, thereby setting the size of the voxel array, and optionally, wherein the operator is a symmetric positive semidefinite linear operator.

8. The method according to any one of claims 1 to 7, wherein an example of using the tensor template to identify the target particle within the tomographic image includes: The tensor template is convolved with the tomographic image to obtain a convolutional tensor field, which provides information about the correlation between the tensor template and the tomographic image at all locations and under all rotations within the tomographic image.

9. The method of claim 8, wherein an instance of using the tensor template to identify the target particle within the tomographic image comprises: The optimal correlation value for each location within the tomographic image is determined by computing a scalar map of the convolutional tensor field at each location.

10. The method of claim 9, wherein calculating a scalar map of the convolution tensor field at each location comprises using a norm, optionally a Frobenius norm, to calculate a scalar map of the total energy of the convolution tensor field at each location, and the method further comprises identifying the location of an instance of the target particle within the tomographic image by determining a local maximum of the optimal correlation value.

11. The method according to claim 10, further comprising: Using global optimization techniques, the rotation of the target particle at each of the identified locations can be identified, optionally by identifying the feature vectors and eigenvalues ​​of the convolution tensor field at each of the identified locations, or by sampling the rotation space of the target, represented as a quaternion, and calculating the inner product between a selected tensor of the convolution tensor field and the quaternion.

12. The method of claim 11, wherein an example of using the tensor template to identify the target particle within the tomographic image includes: The normalized correlation value of each identified instance of the target particle in the tomographic image is calculated by using the position and rotation identified for each instance of the target particle within the tomographic image as input to calculate the non-tensor template of the target particle and the local normalized cross-correlation of the tomographic image. Optionally, The method further includes generating each non-tensor template based on the model of the target particle, as a template describing the shape of the target particle at the corresponding identified rotation.

13. The method according to any one of claims 1 to 12, further comprising: The tomographic image is segmented into blocks; as well as For each block in turn: The tensor template is convolved with the block of the tomographic image to obtain a convolution tensor field, which provides information about the correlation between the tensor template and the tomographic image at all positions and under all rotations within the block of the tomographic image. Determine the optimal correlation value at each of the locations within the block of the tomographic image, thereby calculating a scalar map from the convolutional tensor field; The location of an instance of the target particle within that block of the tomographic image is identified by determining the local maximum of the optimal correlation value; Global optimization techniques are used to identify the rotation of the target particle at each of the identified positions; and optionally, Before calculating the normalized correlation value for each identified instance of the target particle within the tomographic image, the rotation and position identified in each block of the tomographic image are combined.

14. The method of claim 13, wherein segmenting the tomographic image into blocks comprises segmenting the tomographic image into continuous units and defining each block as comprising a unit plus an overlapping portion extending into the adjacent unit, and optionally, the size of the overlapping portion extending into the adjacent unit is less than the maximum size of the target particle and greater than half the test size of the target particle.

15. A computer program comprising computer program instructions that, when executed by a computer processor, cause the computer processor to perform the method according to any of the preceding claims; A computer-readable medium on which such a computer program is stored; or a computer system comprising a computer processor and a computer memory in which such a computer program is stored.

Images