Tensorial template matching

Abstract

A method of template matching a tensorial template of a target particle to identify instances of that target within a tomogram comprises obtaining a tensorial template of the target particle, wherein the tensorial template describes the shape at all rotations of the target particle with a tensor field, and using the tensorial template to identify instances of that target particle within a tomogram.

Description

[0001] TENSORIAL TEMPLATE MATCHING

[0002] FIELD OF THE INVENTION

[0003] The present invention relates to an improved method of template matching in image analysis. In particular, the present invention relates to template matching of target molecules in cellular cryogenic electron tomography that may be performed faster than before.

[0004] BACKGROUND OF THE INVENTION

[0005] The present invention relates to template matching where a scan is performed of an image to determine matches to a target template. The present invention can be applied to 2D images and 3D images (volumetric images or tomograms), although this invention is considerably advantageous in the case of 3D images. The present invention may be used in template matching of a target that is present many times in an image, for example a particle within an image, e.g., ribosomes or proteasomes. More generally, template matching is widely used in the industry e.g., for image registration, alignment, object detection, counting, tracking, defect detection, volume reconstruction, etc. It is often used in such applications due its ease of use and interpretation.

[0006] A particular application of the present invention is cellular cryogenic electron tomography (Cryo-ET) which is a label-free imaging technique that provides 3D images (or tomograms) of organelles and protein complexes at nanometre resolution. This is done by opening windows into the cell with focused ion beam milling of a cryogenically frozen cell. A series of 2D images is taken of this thinned cellular sample from varying angles. The 2D images are then reconstructed into a 3D tomogram. This technique provides 3D snapshots of proteins at work within their functional environment, and provides a window into understanding how they, and other molecules, work together to carry out major processes in a cell. Cryo-ET serves two main purposes, namely protein structures discovery and functional cellular behaviour. In both cases, particles of interest are automatically picked to get their position, orientation, and class, usually using template matching.

[0007] Template matching is a brute force algorithm that computes correlation with a template of a target particle at every voxel of a tomogram. Voxels showing the highest correlation values correspond to instances of template positions in the tomogram having the best match to a particle recorded in the tomogram. In addition, because every instance of a template position may have any rotation of the target particle, the whole space of rotations of the target particle must be sampled and the correlation must be computed for all the rotations to find the optimal rotation at every voxel in the tomogram. This approach is computationally intensive and requires a significant amount of time to be computed even after algorithmic optimization, usage of the Fourier transform, parallelization, and down-sampling of the data. This leads to a strong bottleneck in the Cryo-ET analysis workflow.

[0008] SUMMARY OF THE INVENTION

[0009] There is described herein a method of generating a tensorial template for template matching to determine instances of a target within an image; the method comprising: generating a tensorial template of the target, wherein the tensorial template describes the shape over all rotations of the target with a tensor field.

[0010] Optionally, the tensorial template comprises a tensor field of multiple tensors. The tensors may be of any order and any dimension according to a compromise between rotation accuracy and computational complexity required for the targeted application. For example, tensors of order 4 and dimension 4 may be used.

[0011] The target may be described by a model comprising an array of elements, for example pixels in a two-dimensional image or voxels in a three-dimensional image. The tensorial template may comprise a tensor per element of the model. The templates may be generated for different rotations of the target, each template comprising the array of elements. Optionally, a tensor power is generated for each rotation of the target and the tensor power may comprise components that describe the shape of the target at that rotation. Each tensor power may be calculated as a tensor power of a quaternion that describes the shape of the target at that rotation. Each tensor power may be calculated as a fourth power multiplication of the four components of a quaternion.

[0012] Preferably, all tensors and all tensor powers have the same number of components, and each component of each tensor that represents one 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 a value of the one element of the template of the target for the one rotation. The value of each element may denote an intensity of that element. For example, the image may comprise a surface model of the target that used intensity values for each element to describe the shape of the surface of the target.

[0013] Optionally, the method comprises masking the model using an operator that bounds the model of the target prior to generating the components of the tensors thereby setting the size of the array of elements. The operator may be a symmetric positive semidefinite linear operator.

[0014] According to a first aspect of the present invention, there is provided a method of tensorial template matching to identify instances of a target within an image, the method comprising: generating a tensorial template of the target according to any of the methods described above; and using the tensorial template to identify instances of that target within an image. According to a second aspect, there is provided a method of template matching a template of a target to identify instances of that target within an image, the method comprising: obtaining a tensorial template of the target, wherein the tensorial template describes the shape at all rotations of the target with a tensor field; and using the tensorial template to identify instances of that target within an image. Obtaining a tensorial template of the target may comprise generating the tensorial template or retrieving the tensorial template from a library of saved tensorial templates.

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

[0016] Optionally, using the tensorial template to identify instances of the target within an image comprises: convolving the tensorial template with the image to obtain a convolved tensor field that provides information on the correlation of the tensorial template with the image at all positions and at all rotations within the image. Using the tensorial template to identify instances of the target within an image may comprise determining the best correlation value for each of the positions within the image by calculating a scalar map of the convolved tensor field at each position. The values populating the scalar map may be a measure of the best correlation value at each position. Each position may correspond to an element of the image, for example pixels in a two-dimensional image or voxels in a three-dimensional image. Calculating a scalar map of the convolved tensor field at each position may comprise calculating a scalar map of the total energy of the convolved tensor field at each position. The method may comprise calculating a scalar map of the total energy of the convolved tensor field at each position using a norm of the convolved tensor field at each position. A Frobenius norm may be used. Calculating the Frobenius norm of the convolved tensor field at each position may comprise calculating the square root of the absolute squares of the components of the tensors of the tensorial template at each position.

[0017] The method may further comprise identifying positions of instances of the target within the image by determining local maxima of the best correlation values. Identifying rotations of the target at each of the identified positions may comprise using a global optimisation technique.

[0018] Identifying rotations of the target at each of the identified positions may comprise identifying eigenvectors and eigenvalues of the convolved tensor field at each of the identified positions. Identifying rotations of the target at each of the identified positions by identifying eigenvectors and eigenvalues of tensors in the convolved tensor field at each of the identified positions may comprise unfolding tensors in the convolved tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues. Performing operations on the matrix to determine the eigenvectors and eigenvalues comprises using several calls to shifted symmetric higher order power method restarted at different random rotations. Identifying rotations of the target at each of the identified positions using a global optimisation technique may comprise sampling the space of rotations of the target expressed as quaternions and calculating the inner product between the selected tensors of the convolved tensor field and the quaternions. This provides a list of rotations that can be ranked by the value of this product. A list of peak positions with their most promising rotations (highest inner products) per peak may be returned.

[0019] Optionally, the method may further comprise identifying positions and rotations of the target using a global optimisation technique. The positions and rotations may be identified in a single step, rather than identifying the positions and the rotations sequentially as may be done for the methods described above. Positions and rotations of the target may be identified by identifying eigenvectors and eigenvalues of the convolved tensor field. Identifying positions and rotations of the target by identifying eigenvectors and eigenvalues of tensors in the convolved tensor field may comprise unfolding tensors in the convolved tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues. Performing operations on the matrix to determine the eigenvectors and eigenvalues may comprise using several calls to shifted symmetric higher order power method restarted at different random rotations. Identifying rotations of the target at each of the identified positions may comprise sampling the space of rotations of the target expressed as quaternions and calculating the inner product between the selected tensors of the convolved tensor field and the quaternions. This may provide a list of rotations that can be ranked by the value of this product. A list of peak positions with their most promising rotations (highest inner products) per peak may be returned.

[0020] Preferably, using the tensorial template to identify instances of that target within an image comprises calculating normalised correlation values for each identified instance of the target within the image by local normalised cross-correlation of a non-tensorial template of the target with the image using, as inputs, the position and rotation identified for each instance of the target within the image. Advantageously, this step produces normalised cross-correlation values that allow better comparison of instances of identified positions and rotations. While calculating local normalised cross-correlation values is usually very time consuming, it is not in this instance as only the values for the identified positions and rotations need to be calculated, which massively reduces the computational space (from all positions and all rotations as is performed in conventional local normalised cross-correlation). The method may comprise generating each non-tensorial template from a model of the target to be a template describing the shape of the target at the corresponding identified rotation.

[0021] Optionally, the method further comprises providing a list of the identified instances of the target within the image. The list may comprise, for each instance of the target, the calculated local normalised cross-correlation value, the identified position and the identified rotation. The list may be limited according to a preference selected by a user.

[0022] The method may further comprise dividing the image into blocks. For each block in turn, the method may comprise convolving the tensorial template with that block of the image to obtain a convolved tensorial field that provides information on the correlation of the tensorial template with the image at all positions and all rotations within that block of the image. Then, the best correlation value for each of the positions within that block of the image may be determined thereby calculating a scalar map from the convolved tensor field. Then, the positions of instances of the target within that block of the image may be identified by determining local maxima of the best correlation values. Finally, rotations of the target at each of the identified positions may be identified using a global optimisation technique.

[0023] Hence, the previously described steps of convolving the tensorial template, determining correlation values for each position, and identifying positions and rotations of each instance of the target in the image as a whole is broken down into a series of routines performed on each block in turn. This can be advantageous where computer memory is limited.

[0024] Dividing the image into blocks may comprise dividing the image into contiguous units and defining each bock to comprise one unit plus an overlap portion extending into the adjacent units. The overlap may be a unform distance and the size of the overlap portion extension into adjacent units may be less than the greatest dimension of the target and may be more than half the test dimension of the target. The method may comprise, for each block, keeping positions identified within the unit of that block and discarding any positions identified within the overlap portion of that block.

[0025] The method may further comprise merging the rotations and positions identified in each block of the image prior to calculating the normalised cross correlation values calculated for each identified instance of the target within the image.

[0026] Optionally, the target is a particle and the image is a tomogram. The tomogram may be obtained using cellular cryogenic electron tomography. Alternatively, the target may be a defect and the image may be an image of a manufactured part. The defect may be a defect in a semiconductor material and the manufactured part may be a semiconductor part such as a chip or wafer.

[0027] There is also described a method of counting objects within an image, comprising any of the method described above, and providing the number of identified instances of the target in the image as an output. There is also described a method of tracking an object within a succession of images, comprising performing any of the methods described above to identify an instance of the target in each image of the succession of images. The present invention also resides in a computer program comprising computer program instructions that, when executed by a computer processor, cause the computer processor to perform any of the methods described above, in a computer readable medium having such a computer program stored thereon and in a computer system comprising a computer processor and computer memory having such a computer program stored therein.

[0028] The following clauses describe further features.

[0029] 1. A method of generating a tensorial template for template matching to determine instances of a target within an image; the method comprising: generating a tensorial template of the target, wherein the tensorial template describes the shape over all rotations of the target with a tensor field.

[0030] 2. The method of clause 1, wherein the tensorial template comprises a tensor field of multiple tensors.

[0031] 3. The method of clause 2, wherein the target is described by a model comprising an array of elements, and the tensorial template comprises a tensor per element of the model.

[0032] 4. The method of clause 3, wherein templates are generated for different rotations of the target, each template comprising the array of elements.

[0033] 5. The method of clause 4, wherein a tensor power is generated for each rotation of the target and wherein the tensor power comprises components that describe the shape of the target at that rotation.

[0034] 6. The method of clause 5, wherein each tensor power is calculated as a tensor power of a quaternion that describes the shape of the target at that rotation.

[0035] 7. The method of clause 6, wherein each tensor power is calculated as a fourth power multiplication of the four components of a quaternion.

[0036] 8. The method of any of clauses 5 to 7, wherein all tensors and all tensor powers have the same number of components, and each component of each tensor that represents one 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 a value of the one element of the template of the target for the one rotation.

[0037] 9. The method of clause 8, wherein the value of each element denotes an intensity of that element.

[0038] 10. The method of any of clauses 3 to 9, comprising masking the model using an operator that bounds the model of the target prior to generating the components of the tensors thereby setting the size of the array of elements.

[0039] 11. The method of clause 10, wherein the operator is a symmetric positive semidefinite linear operator. 12. A method of tensorial template matching to identify instances of a target within an image, the method comprising: generating a tensorial template of the target according to the method of any of clauses 1 to 11; and using the tensorial template to identify instances of that target within an image.

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

[0041] 14. The method of clause 13, wherein obtaining a tensorial template of the target comprises generating the tensorial template or retrieving the tensorial template from a library of saved tensorial templates.

[0042] 15. The method of any of clauses 12 to 14, wherein using the tensorial template to identify instances of the target within an image comprises: convolving the tensorial template with the image to obtain a convolved tensor field that provides information on the correlation of the tensorial template with the image at all positions and at all rotations within the image.

[0043] 16. The method of clause 15, wherein using the tensorial template to identify instances of the target within an image comprises: determining the best correlation value for each of the positions within the image by calculating a scalar map of the convolved tensor field at each position.

[0044] 17. The method of clause 16, wherein calculating a scalar map of the convolved tensor field at each position comprises calculating a scalar map of the total energy of the convolved tensor field at each position.

[0045] 18. The method of clause 17, comprising calculating a scalar map of the total energy of the convolved tensor field at each position using a norm, optionally the Frobenius norm, of the convolved tensor field at each position and, optionally, wherein calculating the Frobenius norm of the convolved tensor field at each position comprises calculating the square root of the absolute squares of the components of the tensors of the tensorial template at each position.

[0046] 19. The method of any of clauses 16 to 18, further comprising: identifying positions of instances of the target within the image by determining local maxima of the best correlation values.

[0047] 20. The method of clause 19, further comprising: identifying rotations of the target at each of the identified positions using a global optimisation technique.

[0048] 21. The method of clause 19 or 20, further comprising: identifying rotations of the target at each of the identified positions by identifying eigenvectors and eigenvalues of the convolved tensor field at each of the identified positions. 22. The method of clause 21, wherein identifying rotations of the target at each of the identified positions by identifying eigenvectors and eigenvalues of tensors in the convolved tensor field at each of the identified positions comprises unfolding tensors in the convolved tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0049] 23. The method of clause 22, wherein performing operations on the matrix to determine the eigenvectors and eigenvalues comprises using several calls to shifted symmetric higher order power method restarted at different random rotations.

[0050] 24. The method of clause 20, further comprising: identifying rotations of the target at each of the identified positions by sampling the space of rotations of the target expressed as quaternions, and calculating the inner product between the selected tensors of the convolved tensor field and the quaternions.

[0051] 25. The method of clause 15, further comprising: identifying positions and rotations of the target using a global optimisation technique.

[0052] 26. The method of clause 27, further comprising: identifying positions and rotations of the target by identifying eigenvectors and eigenvalues of the convolved tensor field.

[0053] 27. The method of clause 26, wherein identifying positions and rotations of the target by identifying eigenvectors and eigenvalues of tensors in the convolved tensor field comprises unfolding tensors in the convolved tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0054] 28. The method of clause 27, wherein performing operations on the matrix to determine the eigenvectors and eigenvalues comprises using several calls to shifted symmetric higher order power method restarted at different random rotations

[0055] 29. The method of clause 19 or 20, further comprising: identifying rotations of the target at each of the identified positions by sampling the space of rotations of the target expressed as quaternions, and calculating the inner product between the selected tensors of the convolved tensor field and the quaternions.

[0056] 30. The method of any of clauses 20 to 29, wherein using the tensorial template to identify instances of that target within an image comprises: calculating normalised correlation values for each identified instance of the target within the image by calculating the local normalised cross-correlation of a non-tensorial template of the target with the image using, as inputs, the position and rotation identified for each instance of the target within the image.

[0057] 31. The method of clause 30, wherein calculating normalised correlation values for each identified instance of the target within the image comprises calculating the local normalised cross-correlation of the non-tensorial template of the target with the image for only the positions and rotations identified for the instances of the target within the image.

[0058] 32. The method of clause 30 or 31, further comprising generating each non-tensorial template from a model of the target to be a template describing the shape of the target at the corresponding identified rotation.

[0059] 33. The method of any of clauses 20 to 32, further comprising providing a list of the identified instances of the target within the image; and, optionally, wherein: the list comprises, for each instance of the target, the calculated local normalised cross-correlation value, the identified position and the identified rotation; and / or the list is limited according to a preference selected by a user.

[0060] 34. The method of any of clauses 12 to 33, further comprising: dividing the image into blocks; and, for each block in turn: (i) convolving the tensorial template with that block of the image to obtain a convolved tensor field that provides information on the correlation of the tensorial template with the image at all positions and at all rotations within that block of the image; (ii) determining the best correlation value for each of the positions within that block of the image thereby calculating a scalar map from the convolved tensor field; (iii) identifying positions of instances of the target within that block of the image by determining local maxima of the best correlation values; and (iv) identifying rotations of the target at each of the identified positions using a global optimisation technique.

[0061] 35. The method of clause 34, wherein dividing the image into blocks comprises dividing the image into contiguous units and defining each bock to comprise one unit plus an overlap portion extending into the adjacent units.

[0062] 36. The method of clause 35, wherein the overlap is a unform distance and, optionally, the size of the overlap portion extension into adjacent units is less than the greatest dimension of the target and more than half the test dimension of the target.

[0063] 37. The method of clause 36 comprising, for each block, keeping positions identified within the unit of that block and discarding any positions identified within the overlap portion of that block.

[0064] 38. The method of any of clauses 34 to 37, further comprising merging the rotations and positions identified in each block of the image prior to calculating the normalised correlation values calculated for each identified instance of the target within the image.

[0065] 39. The method of any preceding clause, wherein the target is a particle and the image is a tomogram, optionally a tomogram obtained using cellular cryogenic electron tomography.

[0066] 40. The method of any of clauses 1 to 38, wherein the target is a defect and the image is an image of a manufactured part. 41. The method of 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.

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

[0068] 43. A method of tracking an object within a succession of images, comprising performing the method of any of clauses 12 to 41 to identify an instance of the target in each image of the succession of images.

[0069] 44. A computer program comprising computer program instructions that, when executed by a computer processor, cause the computer processor to perform the method of any preceding clause, a computer readable medium having such a computer program stored thereon or as computer system comprising a computer processor and computer memory having such a computer program stored therein.

[0070] 101. A method of template matching a tensorial template of a target particle to identify instances of that target within a tomogram, the method comprising: obtaining a tensorial template of the target particle, wherein the tensorial template describes the shape at all rotations of the target particle with a tensor field; and using the tensorial template to identify instances of that target particle within a tomogram.

[0071] 102. The method of clause 101, wherein obtaining a tensorial template of the target particle comprises retrieving the tensorial template from a library of saved tensorial templates.

[0072] 103. The method of clause 101, wherein obtaining a tensorial template of the target particle comprises generating the tensorial template, wherein the tensorial template describes the shape over all rotations of the target particle with a tensor field.

[0073] 104. The method of clause 103, wherein the tensorial template comprises a tensor field of multiple tensors.

[0074] 105. The method of clause 104, wherein the target particle is described by a model comprising an array of voxels, and the tensorial template comprises a tensor per voxel of the model.

[0075] 106. The method of clause 105, wherein templates are generated for different rotations of the target particle, each template comprising the array of voxels.

[0076] 107. The method of clause 106, wherein a tensor power is generated for each rotation of the target particle and wherein the tensor power comprises components that describe the shape of the target particle at that rotation.

[0077] 108. The method of clause 107, wherein each tensor power is calculated as a tensor power of a quaternion that describes the shape of the target particle at that rotation. 109. The method of clause 108, wherein each tensor power is calculated as a fourth power multiplication of the four components of a quaternion.

[0078] 110. The method of any of clauses 107 to 109, wherein all tensors and all tensor powers have the same number of components, and each component of each tensor that represents one 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 a value of the one voxel of the template of the target particle for the one rotation.

[0079] 111. The method of clause 110, wherein the value of each voxel denotes an intensity of that voxel.

[0080] 112. The method of any of clauses 105 to 111, comprising masking the model using an operator that bounds the model of the target particle prior to generating the components of the tensors thereby setting the size of the array of voxels.

[0081] 113. The method of clause 112, wherein the operator is a symmetric positive semidefinite linear operator.

[0082] 114. The method of any of clauses 101 to 113, wherein using the tensorial template to identify instances of the target particle within the tomogram comprises: convolving the tensorial template with the tomogram to obtain a convolved tensor field that provides information on the correlation of the tensorial template with the tomogram at all positions and at all rotations within the tomogram.

[0083] 115. The method of clause 114, wherein using the tensorial template to identify instances of the target particle within the tomogram comprises: determining the best correlation value for each of the positions within the tomogram by calculating a scalar map of the convolved tensor field at each position.

[0084] 116. The method of clause 115, wherein calculating a scalar map of the convolved tensor field at each position comprises calculating a scalar map of the total energy of the convolved tensor field at each position.

[0085] 117. The method of clause 116, comprising calculating a scalar map of the total energy of the convolved tensor field at each position using a norm, optionally the Frobenius norm, of the convolved tensor field at each position.

[0086] 118. The method of clause 117, comprising calculating the scalar map using the Frobenius norm, wherein calculating the Frobenius norm of the convolved tensor field at each position comprises calculating the square root of the absolute squares of the components of the tensors of the tensorial template at each position.

[0087] 119. The method of any of clauses 115 to 118, further comprising: identifying positions of instances of the target particle within the tomogram by determining local maxima of the best correlation values. 120. The method of clause 119, further comprising: identifying rotations of the target particle at each of the identified positions using a global optimisation technique.

[0088] 121. The method of clause 1119 or 20, further comprising: identifying rotations of the target particle at each of the identified positions by identifying eigenvectors and eigenvalues of the convolved tensor field at each of the identified positions.

[0089] 122. The method of clause 121, wherein identifying rotations of the target particle at each of the identified positions by identifying eigenvectors and eigenvalues of tensors in the convolved tensor field at each of the identified positions comprises unfolding tensors in the convolved tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0090] 123. The method of clause 122, wherein performing operations on the matrix to determine the eigenvectors and eigenvalues comprises using several calls to shifted symmetric higher order power method restarted at different random rotations.

[0091] 124. The method of clause 120, further comprising: identifying rotations of the target particle at each of the identified positions by sampling the space of rotations of the target expressed as quaternions, and calculating the inner product between the selected tensors of the convolved tensor field and the quaternions.

[0092] 125. The method of clause 115, further comprising: identifying positions and rotations of the target using a global optimisation technique.

[0093] 126. The method of clause 127, further comprising: identifying positions and rotations of the target particle by identifying eigenvectors and eigenvalues of the convolved tensor field.

[0094] 127. The method of clause 126, wherein identifying positions and rotations of the target particle by identifying eigenvectors and eigenvalues of tensors in the convolved tensor field comprises unfolding tensors in the convolved tensor field into a matrix and performing operations on the matrix to determine the eigenvectors and eigenvalues.

[0095] 128. The method of clause 127, wherein performing operations on the matrix to determine the eigenvectors and eigenvalues comprises using several calls to shifted symmetric higher order power method restarted at different random rotations.

[0096] 129. The method of clause 119 or 120, further comprising: identifying rotations of the target particle at each of the identified positions by sampling the space of rotations of the target expressed as quaternions, and calculating the inner product between the selected tensors of the convolved tensor field and the quaternions.

[0097] 130. The method of any of clauses 120 to 129, wherein using the tensorial template to identify instances of the target particle within the tomogram comprises: calculating normalised correlation values for each identified instance of the target particle within the tomogram by calculating the local normalised cross-correlation of a non-tensorial template of the target particle with the tomogram using, as inputs, the position and rotation identified for each instance of the target particle within the tomogram.

[0098] 131. The method of clause 130, wherein calculating normalised correlation values for each identified instance of the target particle within the tomogram comprises calculating the local normalised cross-correlation of the non-tensorial template of the target particle with the tomogram for only the positions and rotations identified for the instances of the target particle within the tomogram.

[0099] 132. The method of clause 130 or 131, further comprising generating each non-tensorial template from a model of the target particle to be a template describing the shape of the target particle at the corresponding identified rotation.

[0100] 133. The method of any of clauses 120 to 132, further comprising providing a list of the identified instances of the target particle within the tomogram; and, optionally, wherein: the list comprises, for each instance of the target particle, the calculated local normalised cross-correlation value, the identified position and the identified rotation; and / or the list is limited according to a preference selected by a user.

[0101] 134. The method of any of clauses 101 to 133, further comprising: dividing the tomogram into blocks; for each block in turn: (i) convolving the tensorial template with that block of the tomogram to obtain a convolved tensor field that provides information on the correlation of the tensorial template with the tomogram at all positions and at all rotations within that block of the tomogram; (ii) determining the best correlation value for each of the positions within that block of the tomogram thereby calculating a scalar map from the convolved tensor field; (iii) identifying positions of instances of the target particle within that block of the tomogram by determining local maxima of the best correlation values; and (iv) identifying rotations of the target particle at each of the identified positions using a global optimisation technique.

[0102] 135. The method of clause 134, wherein dividing the tomogram into blocks comprises dividing the tomogram into contiguous units and defining each block to comprise one unit plus an overlap portion extending into the adjacent units.

[0103] 136. The method of clause 135, wherein the overlap is a unform distance and, optionally, the size of the overlap portion extension into adjacent units is less than the greatest dimension of the target particle and more than half the test dimension of the target particle.

[0104] 137. The method of clause 136 comprising, for each block, keeping positions identified within the unit of that block and discarding any positions identified within the overlap portion of that block. 138. The method of any of clauses 134 to 137, further comprising merging the rotations and positions identified in each block of the tomogram prior to calculating the normalised correlation values calculated for each identified instance of the target particle within the tomogram.

[0105] 139. A computer program comprising computer program instructions that, when executed by a computer processor, cause the computer processor to perform the method of any of clauses 101 to 138, a computer readable medium having such a computer program stored thereon or as computer system comprising a computer processor and computer memory having such a computer program stored therein.

[0106] LIST OF FIGURES

[0107] In order that the invention can be more readily understood, reference will now be made by way of example only, to the accompanying drawings in which:

[0108] Figure 1A shows a cryo-tomogram of a Chlamydomonas cell and Figure IB shows surface models of the templates used to detect and identify the proteins shown in Figure 1A;

[0109] Figure 2 is a block diagram of a method of identifying molecules in a tomogram;

[0110] Figure 3 is a schematic representation of the method of identifying molecules in a tomogram of Figure 2;

[0111] Figure 4 is a flow diagram of a method of generating a tensorial template; and

[0112] Figure 5 is a schematic representation of an alternative method of identifying molecules in a tomogram.

[0113] DETAILED DESCRIPTION

[0114] Figure 1A shows a cryo-tomogram 10 of a Chlamydomonas cell that includes structures 15 like Golgi apparatus 151, Endoplasmic Reticulum 152 and nuclear membrane 153. The tomogram 10 also contains many large molecules 15, most of which are ribosomes 154, 155, and proteasomes and CDC48-ATPase 156. These molecules 15 have been detected by template matching. Figure IB shows surface models of templates 20 for particular rotations of large molecules 15 like those shown in Figure 1A. These templates 20 are derived from computer-generated models of the large molecules 15 and can be generated as a representation of a particular molecule 15 having a particular rotation. These templates 20 can be used to detect and identify instances of the large molecules 15 in the tomogram 10 of Figure 1A using the conventional template matching described above.

[0115] Figures 2 to 4 show a method 200 of identifying target molecules within a tomogram 10, such as at least some of the large molecules 15 of Figure IB within the tomogram 10 of Figure 1A. In Figure 3, the tomogram 10 is shown as a two-dimensional representation for the sake of simplicity and clarity. The method 200 begins by obtaining a tensorial template 222 of a target molecule. The tensorial template 222 describes the target molecule for all rotations with a tensor field. In traditional template matching, the templates 20 are generated of the target molecule in which the model of the target molecule is reduced to a finite but considerable number of snapshots of the target molecule in different rotations (e.g., tens to hundreds of thousands depending on the targeted accuracy). The tensorial template 222 described herein replaces traditional templates 20 with a tensor field of limited tensor order embedding all rotations of the target molecule.

[0116] A model of a target molecule of interest is obtained at step 202 and a tensorial template 222 is calculated for that target molecule at step 204. The tensorial template 222 describes the shape of the target molecule over all rotations using a tensor field that is generated as follows.

[0117] Many different rotations of the target molecule are sampled, with a mathematical description of each rotation being generated. These descriptions are then combined mathematically using an integration to derive each tensor of the tensorial template 222 that describes the shape of the target molecule over all rotations. A tensor is calculated for each voxel of model of the target molecule, hence forming the tensor field.

[0118] Each tensorial template 222 is a tensor field that embeds all rotations of a target molecule with tensors that have as many independent components as a tensor of order n and dimension d, e.g. 35 for a tensor of order 4 and dimension 4 (the number of components of the tensor depends upon 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 of the sampled rotations, a tensor power is found as the n tensor power of q, where n is the order of the tensor and q is a vector of dimension d representing the rotation. In this embodiment, the quaternions are used as the vectors. Each quaternion is transformed into the tensor power that describes the corresponding rotation through a 4th power multiplication of the four components of each quaternion, which produces tensor powers with a dimension of 4. This dimension is preserved in the calculations that are performed to determine each tensor of the tensorial template 222. This leaves a choice in the order used for the tensor powers and hence in the tensors. This is a compromise in that the higher order of tensors used, the better accuracy that may be achieved. However, the computational cost also increases with the order. Using 4th-order, 4 dimension tensors provides a good compromise and this results in a tensorial template 222 with 35 linearly independent components (i.e. 35 distinct combinations of the four components of the quaternion).

[0119] To combine the tensor powers derived from each individual quaternion that each describe a single rotation of the target molecule, an integration is performed for each voxel of the model of the target molecule over all the sampled rotations, thereby deriving a voxel-by-voxel representation of the target molecule over all rotations. To do this, for each voxel, each of the 35 components of the tensor powers are combined with the corresponding rotated version of the model of the target molecule to provide a corresponding component in a tensor of the tensorial field. Consequently, the tensor of each voxel is created, and each tensor has as many components as the n order tensor describing the rotations (35 in this example). Specifically, each component of each tensor is calculated as the integral over all rotations of the product between the corresponding component of the tensor power, generated by the n power of the quaternion defining the rotation, and the corresponding intensity value of the voxel of the corresponding rotated model of the target molecule.

[0120] An example of a method 400 of calculating the tensorial template 222 is shown in Figure 4. In this example, the method starts at step 252 where a tensorial field is created. The tensorial field is provided with a tensor for each voxel, and each tensor has 35 components that are populated according to the following method steps that comprises three nested loops.

[0121] The first loop 402 is entered at step 254 when a first sampled rotation of the target molecule is selected. That sampled rotation is represented by one of the quaternions and, at step 256, the quaternion for the current rotation is translated into a corresponding tensor power. As described above, this tensor power will have 35 components. At step 258, a template 20 of the target molecule having the current rotation is generated from the model of the target molecule.

[0122] Next, the second loop 404 is entered at step 260 where the tensor of the first voxel is selected. The third loop 404 populates the components of this tensor in turn and, initially, a first of the 35 components of this tensor is selected at step 262. Then, the product of the corresponding component of the tensor power for the first sampled rotation and the intensity of the first voxel of the template 20 of the target molecule for the first sampled rotation is calculated at step 264. The value of this product is a scalar and, at step 266, that value is stored as the first component of the tensor of the first voxel that was selected at step 260.

[0123] With the first component of the tensor populated, the process is then repeated for each of the other 34 components of the tensor by iterating through multiple passes of loop 406, via loop path 268 and 269 and steps 262 to 266. That is, a new component of the tensor is selected at step 262, and 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 at step 264. The product is then added to the current component of the tensor at step 266. This is repeated for all the 35 components of the tensor such that each component of the tensor is populated with an initial scalar value. With this done, the test at step 268 determines that all components have been populated, and so the next tensor may be selected.

[0124] The next tensor is selected at step 260, and the third loop 406 repeated for the new tensor such that all the components of that tensor are populated. Each component of the new tensor is calculated as the method iterates through loop 406. Loop 404 ensures that this is repeated for each tensor, until the tensors for all the voxels have been populated with initial scalar values corresponding to the products of the corresponding component of the tensor powers and the intensity of the corresponding voxel for the first rotation of the target molecule. When the test at step 270 determines that all tensors have been processed once, we have a set of tensors with all tensors having values added for all 35 components of each tensor. These values correspond to the first sampled rotation of the target molecule though. Repeated iterations through loop 401 repeat this process for each sampled rotation of the target molecule.

[0125] Each time the method 400 returns to step 254, via loop path 272 and 273, a new rotation is selected. Then, at steps 256 and 258, the tensor power and representation of the target molecule at the new rotation are generated. The method then proceeds to iterate through the components of each tensor, one tensor at a time via loops 406 and 404. At each pass through 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 comprises 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 result in the scalar values stored in the tensors being updated to add the scalar values being determined. Therefore, as the method 400 completes over all iterations of the sampled rotations, each component of each tensor becomes the integration of all the products of the corresponding components of the tensor powers with the intensity of the corresponding voxel over all rotations. It is this set of tensors that make up the tensor field.

[0126] In summary: the tensorial template 222 is the tensorial field comprising the set of tensors for the voxels of the model of the target molecule; each voxel has an associated tensor; each tensor comprises the same number of components as the tensor powers, 35 components in this example; each component of each tensor is a scalar value; and each scalar value equals the integration of each of the products of the corresponding component of the tensor power and the intensity of the corresponding voxel in each rotated model of the target.

[0127] Before computing the integral, the model of the target molecule may be processed by a symmetric positive semidefinite linear operator which includes a spherical mask that defines the local neighbourhood of the model. This limits the number of voxels and provides a further advantage in a subsequent step of the method 200 that is described below.

[0128] Although its computation is not instantaneous (a few minutes to few hours, depending on the number of rotations sampled and the number of voxels in the model of the target molecule), it need be performed only once per target molecule (although it may be updated from time to time as the model of the target molecule is refined). Once calculated, the tensorial template 222 may be saved in a library for future use. Hence, as shown in an alternative path in Figure 2, the method 200 may start at step 206 where a tensorial template 222 of a target molecule 20 is simply retrieved from a library. Although not depicted in Figures 2 and 3, the model of the target molecule may be masked with a sphere enclosing tightly but completely the target molecule before computing the tensorial template 222. Masking the target molecule model helps to produce a functional tensorial template 222.

[0129] At step 208, the tensorial template 222, however obtained, is convolved with the tomogram 10 to produce a convolved tensor field 226 that contains information about the correlation of the tensorial template 222 in the whole space of rotation used to generate the tensorial template 222 with each voxel of the tomogram 10. Although not depicted in Figures 2 and 3, where a spherical template mask is used to produce the tensorial template 222, this also focusses the correlation to a neighbourhood approximating the target molecule size at every voxel of the tomogram 10. This convolved tensor field 226 may then be used to determine the best image matches for instances of the target molecule in the tomogram 10, both in terms of voxel position within the tomogram 10 and in terms of the rotation of the target molecule. To identify instances of the target molecule, the following multi-step process is performed.

[0130] In a first step 210, a scalar value is calculated for each voxel for the tomogram 10 from the convolved tensor field 226. These correlation values are calculated using a scalar valued function that has a critical value that indicates the best correlation position. The result of this calculation is a scalar map 228 that provides an estimation of the value of the best correlation of the tensorial template 222 at each voxel. In Figure 3, the scalar map 228 is shown as a two-dimensional representation for the sake of simplicity and clarity, although it will be understood that the scalar map 228 is three dimensional as it provides the estimated best correlation value at each voxel. A norm may be used to calculate the scalar values , for example a Frobenius norm a shown at step 210 of Figure 2. The Frobenius norm of a tensor is computed as the square root of the sum of the absolute squares of its components, so its computation is very efficient.

[0131] Step 210 produces a scalar map having a scalar value for each voxel of the tomogram. To determine the best candidates for instances of the target molecule in the tomogram 10, peak positions of peak values within the scalar map 228 are determined. This is found by determining the positions of the local maxima of the correlation values in the scalar map 228 using any conventional peak finding technique.

[0132] Next, at step 214, the rotation of the target molecule that produces the best correlation value is determined for each of the potential peak positions found in step 212. This may be performed using an eigen decomposition or by sampling the space of rotations defined by the quaternions, expanding the corresponding quaternions to n-order tensors and computing their inner product with respect to the convolved tensor field. Step 214 of Figure 2 shows an example that uses an eigen decomposition. As noted above, the convolved tensor field 226 retains the information about all rotations, and eigenvector and eigenvalues can be used to determine the best match positions and corresponding rotations, i.e., those that produce the highest correlation values. This is performed by unfolding the tensors of the convolved tensor field 226 into a matrix and determining eigenvalues and eigenvectors from the matrix with each eigenvalue corresponding to the best correlation value. The eigenvector and eigenvalue for each identified position can be determined using a shifted symmetric higher-order power method. Since the symmetric higher- order power method is an iterative eigen decomposition that only allows for local searches, its precision is highly related with the starting point, which may be taken randomly, to ensure that the whole space of rotations is covered. The actual optimum is always found because this process is repeated many times, for example using 1000 uniformly distributed quaternions as starting points.

[0133] This technique can be used to identify the best rotations and the best peak positions, which may suggest step 212 is unnecessary. However, determining eigenvalues and eigenvectors is computationally expensive and so very time consuming if applied to all voxels in a typical tomogram 10. To avoid this problem, the much quicker technique of using the Frobenius norm is first used to determine the peak positions (this technique cannot identify the rotations) at steps 210 and 212. The Frobenius norm is related to its largest eigenvalue and so can be used as an approximation for finding the peak positions. Hence the Frobenius norm approximates the correlation value at each voxel, from which the peak positions are determined that correspond to the target molecule locations in the tomogram 10. Hence, steps 210 and 212 use the Frobenius norm and peak finding to identify peak positions that are then fed into the more refined technique of using eigenvectors to determine the best rotations of the target molecule at step 214. This produces much quicker execution times.

[0134] Alternatively, to ensure that the optimal rotation is not missed during the computation of the eigenvector, the space of rotations of the target molecule expressed as quaternions can be sampled. Afterwards, the inner product between the selected tensors of the convolved tensor field and the quaternions is computed so rotations can be ranked by this value. Finally, a list of peak positions with their most promising rotations (highest inner products) per peak is returned. Despite the necessity of sampling the space of rotations, this process is faster than computing standard cross-correlations because computing the inner product between a tensor and a quaternion is less costly than computing conventional local normalised cross-correlations. Moreover, this computation is suitable for further acceleration using high performance computing techniques.

[0135] As noted above, a symmetric positive semidefinite linear operator may be used when calculating the tensorial template 222 at step 204. This ensures that the inner product between the tensors and the unit quaternions expanded to n-order tensors attain the global maximum for the optimal template rotation. This global maximum is found in step 214 either by sampling the space of rotations defined by the quaternions or by unfolding tensors of the convolved tensor field 226 into matrices and finding their eigenvalues and eigenvectors. Thus, use of the symmetric positive semidefinite linear operator ensures the necessary information to provide correlation values that reflect the closeness of the match of the tensorial template 222 at any particular voxel in the tomogram 10.

[0136] Step 214 can produce a list of potential peaks as a peak position and best rotation (or several candidates providing the best rotations, as described in the preceding paragraph) against a correlation value, as shown at 215 in Figure 3. However, a consequence of using tensors is that the correlation values are not normalised across the entire tomogram 10. There is likely to be a variation across the tomogram 10, for example due to variations in illumination for various positions within the tomogram 10. A lack of normalised correlation values limits any attempt to rank the potential peaks. In addition, when the list contains more than one rotation per peak, providing normalised correlation values is required to allow a meaningful ranking of the best correlation values which in turn allows the best rotation per peak to be identified and selected.

[0137] Step 216 provides normalised cross-correlation values for each potential peak by utilising the conventional technique of local normalised cross-correlation, albeit in an extremely limited way. Advantageously, in this method 200, the calculation is required only for each potential peak position and best rotation value (or for each peak position and the limited list of best rotations, where provided). That is, a local normalised cross-correlation calculation is performed by comparing a conventional template 20 of the target molecule against the tomogram 10: however, it is performed by comparing the template 20 of the target molecule at the position for the target molecule found in step 212 and for only the single rotation of the target molecule corresponding to the best rotation found at step 214 (or for the select few best rotations identified). This may be contrasted to prior art techniques where local normalised cross- correlation is calculated for every single rotation of the target molecule (which is often around 30,000 rotations) at every single voxel. Hence, computational times are very much reduced and / or a larger number of best rotations may be used to increase accuracy.

[0138] Consequently, at step 216, a conventional template 20 is generated from the model of the target molecule: the template 20 provides a snapshot of the target molecule at the identified best rotation (or a template 20 for each of the best rotations where more than one rotation is provided). The local normalised cross-correlation value is then calculated for the template(s) at the identified position. Where more than one rotation is provided in step 214, the cross-correlation of step 216 provides a normalised cross- correlation score for each rotation. This allows the cross-correlation score of each of the best rotations to be compared in a meaningful way: the rotation achieving the highest cross-correlation score may then be selected and provided with the identified position.

[0139] At step 218, a list 219 of potential peak positions is provided. This list 219 contains, for each instance of the target molecule, the identified position, best rotation and normalised cross-correlation value (although, in some embodiments, the best rotation and / or normalised cross-correlation value may be omitted). The list 219 may be provided in different ways and with or without user input. For example, the list 219 may be ranked according to the normalised cross-correlation values.

[0140] Also, the numbers of potential peaks returned in the list 219 may be set by a user, as indicated at step 217. Further parameters set by the user at step 217 may include an absolute peak threshold, a relative peak threshold, and a minimal distance between peaks.

[0141] The number of potential peaks returned in the list 219 set by a user vary according to the intended use. As a first example, all instances of a particular target molecule may be wanted, in which case all potential peaks may be returned. Alternatively, to limit false identifications, a threshold normalised cross- correlation value may be set such that only potential peaks meeting or exceeding this cross-correlation value are included in the list of potential peaks 219. As a second example, a user may want to refine the structure of a target molecule through averaging the images of aligned instances of the target molecule found. This may be done by aligning each instance of the target molecule using the best rotation values and then averaging the aligned images of the target molecule. In this example, it is likely to be better to use fewer potential peaks, namely those peak positions with the best normalised cross-correlation scores as these instances will correspond best to the already determined target molecule structure (hence allowing that determined structure to be refined further).

[0142] An alternative method 500 is shown in Figure 5 that uses block processing to reduce memory requirements. A tomogram 10 is divided into three-dimensional blocks 502 that are processed separately, before the results are merged. In the embodiment shown in Figure 5, the tomogram 10 is cubic (although illustrated as a square for simplicity) and is divided into eight contiguous cubic blocks. As will be appreciated, other arrangements may be used depending on the size and shape of the tomogram 10 and the memory available. Figure 5 shows one of the blocks 502, again as a two-dimensional square within the two-dimensional visualisation of the figure. The method 500 corresponds to method 200 already described, although implemented on a block-by-block basis. Hence, only the differences will now be explained.

[0143] Step 208 sees the tensorial template 222 convolved with part of the tomogram 10 to identify peaks, and hence molecules, within one of the blocks 502. To account for the fact that one or more target molecules may straddle a border between adjacent blocks 502, the tensorial template 222 is convolved with a portion 505 of the tomogram 10 that includes the block 502 and a halo 504 around the borders of the block 502 to overlap with adjacent blocks 502. Hence, the halo 504 extends into the adjacent blocks 502. The width of the halo 504 is selected to be around the greatest dimension of the target molecule (i.e., height, width, depth), for example between half of and all the greatest dimension as this ensures any instance of a target molecule is enclosed within one block 502 and its halo 504.

[0144] The portion 505 is then analysed in the same way as was described for method 200. The convolution at step 208 produces a convolved tensor field 226 that is processed at step 210, for example using the Frobenius norm, to produce a scalar map 228. Step 212 sees the peaks identified according to the local maxima in the correlation values present in the scalar map 228. Only peaks within the block 502 are identified or, if peaks in both the block 502 and the halo 504 are identified, only the peaks in the block 502 are retained and peaks in the halo 504 are discarded (as they will be identified when the block 502 in which the halo 504 resides is processed). At step 214, the best rotations are determined, for example using eigenvectors and eigenvalues found for each peak position in the list of peaks to determine the best rotations at each peak position, or alternatively to find several rotations for each peak using the inner products as described above.

[0145] Steps 208, 210, 212 and 214 are repeated for all the 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, conventional localised cross-correlation is used to determine the normalised cross-correlation values and these values are computed only for each peak position and best rotation pair (or each rotation for each position where more than single rotation is identified in step 214, in which case the best rotation is also identified from the normalised correlation values). The refined list of peak positions, rotations and normalised cross correlation scores are provided in step 218, as may be guided through user input of a desired number of peaks at step 215.

[0146] A person skilled in the art will appreciate that the above embodiments may be varied in many different respects without departing from the scope of the present invention that is defined by the appended claims.

[0147] Although the above embodiments relate to use of the invention in Cryo-ET, the present invention has many other uses. For example, the present invention can be applied both to 2D images and 3D images. The present invention may find use in object detection in image processing, counting objects in metrology, fiducial tracking for drift correction, fiducial detection for tomogram reconstruction, defects detection in semiconductors, etc. Even within Cryo-ET, the present invention can be used both for protein structures discovery and functional cellular behaviour studies.

[0148] There follows a mathematical presentation of the concepts and proofs that demonstrate the efficacy of the methods described above. 1 Introduction

[0149] A classical problem in image processing and, particularly, in pattern recogni- tion, is to identify if a large image contains copies -and how many, and their locations and orientations- of a small image, named “template” . The resulting algorithms arc generically known as template matching algorithms [1, 2, 3■] . The most classical solution is based on using cross-correlations, although there arc other approaches based, for example, in mctahcuristic algorithms [i] or on deep learning [ 1 , 5, 6]. In this document, we show the mathematical foundations of the cross-corrclation-bascd template matching algorithm (TM in all that fol- lows), and we introduce a new fast algorithm that solves the problem using tensors.

[0150] The main advantages of TM, when compared to the algorithms based on machine learning, arc that TM is a white box model, it is directly applicable when you have just one template and one larger image (not requiring any kind of training, which may be a very difficult task in some applications), and locates rotations with arbitrary precision (Current deep learning based algorithms for template matching in three-dimensional images arc not able to estimate rota- tions accurately [4]).

[0151] On the other hand, a major drawback of TM is its computational cost. TM basic idea is to compute the inner product between the (rotated) template and the (translated) image, and normalize the result. These computations arc made, for each rotation, in the Fourier domain to efficiently address translations [7, 3, •;]. However, this process has to be repeated for every rotation to be investigated, thus the resulting complexity has a dependency with the rotations processed. The computational cost of this process may become restrictive for 3D images since SO(3), the space of rotations of R3. is a (compact) manifold of dimension 3. In application domains such as cryo-clcctron microscopy, there arc required more than ten thousand rotations for achieving an angular precision of a few degrees.

[0152] We propose an algorithm called tensorial template matching, TTM, which integrates into a unique symmetric tensor the information relative to the tem- plate in all rotations. In other words, the tensor template incorporates in a unique object the information about all rotations of the template, thus allowing us to find the position and rotation of instances of the template in any tomo- gram with just a few correlations with the linearly independent components of the tensor. The tensor template is computed only once per template, and, as soon as it is generated, it enables to process any image.

[0153] 2 Classical template matching

[0154] Let us introduce some notation, d-dimcnsional images arc just elements of L2(Hd), which is a Hilbert space with the inner product (ƒ, g) = ƒRdƒ(x)g(x)dx. It is natural to use the inner product to compare two images ƒ, g of the same size. Concretely, we can use that . where 6 is the angle formed by ƒ and g. In particular, ƒ = αg for some positive constant

[0155] Template matching is typically used to study if instances of a “small” image t (the template) is present in a larger image ƒ, c.g. look for instances of an specific macromolecule in a cryo-clcctron tomogram (3D volumetric image) . The size of the image is connected to the set of points where the image docs not vanish, the support of the image. That is, t is meant “small” when the set is small (c.g., is a subset of a small ball D) . Let’s assume that f and t have quite different sizes, so our interest is to compare t (the template, the small image) with just a part of f. In such ease we need to introduce some special operators S ( ) ) that fix our attention in just a part of the domain of f. An interesting example of such operators is where denotes Heaviside’s unit step function and . If the support of the template the ball of radius r centered at the normalized inner product informs about the similarity between t and the restriction of f to ( ) More- over, if we introduce the translation operator and compute the result informs about the similarity between t and the restriction of f to . Of course, it may happen that f contains a copy of a rotated version of t, so rotations arc also necessary for a complete discussion of the problem Thus given we define the operator we define a rotated version of t,

[0156] The normalized inner product informs about the similarity of tp and the restriction of f to Dx(r). It is important to notice that

[0157] Thc operator Srdefined by (1) has some special properties. Concretely, it is symmetric, semidefinite positive and commutes with rotations Recall that, given (X, (•, -)x) a (real) Hilbert space (in all that follows, we will also use the notation x • y to denote (x, y), if this simplifies computations.), an operator S : X — > X is named:

[0158] • Symmetric (also named self-adjoint) if

[0159] • Semidefinite positive, if • Definite positive, if it is semidefinite positive and = 0 implies f = o.

[0160] If S : X — > X is a symmetric semidefinite positive operator (SSP, in all what follows), then X becomes a pre- Hilbert space with the inner product ff, g)s ■= ff, S(g))x (3) and the seminorm llƒlls = VWJTs (4)

[0161] Observe, for example, that if S is given by (1), then ||ƒ||s = 0 means that ƒ|D0(r) = 0 almost everywhere.

[0162] Theorem 2.1 and consider t

[0163] Moreover, if f, g G L2(Rd), || g\\s yt Q, the following are equivalent statements:

[0164] W (f, g)s = \\f\\s\\g\\s -

[0165] Proof. As S is SSP, we have that, for all a e R,

[0166] 0 < (ƒ + ag, f + ag)s = llƒlls + Mƒ, g}s + a2||.g||||. (5)

[0167] Hence, if ||.g||s 0, the only way that the quadratic polynomial (in a) above is nonnegative everywhere is that

[0168] Kf, g)2s — 4||ƒ||||||.?|||| < o, which is equivalent to

[0169] |(ƒ, .?>s| < ||ƒ||s||.?||s- (6)

[0170] On the other hand, if ||.g||s = 0, the only way to satisfy (5) is (f, g}s = 0, in whose ease (6) also holds. This proves (a).

[0171] Let us now demonstrate (b) O (c) whenever 0. Indeed, (c) is equiv- alent to

[0172] 0 which holds if and only if tf, g)s = \\f\\s\\g\\s- Thus (&) •<=> (c). >

[0173] Note that, if ƒ, ( ) arc two images, a > 0, and we take S = Srgiven by ( i ), then means that ƒ has a match with tRin the unit ball centered at x. Indeed, there arc many ways to define operators S with the property that \\ means that f = g in a neighbourhood of 0, so that 0 means that f has a match with a rotated version of t in a neighbourhood of 0. Although arbitrary SSP operators may not enjoy this property, they allow the creation of a general way to deal with this kind of operators.

[0174] Thus, in all what follows, we assume that is an SSP operator and that || ||| , where 1 (T) = 1 is the constant image.

[0175] Notice that 1 is not an element of ), but this can be managed in several ways. In fact, although we use to denote the space of d-dimcnsional images, in practice we only consider images f with compact support K. Then, when we compute (ƒ, 1), we mean ƒR K- Moreover, since our interest is on operators S that vanish on functions vanishing outside of a certain neighbourhood D of 0, by (1, 1)s we mean àIXD , 1%D)S -

[0176] Rotations and composition of operators play an important role in embod- iments of the present invention set forth therein. Thus, it is natural to ask how the composition of rotations acts on the images. This is, indeed, a simple computation:

[0177] Hence (7) and

[0178] (8)

[0179] Given an image ƒ, we consider its projection onto the space of images which arc S— orthogonal to the constant image 1,

[0180] Remark 2.2 These projections are important to study invariant properties with respect to constant brightness changes in the images of translations and rota- tions. Note that there is no “real” difference between an image f and the images of the form f , When we modify the constant a, what we observe is a uniform change in the density or the brightness, but not the apparition of new structures or forms, in the image f . Thus, f and its projection Ps(ƒ) essentially represent the very same image, since for certain

[0181] Given two images ƒ, t, we have that Hence since Consequently, if and , there arc two constants such that

[0182] Assume that S commutes with rotations, and take fixed. Then, for each we have that since . Moreover

[0183] (just take

[0184] Hence and this, in conjunction with (10), implies that

[0185] Hence is an S-orthogonal decomposition of tR, which means that the constant β that multiplies 1 in the -orthogonal decomposition of tRdocs not depend on R, and

[0186] In particular, for each x G EC the problems: • Maximize over rotations R.

[0187] • Maximize over rotations R.

[0188] • Maximize over rotations R. arc equivalent.

[0189] Let us define:

[0190] Lemma 2.3 If S is an SSP operator that commutes with rotations, the param- eter 5 that appears in the S-orthogonal decomposition

[0191] • Maximize over rotations R.

[0192] • Maximize over rotations R. are equivalent.

[0193] Proof. Wc know that so that wc only need to prove that docs not depend on R. Indeed,

[0194] (Make the change of variable z = Ry )

[0195] (since S commutes with

[0196] Wc can now state and demonstrate the following:

[0197] Theorem 2.4 (Classical template matching) Let S be a SSP operator which commutes with rotations and let be fixed. Then the following are equiv- alent problems:

[0198] (a) Maximize over rotations R. Maximize over rotations R.

[0199] Maximize overrotations R.

[0200] (d) Maximize over rotations R.

[0201] Moreover, if and S also has the property that for a certain neighborhood which contains the supports of all the rotated templates t,Q with then a match between f and tRin x is got whenever any one of the following claims hold:

[0202] Finally, the normalized correlations described in (a*), (b*), (c*), and (d*) do not change when we substitute

[0203] Proof. The equivalences and ( ) ( ) have been already shown. The following identities demonstrate

[0204] ( just take Rz = y and use that det R = 1 )

[0205] (since S commutes with OR-1 )

[0206] The other claims arc a direct consequence of Theorem 2.1.

[0207] □

[0208] In all that follows, wc assume that S is an SSP operator that commutes with rotations and t is normalized in the sense that Then forall rotation R. Consequently, and attains its maximum (= 1) if and only if there is a perfect match between ƒ and tRin x. Moreover, if we define consider the cross-correlation of functions ƒ ( ) which is defined by then

[0209] A perfect match is, in general terms, never attained. This is so because the desired image, represented by the template t, is usually supported on a strict subset fl of the domain D were the operator S is able to distinguish functions. Thus, the image f may well contain a copy of the image represented by tp but in the neighbourhoods of the support of tp, f will contain some information which is not present in tp. In addition, f is usually corrupted by noise and distortions. This means that the normalized correlations described in items (a*) — (d*) of Theorem 2.4, will never equal 1. Consequently, a threshold should be introduced in order to decide if a match has (or has not) been produced.

[0210] In order to find the rotation which maximizes c(x, R), the cross-correlation ƒ should be computed for a huge amount of rotations R, which makes classical matching an inefficient approach for template matching. Indeed, for d = 3, the size of the set of rotations R used to sample SO(3) well enough to guarantee a reliable result varies between 104and 107rotations.

[0211] Due to numerical reasons, high frequencies may be altered during rotation transformation. Thus, in practice, we do not apply the operator S to the original images ƒ, t but to a filtered version of them that eliminates these high frequen- cies. Concretely, we apply an isotropic (i.c. rotation invariant) low-pass filter h to both images and, after that, we apply the template matching algorithm to the resulting images. The idea behind this is that, if there is a match between ƒ and t, there will be a match between f f h too. The operator S results from applying a rotationally symmetric mask m ) p(|| ||) to the given image. Thus, we substitute Then we apply the classical (or tensor) matching algorithm to the pair of images f, t using the SSP operator (f) ( ) ( )f( ) Usually, the mask m equals 1 within a certain radius around 0 and equals 0 outside a sightly larger radius. In between these radii the mask takes values between 0 and 1. Under these restrictions, it is clear that the operator S is SSP and commutes with rotations. Moreover, if , we have that ƒ where D is a ball of positive radius centered at 0. Let us compute the inner product

[0212] (since every filter is translation invariant, and h is isotropic) where and we use • to denote the standard product of real functions. This means that we would have the same effect just considering the template matching algorithm associated with the operator applied to the images f, t. Moreover, the following holds:

[0213] Lemma 2.5 (14) with h defining an isotropic filter and m a rotationally symmetric mask as de- scribed above. Then S is SSP.

[0214] Proof. For the proof, we use the following (well-known) formulae: For functions , we have that and

[0215] Let us now consider the product

[0216] (since ). This proves that S is semidefinite positive. Let us show the symmetry:

[0217] (since g * h = h * g and • is the standard product of functions)

[0218] □

[0219] In all that follows, we assume that the SSP operator S is of the form (14) with h, m verifying the hypotheses of Lemma 2.5. Thus, the template matching algorithm is applied with this operator and a fast computation of c(x, R) is needed. A direct computation leads to:

[0220] (since h is isotropic, and m is rotationally symmetric). Moreover,

[0221] Now, using the definition of S (and imposing h * 1 = 1), we can simplify the computation as follows: Note that the FFT algorithm can be used to compute the inner products ap- pearing at the end of the formula above, which helps to fasten the algorithm. Indeed, if ƒ, g arc two images, , so that (15)

[0222] Moreover, the following identities also hold: and

[0223] Thus,

[0224] Algorithm 1 , devoted to classical template matching, has been developed using the formulae above.

[0225] An important tool we will use in the embodiments of the present invention set forth herein is the set of quaternions. In particular, we will use that rotations can be parametrized by unit quaternions (which can be identified with the unit 3-sphcrc S3) , as well as the following formulae (see, c.g. [ 10, 1 1 ]):

[0226] •

[0227] • Given , we identify x with a pair (a, v) where ( ) , and call a, the real part of . Then, if we have that

[0228] Consequently, if have norm 1, then

[0229] We end this section with a well-known result about composition of SSP operators that will be used in the proof of the main theorem described here:

[0230] Lemma 2.6 are semidefinite positive symmetric and commute, then TS is also semidefinite positive symmetric.

[0231] Proof. If S, T arc SSP, then they have (unique) positive square roots which arc symmetric operators too. Moreover, the composition is also symmetric. Since T, S commute, their positive square roots commute. Hence and the square of any symmetric operator is positive. > Algorithm 1 Classical template matching with rotations

[0232] Load and Fourier transform of the image

[0233] Load data into f

[0234] Load mask and pre-process the image

[0235] Load the mask into m

[0236] Load and pre-process the template

[0237] Load template into t

[0238] Compute the output 3 Tensor template matching

[0239] This section introduces a tensorial template matching (TTM) algorithm. The purpose is to handle translations and rotations efficiently at the same time. First, we introduce some background necessary to understand further mathe- matical developments. Second, we present the main theorem for TTM, which allows us to determine the optimal rotation of the template, t, on every match in the image f without sampling the SO(3) by computing some tensors. Finally, we explain how to determine match positions (template translations) directly from the computed tensors.

[0240] 3.1 Tensor background

[0241] A tensor of order n and dimension d is just an array of the form where all the entries arc real numbers. The tensor A is named symmetric if ( )every permutation (the set of permutations of {I, - - - , n}). Wc denote by ( ) the set of symmetric tensors of order n and dimension d. An important example of symmetric tensor of order n is the so called n-th tensor power of a vector , which is defined as

[0242] It is well known that ) ( ) arc real vector spaces with the natural operations (pointwisc sum and multiplication by a scalar), and that dim For example, dim . Moreover, every symmet- ric 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 minimal number of tensor powers used to represent the tensor with their sum

[0012] .

[0243] The map (•, •) given by defines an inner product. It is also usual to denote A • B = (A, B). Moreover, with this notation , , y arc d-dimcnsional vectors, a direct application of the multinomial theorem shows that

[0244] Moreover, if ( ) we can also consider the inner product which can be seen as an homogeneous polynomial in d variables, of degree n, which justifies using the notation Moreover, if k < n, ( ) denotes the symmetric tensor whose components arc

[0245] In particular, is a vector whose z-th component is

[0246] Indeed, if then where denotes the gradient of the function

[0247] Note that the vector x can be chosen from Cdin the definitions above. This justifies the following definition (see

[0013] ): Given ) ) We say that is a B-cigcnvaluc of A and is its associated eigenvector (equivalently, that (A, u) is a B-cigcnpair of A) if 1.

[0248] Using gradients, we can rewrite the equation as

[0249] Hence u is a B-cigcnvcctor of A if and only if it is a critical point of the following optimization problem:

[0250] Two particularly important eases arc the B-cigcnvcctors and the Z-cigcnvcctors

[0251] The optimization problem associated with finding Z-cigcnvcctors of a given symmetric tensor is particularly important for us since the tensor matching algorithm we propose is reduced to one of these problems in each position, and, fortunately, there arc good iterative algorithms to approximate the solutions of (25) (see c.g. [14 , 15]) . These algorithms have a linear rate of convergence. In section 3.4 we show an heuristics that can be used to select the positions where a match is probable, so that solving (25) is necessary. 3.2 Defining of the Tensor template

[0252] In all that follows, denotes an SSP operator which com- mutes with rotations and the template is assumed to be normalized by In section 2 we proved that, for each attains its maximum value on rotation R (and this value equals 1) if and only if there is a match between f and t at (i.c., a match between Let us define the symmetric tensor n( ) ( ), where d! is the number of parameters used to describe the rotations SO(d) (in particular, for d = 3, we get d' = 4) by the formula:

[0253] This means that for all Hence where (28) is a tensor template (or tensorial needle).

[0254] It is of fundamental importance to observe that T(z) is computed only once and contains a reduced number of components, since dim (in particular, dim . Indeed, this is one of the main rea- sons why the tensor template matching algorithm we introduce is fast. The other reason is that rotations R defining a match between f and tRat x arc Z-cigcnvcctors of the symmetric tensor Cn(x), which is really remarkable be- cause the power method used in [ 15, 1 4] for the solution of the corresponding optimization problem docs not require using myriads or even millions of rota- tions -as is the ease with classical matching algorithms- but just a reduced set of them: one by iteration. 3.3 Finding the correct rotation

[0255] Let us state the main result of this document:

[0256] Theorem 3.1 Let f, ( ), , be given. If there is a match between f and tRat x, the function ( ) ( ) , defined on rotations of . when parametrized by unit quaternions Q, attains its global maximum (and it is equal to 1) at Q = R.

[0257] Proof. Wo prove the result as a consequence of Theorem 2.4. Thus, our main goal is to represent < in terms of a scalar product for some SSP operator S , which would guarantee that if there is a match between f and tRat x, then p( ) attains its global maximum (and is equal to 1) at

[0258] Let us compute

[0259] (by definition of c(x, Rf) . Hence, dividing by w(x), we get:

[0260] (where S'(t)(z)(R) = (OR-1 ° )(t)(z)) where K(R) = (Re(R))nand denotes the convolution of functions

[0261] In other words, (29)

[0262] Here, must be interpreted as a function defined on with values on (indeed, it is an element of where

[0263] Indeed, in general, every element can be interpreted, for each , as an element of L2(SO(3)) just making t(z)(R) = tp(z) . Consequently, is an element of L2(R3),

[0264] It follows that

[0265] (in the last equality, set R = QP and use that ) where Iddenotes the identity rotation

[0266] Let us now denote by the operator given by and let

[0267] (since S, S2commute with rotations)

[0268] Thus, the proof ends as soon as we demonstrate that S' is an SSP operator, and it is for this that we need to use Lemma 2.6. Indeed, S' = S2o S is a com- position of operators, S is, by hypothesis, symmetric semidefinite positive, and S, S2commute because S commutes with rotations and S2is defined in terms of convolution in SO(3). Thus, Lemma 2.6 implies that S' is SSP whenever S2is SSP.

[0269] To prove that S2is symmetric semidefinite positive, we use the properties of the convolution on SO(3) when interpreted as a hyperspherieal convolution on S3, the unit sphere of 4. Recall that if denotes the (unit) sphere of R , then SO(d) acts transitively on (which means that, given i,2there is a rotation such that which makes of a homogeneous space and allows to introduce the convolution of functions defined on S as follows: where is the north pole of the sphere and ƒ Thus, if we use that the elements of SO (3) arc parametrized by quaternions of norm 1, which that can be identified with the elements of the sphere in fourth-dimensional space, then assuming that the north pole of S3is given precisely by the identity rotation Id, the convolution of can be interpreted as a hyperspherieal convolution on S3:

[0270] Now, as it is well known, ( ) is a Hilbert space and the so-called hyper- spherieal harmonics, . form an orthonormal basis of this space. Thus, every function admits a Fourier expansion Moreover, in

[0016] , it was proven that, if , then

[0271] It follows that, given a template t, for each , the map ( )( ) ) belongs to L2(S3) (here the rotations R arc parametrized as unit quaternions, so that and and

[0272] We need the following Lemma, whose proof is included in section 4:

[0273] Lemma 3.2 0 for all t.

[0274] Then

[0275] This means that convolution with K, which is an operator ( ( )) , is semidefinite positive. Moreover, it is well known that this operator is symmetric (and we will use both things in our computations bellow).

[0276] In order to prove that S2 is SSP, we introduce the operator ( ( ), ( )) defined by ( )( )R, as well as the operator ( ( ) ( )) defined

[0277] Then where ƒ On the other hand,

[0278] Thus, the volume of 50(3), then where we have used that and that so that It follows that

[0279] Thus, S2(f) is scmidcfinitc positive.

[0280] Moreover, the same type of computation shows that which proves that S2 is symmetric. This ends the proof of the theorem. Note that Theorem 3. 1 connects the problem of finding, at a given position x, the rotation R which gives a match between f and tn;at x with the problem of finding the dominant Z-cigcnvaluc-cigcnvcctor pair by solving ( and n e ven.

[0281] 3.4 Finding the correct position

[0282] In the previous subsection we showed how the tensor-based correlation function can be seen as using a slightly different degenerate inner product (based on S' from the proof of Theorem 3.1 , rather than S). Unfortunately, this was not taken into account in the normalization: not for the template t, but even more crucially, also not for the image f. So what is the impact of this? First of all, observe that the operation that is missing from S in the normalization is effectively a kind of convolution, and that it should at most scale the constant component of an image. As a result, if an image is S-orthogonal to 1, it will also be S'-orthogonal to 1. However, norms will be affected.

[0283] For the template t this means that the normalization is off by a certain factor, but this factor is the same everywhere, so it docs not affect the ratios between responses at different positions. For the image ƒ, the impact is less benign though, as g will differ from 11 ( ) 11 g, in a nonuniform way.

[0284] So when will this shortcoming cause a false positive? For this to happen, the normalization factor used at a non-match position would have to be much higher than the “correct” normalization factor, andƒor the normalization factor would have to be (way) too low at a match position. Since the difference between S and S' is essentially a smoothing operation, and the normalization factor is the reciprocal of the norm of the projected image, the image would thus have to be (very) smooth at the non-match position, while exhibiting a lot of high frequency energy around the match position. Such a situation would not be impossible, but would at the very least be unusual in the context of a typical application like the analysis of electron microscopy images.

[0285] Although we can find the spatial positions of peaks by running an algorithm to find the dominant Z-cigcnvaluc-cigcnvcctor pair for each and every voxel, this is fairly expensive using the current decomposition algorithms for higher degree tensors. However, the Frobenius norm of a tensor is related to its spectral norm, and in practice it turns out it can be used as an excellent proxy for finding the spatial locations of peaks. Indeed, we know that ) ( ) Now, if denotes the spectral norm of tensor || || denotes its Frobenius norm, it is well-known that the largest singular value of T equals its spectral norm, and that

[0286] (scc c.g.,

[1718] ).

[0287] In fact, the connection between is stronger than just this inequality. As is well-known, every tensor is a finite sum of tensors of rank 1 (indeed, if the tensor is symmetric, the tensors of rank one can also be chosen symmetric)

[0012] . Moreover, if is a tensor of rank 1 satisfying ( ) then (see, c.g., [li]) and

[0288] Thus,

[0289] \

[0290] Hence if Ei(T) is preserved, an increment on the size of || ||CT(|| || p, respec- tively) is translated into an increment on the size of , respectively) .

[0291] Moreover, in 1938 in

[0020] was demonstrated that, for any symmetric tensor T,

[0292] Thus, large implies large spectral norm of T, and the spectral norm of ) is strongly connected to the optimization problem solved in Theorem 3. i , which justifies using the Frobenius norm of as a parameter to select positions x where a match is possible.

[0293] For each position x identified as a potential peak, the SS-HOPM algorithm (see [14, 13, 19)] for precise definition and implementation of this algorithm) is used to find the exact dominant Z-cigcnvaluc and its associated Z-cigcnvector. which is the rotation R candidate to give a match at x. 3.5 Integrating the template

[0294] Properly integrating the noodle over all rotations is a critical part of the proce- dure. In this work, this is accomplished by summing over a largo sampling sot of rotations. In the traditional approach, the exact distribution of rotations is not terribly critical, as wc only retain the maximum response. However, when using the samples for integration, it is imperative that they arc as uniformly dis- tributed 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 (hypcr)sphcrc (except for a circle). Here, a two-step approach is used: first, the approach suggested by [ ; : ] is used (with some tweaks) to generate points on a sphere rep- resenting rotation axes, and then, those axes arc combined with (ID) rotations distributed in such a way that the overall distribution on the 3-sphcrc is still uniform.

[0295] 4 Proof of Lemma 3.2

[0296] Let us start recalling the formulae associated to Fourier expansions in hypcr- sphcrical harmonics on the sphere S3. The parametrization of the sphere wc consider is the following one: where is identified with the unit quaternion Q j , which represents a rotation of three dimensional euclidean space . The volume element (used for integration on S3and, henceforth, also in SO(3)) is then given by f f p

[0297] Then every function can be decomposed as where denotes the orthonormal basis of L (S ) formed by the hy- pcrsphcrical harmonics and arc the Fourier coef- ficients of f in this basis. Wc want to prove that 0 for all . Now, and where ^ is a positive constant and denotes the Gegenbauer polynomial of degree £, which appears as the £-th Taylor coefficient in the expansion: (1 — It is well-known that Chebyshev’s polynomial of second kind) and that Hence

[0298]

[0299] To estimate the integral above, we need to use a few trigonometric formulas, as well as the hypothesis that n is even. Concretely, n even implies that nƒ2 is an integer and . Moreover, for odd, we have that

[0300] This makes the integral equal to 0 for

[0301] Assume Then (for the last line, just set k = n / 2 — .s). Moreover, it is well-known that so that, by a direct substitution in the formula defining we get

[0302] We can now use that to claim that

[0303] The parity of and n implies that all factors that appear multiplying the variable 0 inside of the eosine functions arc even numbers. This makes the corresponding integrals (on [0, TT]) equal to 0, except in the ease that the factor itself is 0. In such ease, cos(0) = 1 implies that only the eosine functions that appear with a minus sign in front of them in the formula can contribute with a negative number to the integral. Now clearly always since , 0 implies , so that / which is impossible since the sum’s range goes from k = 1 to k = n / 2. This means that the term — cos never contributes with a negative number to the sum. On the other hand, if the eosine function with factor + + contributes, which means that , ) / / In particular, taking we have that / so that cos(0(£ n)) = cos ( ( )) ( ) 1 and the corresponding term effectively appears in the sum. In particular, adding these two terms of the sum we get since n even and This ends the proof of Lemma 3.2

[0304] □

[0305] 5 Conclusions

[0306] We have exposed the maths of classical template matching with rotations. More- over, an alternative to the classical algorithm, named tensorial template match- ing (or TTM), has been shown. TTM integrates the information relative to all rotated versions of a template t into a unique symmetric tensor template T, which is computed only once per template. The main theorem, Theorem 3.1, shows that finding an exact match between an image f and a rotated version tn °f the template t at a given position x is equivalent to finding a best rank 1 approximation (in the Frobenius norm) to a certain tensor . The resulting algorithm has reduced computational complexity when compared to the classical one. TTM finds the position and rotation of instances of the template in any to- mogram with just a few correlations with the linearly independent components of T. In particular, Cryo-clcctron tomography (3D images) for macromolecular detection requires 7112, 45123 and 553680 rotations to achieve an accuracy of 13°, 7° and 3° respectively

[0022] . Therefore, and considering 4-dcgrcc tensors (35 linearly independent components) , the potential speed-up of our approach with respect to TM is 203x, 1239x and 184560x in these eases.

[0307] References

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

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

[0310] [3] R. Gonzalez and R. Woods. Digital image processing, / Th Global Edition. Pearson Education, 2017.

[0311] [4] G. Corona, O. Maciel-Castillo, J. Morales-Castaneda, A. Gonzalez, and E. Cuevas. A new method to solve rotated template matching using meta- hcuristic algorithms and the structural similarity index. Mathematics and Computers in Simulation (MATCOM), 206(C):130-146, 2023. [5] L. Lamm, R.D. Righctto, W. Wictrzynski, M. Pogo, A. Martinez-Sanchez, T. Pong, and B.D. Engel. Membrain: A deep learning-aided pipeline for detection of membrane proteins in cryo-clcctron tomograms. Computer Methods and Programs in Biomedicine, 224:106990, 2022.

[0312] [6] E. Mocbcl, A. Martinez-Sanchez, L. Lamm, R.D. Righctto, W. Wictrzyn- ski, S. Albert, D. Lariviere, E. Fourmentin, S. Pfeffer, J. Ortiz, W. Baumcis- tcr, T. Peng, B.D. Engel, and C. Kcrvrann. Deep learning improves macro- molecule identification in 3d cellular cryo-clcctron tomograms. Nature Methods, 18(11) : 1386— 1394, 2021.

[0313] [7] J.P. Lewis. Fast template matching. In Vision interface, volume 95, pages 15-19. Quebec City, QC, Canada, 1995.

[0314] [8] J. Bohm, A.Ss Frangakis, R. Hcgcrl, S. Nickell, D. Typke, and W. Baumcis- tcr. Toward detecting and identifying macromolecules in a cellular context: template matching applied to electron tomograms. Proceedings of the Na- tional Academy of Sciences, 97(26):14245-14250, 2000.

[0315] [9] A.M. Roseman. Particle finding in electron micrographs using a fast local correlation algorithm. Ultramicroscopy, 94(3-4) :225-236, 2003.

[0316]

[0010] H. D. Ebbinghaus, H. Hermes, F. Hirzcbruch, M. Kocchcr, K. Mainzer, J. Ncukirch, A. Prestel, and R. Remmert. Numbers. Springer, 1991.

[0317]

[0011] Lev Pontryagin. Generalization of numbers. CrcatcSpacc, 2010.

[0318]

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

[0319]

[0013] C.F. Cui, Y.-H. Dai, and J. Nic. All real eigenvalues of symmetric tensors. SIAM Journal on Matrix Analysis and Applications, 35(4) :1582- 1601 , 2014.

[0320]

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

[0321]

[0015] T.G. Kolda and J. Mayo. Shifted power method for computing tensor cigcnpairs. SIAM J. Matrix Anal. Appl., 32:1095-1124, 2010.

[0322]

[0016] I. Dokmanic and D. Pctrinovic. Convolution on the n-sphcrc with applica- tion to pdf modeling. IEEE transactions on signal processing, 58(3) : 1157— 1170, 2009.

[0323]

[0017] 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.

[0324]

[0018] K. Kozhasov and J. Tonclli-Cucto. Probabilistic bounds on best rank-one approximation ratio. Arxiv, 2022.

[0325]

[0019] P.A. Regalia and E. Kofidis. The higher-order power method revisited: convergence proofs and effective initialization. In 2000 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings ( Cat. No.00CH37100), volume 5, pages 2709-2712 vol.5, 2000.

[0020] S. Banach. Uber homogcnc polynomc in (L2). Studio, Mathematica, 7(1):36- 44, 1938.

[0326]

[0021] C.G. Koay. A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere. J Comput Sei., 2(4):377- 381, 2011.

[0327]

[0022] M.L. Chaillot, G. van dor Schot, I. Gubins, S. Root, R.C. Vcltkamp, and F. Forster. Extensive angular sampling enables the sensitive localization of macromolecules in electron tomograms. International Journal of Molecular Sciences, 24(17):13375, 2023.

Claims

CLAIMS1. A method of template matching a tensorial template of a target particle to identify instances of that target within a tomogram obtained using electron tomography, such as cellular electron tomography, the method comprising: obtaining a tensorial template of the target particle, wherein the tensorial template describes the shape at all rotations of the target particle with a tensor field; and using the tensorial template to identify instances of that target particle within a tomogram.

2. The method of claim 1, wherein obtaining a tensorial template of the target particle comprises retrieving the tensorial template from a library of saved tensorial templates.

3. The method of claim 1, wherein obtaining a tensorial template of the target particle comprises generating the tensorial template, wherein the tensorial template comprises a tensor field of multiple tensors that describes the shape over all rotations of the target particle with a tensor field.

4. The method of claim 3, wherein: the target particle is described by a model comprising an array of voxels, and the tensorial template comprises a tensor per voxel of the model; templates are generated for different rotations of the target particle, each template comprising the array of voxels; and a tensor power is generated for each rotation of the target particle and wherein the tensor power comprises components that describe the shape of the target particle at that rotation.

5. The method of claim 4, wherein each tensor power is calculated as a tensor power of a quaternion that describes the shape of the target particle at that rotation, optionally a fourth power multiplication of the four components of a 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 that represents one 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 a value of the one voxel of the template of the target particle for the one rotation.

7. The method of any of claims 1 to 6, comprising masking the model using an operator that bounds the model of the target particle prior to generating the components of the tensors thereby setting the size of the array of voxels and, optionally, wherein the operator is a symmetric positive semidefinite linear operator.

8. The method of any of claims 1 to 7, wherein using the tensorial template to identify instances of the target particle within the tomogram comprises: convolving the tensorial template with the tomogram to obtain a convolved tensor field that provides information on the correlation of the tensorial template with the tomogram at all positions and at all rotations within the tomogram.

9. The method of claim 8, wherein using the tensorial template to identify instances of the target particle within the tomogram comprises: determining the best correlation value for each of the positions within the tomogram by calculating a scalar map of the convolved tensor field at each position.

10. The method of claim 9, wherein calculating a scalar map of the convolved tensor field at each position comprises calculating a scalar map of the total energy of the convolved tensor field at each position using a norm, optionally the Frobenius norm, and the method further comprises identifying positions of instances of the target particle within the tomogram by determining local maxima of the best correlation values.

11. The method of claim 10, further comprising: identifying rotations of the target particle at each of the identified positions using a global optimisation technique, optionally by identifying eigenvectors and eigenvalues of the convolved tensor field at each of the identified positions or by sampling the space of rotations of the target expressed as quaternions, and calculating the inner product between the selected tensors of the convolved tensor field and the quaternions.

12. The method of claim 11, wherein using the tensorial template to identify instances of the target particle within the tomogram comprises:calculating normalised correlation values for each identified instance of the target particle within the tomogram by calculating the local normalised cross-correlation of a non-tensorial template of the target particle with the tomogram using, as inputs, the position and rotation identified for each instance of the target particle within the tomogram and, optionally, the method further comprises generating each non-tensorial template from a model of the target particle to be a template describing the shape of the target particle at the corresponding identified rotation.

13. The method of any of claims 1 to 12, further comprising: dividing the tomogram into blocks; and for each block in turn: convolving the tensorial template with that block of the tomogram to obtain a convolved tensor field that provides information on the correlation of the tensorial template with the tomogram at all positions and at all rotations within that block of the tomogram; determining the best correlation value for each of the positions within that block of the tomogram thereby calculating a scalar map from the convolved tensor field; identifying positions of instances of the target particle within that block of the tomogram by determining local maxima of the best correlation values; identifying rotations of the target particle at each of the identified positions using a global optimisation technique; and, optionally, merging the rotations and positions identified in each block of the tomogram prior to calculating the normalised correlation values calculated for each identified instance of the target particle within the tomogram.

14. The method of claim 13, wherein dividing the tomogram into blocks comprises dividing the tomogram into contiguous units and defining each block to comprise one unit plus an overlap portion extending into the adjacent units and, optionally, the size of the overlap portion extension into adjacent units is less than the greatest dimension of the target particle and more than half the test dimension 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 of any preceding claim, acomputer readable medium having such a computer program stored thereon or as computer system comprising a computer processor and computer memory having such a computer program stored therein.

Images