System and method for chamber reconstruction from local volume

Abstract

A plurality of coordinates are mapped to respective first voxels. A value f(0) is assigned to each first voxel, f(n) being a monotonic function defined between 0 and M, inclusive. Respective second values {f(d(v i ))} are assigned to a plurality of second voxels {v i}, each of the plurality of second voxels {v i} being at a distance of 1 ≤ d(v i ) ≤ M voxels from a nearest first voxel. Each voxel of at least some of the second voxels is then iteratively assigned a weighted average of respective values of its immediate neighborhood, any value differing from f(M) by no more than a first threshold being given a higher weight than any other value. A subset of the second voxels is identified, each of the subset having a value differing from f(M) by more than a second threshold. Subsequently, a mesh representing a surface of a volume comprising the first voxels and the subset of the second voxels is generated.

Description

Technical Field

[0001] The present invention relates to the field of computer modeling. Background Art

[0002] Some medical applications require constructing mesh models of anatomical structures such as heart chambers.

[0003] Lorensen, William, and Cline, Harvey E., “Marchingcubes: a high resolution 3D surface construction algorithm” (ACM SIGGRAPHcomputer graphics, Vol. 21, No. 4, ACM, 1987), incorporated herein by reference, provides an algorithm called “marching cubes” that generates triangular models with constant density surfaces from 3D medical data. Summary of the Invention

[0004] According to some embodiments of the present invention, a system comprising a monitor and a processor is provided. The processor is configured to map corresponding coordinates of a plurality of locations within an anatomical cavity to corresponding first voxels, and assign a first value f(0) to each of the first voxels, where f(n) is a monotonic function defined for each integer n between 0 and M, inclusive, where M is a predefined positive integer. The processor is further configured to map the corresponding second value {f(d(v i ))} assigned to multiple second voxels {v i} to evaluate the plurality of second voxels {v i}, the plurality of second voxels {v i Each of the voxels is 1≤d(v i )≤M voxels. The processor is further configured to, after evaluating the second voxel, iteratively re-evaluate the second voxel by iteratively assigning to each of at least some of the second voxels a weighted average of corresponding values ​​of the voxel's immediate neighborhood, wherein any one of the values ​​that differs from f(M) by no more than a first predefined threshold Δ is given a higher weight than any other one of the values. The processor is further configured to, after iteratively re-evaluating the second voxel, identify a subset of the second voxels, each of the subset having a value that differs from f(M) by more than a second predefined threshold. The processor is further configured to generate a mesh representing a surface of a volume including the subset of the second voxels and the first voxels by applying a mesh generation algorithm to the volume, and to display the mesh on a monitor.

[0005] In some embodiments, the volume is composed of the first voxels and a subset of the second voxels.

[0006] In some embodiments, the mesh generation algorithm is a marching cubes algorithm.

[0007] In some embodiments, f(n) is monotonically decreasing over the domain [0, M].

[0008] In some embodiments, f(M)=0 and f(M-1)=Δ.

[0009] In some embodiments, for any integer n0 in the domain [0, M-3], f(n0+1) is a predefined percentage of f(n0).

[0010] In some embodiments, the processor is configured to calculate the value of d(v i )=j to evaluate the second voxel.

[0011] In some embodiments, the processor is configured to evaluate the second voxel using multiple parallel execution threads.

[0012] In some embodiments, the processor is configured to re-evaluate the second voxel using multiple parallel execution threads.

[0013] According to some embodiments of the present invention, a method is also provided, comprising mapping respective coordinates of a plurality of locations within an anatomical cavity to respective first voxels, and assigning a first value f(0) to each of the first voxels, f(n) being a monotonic function defined for each integer n between 0 and M, inclusive, where M is a predefined positive integer. The method further comprises assigning a first value {f(d(v i ))} assigned to multiple second voxels {v i} to evaluate the plurality of second voxels {v i}, the plurality of second voxels {v i Each of the voxels is 1≤d(v i )≤M voxels. The method further comprises, after evaluating the second voxel, iteratively re-evaluating the second voxel by iteratively assigning to each of at least some of the second voxels a weighted average of corresponding values ​​of the voxel's immediate neighborhood, wherein any one of the values ​​that differs from f(M) by no more than a first predefined threshold Δ is given a higher weight than any other one of the values. The method further comprises, after iteratively re-evaluating the second voxel, identifying a subset of the second voxels, each of the subset having a value that differs from f(M) by more than a second predefined threshold. The method further comprises generating a mesh representing a surface of the volume by applying a mesh generation algorithm to the volume including the subset of the second voxels and the first voxels.

[0014] In some embodiments, the anatomical cavity comprises a chamber of the heart.

[0015] According to some embodiments of the present invention, there is also provided a computer software product comprising a tangible, non-transitory computer-readable medium storing program instructions. The instructions, when read by a processor, cause the processor to map corresponding coordinates of a plurality of locations within an anatomical cavity to corresponding first voxels, and assign a first value f(0) to each of the first voxels, where f(n) is a monotonic function defined for each integer n between 0 and M, inclusive, where M is a predefined positive integer. The instructions also cause the processor to assign a corresponding second value {f(d(v i ))} assigned to multiple second voxels {v i} to evaluate the plurality of second voxels {v i}, the plurality of second voxels {v i Each of the voxels in the first voxel is 1≤d(v i )≤M voxels. The instructions further cause the processor, after evaluating the second voxel, to iteratively re-evaluate the second voxel by iteratively assigning to each voxel of at least some of the second voxels a weighted average of corresponding values ​​of the voxel's immediate neighborhood, wherein any one of the values ​​that differs from f(M) by no more than a first predefined threshold Δ is weighted higher than any other one of the values. The instructions further cause the processor, after iteratively re-evaluating the second voxel, to identify a subset of the second voxels, each having a value that differs from f(M) by more than a second predefined threshold. The instructions further cause the processor to generate a mesh representing a surface of the volume by applying a mesh generation algorithm to the volume including the subset of the second voxels and the first voxels.

[0016] The present disclosure will be more fully understood through the following detailed description of embodiments of the present invention in conjunction with the accompanying drawings, in which: BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 is a schematic diagram of a system for generating a mesh model of one or more chambers of a subject's heart according to some embodiments of the present invention;

[0018] Figure 2 is a flowchart of a method for generating a mesh model according to some embodiments of the present invention;

[0019] Figure 3 According to some embodiments of the present invention Figure 2 Schematic diagrams of various aspects of the method;

[0020] Figure 4is a schematic diagram of a technique for iteratively evaluating voxels according to some embodiments of the present invention; and

[0021] Figure 5A -B is a schematic diagram of a technique for iteratively re-evaluating voxels according to some embodiments of the present invention. DETAILED DESCRIPTION

[0022] Overview

[0023] Embodiments of the present invention provide a technique for constructing a mesh model (also referred to herein simply as a "mesh") of a cardiac chamber. According to this technique, as a catheter moves within the chamber (and optionally within an adjacent blood vessel), its position is continuously recorded by a tracking system. Each of these recorded positions, or "points," is added to a "point cloud" representing the chamber. Subsequently, a volume of voxels is constructed from the point cloud, and then a mesh is constructed from the volume.

[0024] A particular challenge addressed by embodiments of the present invention is that gaps in a voxel volume can complicate the construction of a mesh. For example, the marching cubes algorithm cited above in the background section can treat gaps as being outside the volume and thus mesh the interior of the volume. While other algorithms can be used to fill gaps when constructing a mesh, these algorithms can be slow and do not always efficiently handle all of the gaps.

[0025] To address this challenge, embodiments of the present invention provide a processor configured to execute a filling algorithm to fill most or all gaps in a volume. Through this filling, a fast and efficient mesh construction algorithm such as a marching cube algorithm can be used to construct a mesh.

[0026] The filling algorithm operates on a neighborhood of a volume consisting of (i) the original set of voxels originating from the point cloud (designated in this paper as {v c}) and (ii) other voxel sets {v o}, each of the other voxels in the set is related to the voxel belonging to {v c} is M or fewer voxels away. M is chosen to be large enough so that each voxel representing the subvolume of the chamber is closest to the nearest voxel belonging to {v c} is less than M voxels away.

[0027] When executing the fill algorithm, the processor first assigns the value 1 to {v c} and assigns a non-zero value to each voxel in {v o Specifically, the processor converts the value f(d(v i )) assigned to {v o Each voxel v in i, where d(v i ) is v i Distance {v c}, and f(n) is a monotonically decreasing function of n over the domain [0, M] such that f(0) = 1 and f(M) = 0. Advantageously, this assignment of values ​​can be performed on multiple parallel execution threads (e.g., using a graphics processing unit (GPU)) because the values ​​belonging to {v o} voxels can be evaluated in parallel with each other.

[0028] Next, the processor iteratively re-evaluates {v o} that are at a distance of less than M. Specifically, at each iteration, each of these voxels is assigned a weighted average of the values ​​of its immediate neighbors. In computing this weighted average, all non-zero values ​​receive equal weighting, while each zero value receives a significantly higher weighting than a non-zero value. Thus, if the number of iterations is large enough, the processor will zero all, or at least most, of the subvolumes that are outside the chamber. (As mentioned above, M is typically large enough so that no voxels within the chamber have a zero value, and thus voxels within the chamber are not zeroed). Advantageously, this re-evaluation can also be performed on multiple parallel execution threads.

[0029] The processor then {v o} are merged into the volume with those voxels in the volume that have a value greater than a predetermined threshold between 0 and 1. The processor thus fills the gaps in the volume.

[0030] Finally, the processor performs a fast mesh generation algorithm, such as a marching cubes algorithm, on the filled volume. The mesh can then be displayed on a computer screen. Optionally, the mesh can be superimposed with electrophysiological information, such as local activation time (LAT) values, to generate an electroanatomical map.

[0031] Although this description primarily relates to the chambers of the heart, it is noted that the techniques described herein can also be used to model any other anatomical or non-anatomical structure. For example, the techniques described herein can be used in depth sensing applications.

[0032] System Description

[0033] See first Figure 1 , which is a schematic diagram of a system 20 for generating a mesh model representing one or more chambers of a heart 24 of a subject 22, according to some embodiments of the present invention.

[0034] Figure 1A physician 30 is shown operating a catheter 26, the distal end 28 of which is positioned within the heart 24. As the physician 30 moves the distal end 28 of the catheter 26 within the chamber of the heart, the position of the distal end of the catheter is ascertained by a processor 32 belonging to the system 20.

[0035] For example, the distal end of the catheter may include one or more electromagnetic sensors that, in the presence of a generated magnetic field, output signals indicating the respective positions of the sensors. These signals may be received by the processor 32 via an electrical interface 34 (such as a port or socket). Based on these signals, the processor 32 may determine the position of the distal end of the catheter.

[0036] Alternatively or in addition, the distal end of the catheter may include a catheter electrode, and a plurality of electrode patches may be coupled to the body of the subject 22. When a voltage is applied between the catheter electrode and the electrode patch, the corresponding magnitude of the current between the catheter electrode and the electrode patch may be measured. Based on these current magnitudes, the processor may ascertain the position of the distal end of the catheter. Alternatively or in addition, any other suitable technique may be used to track the distal end of the catheter.

[0037] In some embodiments, the system 20 further includes a monitor 36. As the physician operates the catheter 26, the processor 32 can superimpose an icon representing the distal end of the catheter on the image of the subject's heart on the monitor 36, allowing the physician to visually track the position of the distal end. Alternatively or in addition, after the mesh model (as described in more detail below) is generated, the processor can display the mesh model on the monitor 36.

[0038] Generally speaking, the processor 32 can be implemented as a single processor or a group of collaborative networked or clustered processors. In some embodiments, as described herein, the functions of the processor 32 can be implemented in hardware using, for example, only one or more application-specific integrated circuits (ASICs) or field programmable gate arrays (FPGAs). In other embodiments, the functions of the processor 32 are implemented at least in part in software. For example, in some embodiments, the processor 32 is a programmed digital computing device comprising a central processing unit (CPU), a random access memory (RAM), a non-volatile auxiliary storage device (such as a hard drive or CD ROM drive), a network interface, and / or peripherals. Optionally, the processor 32 may also include a GPU. As is known in the art, program code and / or data comprising a software program are loaded into RAM for execution and processing by the CPU and / or GPU, and results are generated for display, output, transmission, or storage. For example, program code and / or data can be downloaded to a computer via a network in electronic form, or alternatively or additionally, they can be provided and / or stored on a non-temporary tangible medium (such as a magnetic memory, an optical memory, or an electronic memory). Such program code and / or data, when provided to a processor, produces a machine, or special purpose computer, configured to perform the tasks described herein.

[0039] although Figure 1 Specific applications are shown, but it is noted that the mesh generation techniques described herein can be used to model the surface of any structure. The structure can include any anatomical or non-anatomical cavity.

[0040] Mesh Generation

[0041] Now see Figure 2 , which is a flow chart of a method 38 for generating a mesh model executed by the processor 32 according to some embodiments of the present invention. Figure 3 , which is a schematic diagram of various aspects of method 38 according to some embodiments of the present invention. (For ease of illustration, this figure shows voxels as two-dimensional squares.)

[0042] The method 38 begins at a receiving step 40 where the processor receives respective coordinates 60 of a plurality of locations within a chamber of the heart 24. For example, as described above with reference to Figure 1 As described above, the coordinates 60 may be received from a position detection routine that detects the position of the distal end of the catheter 26 as the distal end moves within the lumen. Each coordinate 60 may also be referred to as a "point," and Figure 3 The set of coordinates 60 shown in portion A of FIG may be referred to as a “point cloud” 58. Point cloud 58 typically includes hundreds, thousands, or tens of thousands of points, as well as gaps 64 where no points exist.

[0043] Next, at a mapping step 42, the processor maps the coordinates 60 to corresponding first voxels 62, such that each first voxel 62 represents a different corresponding subvolume of the chamber in which at least one of the coordinates 60 exists. The processor also defines a plurality of second voxels 68, each of which represents a different corresponding subvolume of the gap 64 or the exterior of the chamber. The first voxels 62 and the second voxels 68 together define a continuous "voxel cloud" 59, as shown in FIG. Figure 3 As shown in Part B. (Although Figure 3 A specific example is shown, but it is noted that the voxel cloud 59 may have any suitable shape, such as a three-dimensional grid shape, for example, a cube shape. The processor also assigns a first value f(0) to each of the first voxels 62 at an assigning step 46, where f(n) is a function having the following properties.

[0044] For example, the processor may divide the volume including the point cloud 58 into a plurality of sub-volumes, each of which may have any suitable size, such as 0.8 mm × 0.8 mm × 0.8 mm. Each of the sub-volumes may then be represented by a corresponding voxel. Each voxel representing a sub-volume in which at least one of the coordinates 60 exists may be designated as a first voxel 62 and may be assigned a value of f(0).

[0045] exist Figure 3 In section B of FIG5 , first voxels 62 (each located within the boundary 66 of the point cloud 58) are represented by a first pattern, thereby indicating that these voxels have been assigned a value. In contrast, second voxels 68, which at this stage only have a default "placeholder" value such as negative infinity ("-inf"), are represented by a second pattern. As shown in this section of the figure, some second voxels are located within the boundary 66, while other second voxels are located between the boundary 66 and the perimeter 72 of the voxel cloud.

[0046] In general, f(n) is a monotonic function defined for every integer n between 0 and M, inclusive, where M is a predefined positive integer. (In general, M can vary significantly between applications.) For example, f(n) can be monotonically decreasing over the domain [0, M]. Typically, in such embodiments, f(0) is 1, and f(M) is 0.

[0047] In some embodiments, for any integer n0 in the domain [0, M-3], f(n0+1) is a predefined percentage of f(n0). For example, for embodiments in which f(n) is monotonically decreasing over the domain [0, M], f(n0+1) may be equal to p*f(n0), where p is between 0.8 and 0.95 over the domain [0, M-3]. In some such embodiments, f(M-1) is equal to p*f(M-2). Thus, for example, for M=10, f(M)=0, and p=0.9, the values ​​that f(n) may take on over the domain [0, M] are 1, 0.9, 0.81, 0.73, 0.66, 0.59, 0.53, 0.48, 0.43, 0.39, and 0. In other such embodiments, f(M-1) is equal to a decimal Δ, such as any number less than or equal to 0.0001, and f(M)=0.

[0048] After the mapping step 42 and the assigning step 46, the processor evaluates the second voxels 68 (i.e., assigns corresponding values ​​to the second voxels) at an evaluation step 48. Specifically, the processor assigns a second value f(d) to each second voxel 68, where d is the distance (in voxels) of the second voxel to the nearest first voxel.

[0049] In some embodiments, the second voxel 68 is evaluated by evaluating the second voxel for d=j in every j-th iteration throughout the M iterations. In this regard, see now Figure 4 , which is a schematic diagram of a technique for iteratively evaluating a second voxel 68 , according to some embodiments of the present invention.

[0050] like Figure 4 As shown, during the first of the M iterations, a first subset of second voxels 68 is assigned a first value f(1), each of which is a direct neighbor of (i.e., adjacent to) the first voxel 62. Similarly, during the second iteration, a second subset of second voxels is assigned a value f(2), each of which is a direct neighbor of a member of the first subset. Each subsequent iteration then evaluates another subset of the second voxels. (In implementation, typically, during each jth iteration of the total M iterations, all of the second voxels that have not yet been evaluated are processed. Those of the processed voxels that have direct neighbors with a value of f(j-1) are assigned a value of f(j), while the other voxels are not evaluated.)

[0051] In some embodiments, a voxel is called a direct neighbor (or "adjacent") of another voxel if the two voxels share at least one vertex. Thus, a voxel can have up to 26 direct neighbors. (This criterion is presented in Figure 4 The difference is that due to the two-dimensional representation of voxels, Figure 4( 8 rather than 26 direct neighbors are shown.) In other embodiments, two voxels are called direct neighbors of each other only if they share at least one face; thus, a voxel may only have a maximum of six direct neighbors. In other embodiments, other criteria may also be used to determine the direct neighbors of a voxel.

[0052] In some embodiments, multiple parallel execution threads (e.g., running on a GPU) are used to evaluate the second voxels 68. For example, during each of the M iterations, all of the second voxels that have not yet been evaluated may be processed in parallel.

[0053] The effect of evaluating step 48 is shown in Figure 3 , wherein voxels near perimeter 72 are represented by white space, thereby indicating that these voxels have been assigned the value f(M) (which in many embodiments is zero, as described above), while the remainder of the second voxels (including the voxels filling gap 64) are represented by the first pattern, thereby indicating that these voxels have been assigned other values. (Although portion C does not show any unevaluated second voxels, it is noted that in implementations, at least some voxels in voxel cloud 59 may be located at a distance greater than M from the nearest first voxel, such that these voxels may remain unevaluated (i.e., may remain with a default initial value). However, these voxels are not used in method 38.)

[0054] After the evaluation step 48, the processor performs a predefined number of iterations (e.g., between 100 and 200 iterations) across the second voxel at a re-evaluation step 50. During each of these iterations, the processor re-evaluates the second voxel by assigning a weighted average of the corresponding values ​​of the voxel's immediate neighborhood to each voxel of at least some of the second voxels (typically, each of the second voxels that does not have a value f(M)). Typically, the weighted average assigns a higher weight to any value that differs from f(M) by no more than a predefined threshold Δ. (As described above, in some embodiments, f(M-1) = Δ and f(M) = 0.) For example, the weighted average may assign an equal weight of 1 to all values ​​that differ from f(M) by more than Δ, while assigning a weight of, for example, 40 or more to any value that differs from f(M) by Δ or less.

[0055] In this regard, see also Figure 5A -B, which is a schematic diagram of a technique according to some embodiments of the present invention for iteratively re-evaluating second voxels 68. By way of example, each of these figures shows two re-evaluation iterations for a small representative set of second voxels 68, where any voxels that may be outside the set are ignored. Figure 5A -B assumes f(M) = 0 and gives equal weight to all voxel values ​​greater than Δ. For convenience, Figure 5BIn the example, values ​​not greater than Δ are written as 0.

[0056] Figure 5A -B have no difference in the initial voxel values, with only one difference: Figure 5B However, due to this difference, the effect of the revaluation is Figure 5B Specifically, Figure 5A In , the voxel values ​​are "smooth" over the entire set of voxels, but Figure 5B , the higher weight assigned to voxel 68a causes the entire set to become zero valued.

[0057] In general, if M is large enough and the gap 64 is not abnormally large, the portion of the voxel cloud 59 within the boundary 66 (representing the interior of the chamber) is smooth, as shown in FIG. Figure 5A As shown, this portion generally does not include voxels with a value of f(M). On the other hand, the portion of the voxel cloud 59 outside the boundary 66 (representing the interior of the chamber) is generally dominated by f(M) due to the presence of voxels with a value of f(M) in this portion, as shown in FIG. Figure 5B shown.

[0058] In some embodiments, multiple parallel execution threads are used to re-evaluate the second voxel, as described above for evaluation step 48 .

[0059] After the re-evaluation step 50, the processor identifies a subset of second voxels at a subset identification step 52, each of which has a value that differs from f(M) by no more than another predefined threshold α. For example, for an embodiment where f(M)=0, the processor may identify those second voxels whose value is greater than α, where α is, for example, between 0.5 and 0.8. Subsequently, at a volume definition step 54, and as Figure 3 As shown in section D of , the processor defines a volume 74 that includes the first voxel 62 and the identified subset of the second voxel (typically, consisting of these voxels). Generally speaking, the volume 74 includes most or all of the voxels within the boundary 66, while excluding most or all of the voxels outside the boundary 66.

[0060] Typically, the processor defines a table of values ​​for the voxel cloud 59. The processor updates the values ​​in the table as the assigning step 46, the evaluating step 48, and the re-evaluating step 50 are performed. Subsequently, at the volume defining step 54, the processor defines the volume 74 by setting those second voxels that differ from f(M) by more than α to f(0) and setting the remainder of the second voxels to f(M). Thus, for example, those voxels that belong to the volume 74 may be listed in the table as a value of 1, while those voxels that do not belong to the volume 74 may be listed as a value of 0.

[0061] Subsequently, at a mesh generation step 56 , the processor generates a mesh 76 , such as a triangular mesh, representing the surface of the volume 74 by applying a mesh generation algorithm to the volume 74 . Figure 3 Part E of FIG. 7 shows a mesh 76 , which indicates the three-dimensional nature of the mesh by shading some of the triangles in the mesh.

[0062] In some embodiments, the processor uses a marching cube algorithm to generate the mesh 76. This algorithm typically requires as input a function defining the volume whose surface is to be meshed. To meet this requirement, the processor can provide the aforementioned value table indicating those voxels that belong to the volume 74.

[0063] It will be understood by those skilled in the art that the present invention is not limited to the contents specifically shown and described above. On the contrary, the scope of the embodiments of the present invention includes both the combination and sub-combination of the various features described above, as well as variations and modifications that are not within the scope of the prior art that may occur to those skilled in the art when reading the above description. The documents incorporated by reference into this patent application are considered to be an integral part of this application, except that if any term defined in these incorporated documents conflicts with the definition explicitly or implicitly given in this specification, only the definition in this specification should be considered.

Claims

1. A system comprising: Monitor; and a processor configured to: mapping respective coordinates of a plurality of locations within the anatomical cavity to respective first voxels, assigning a first value f(0) to each of said first voxels, f(n) is a monotonic function defined for every integer n between 0 and M, inclusive, where M is a predefined positive integer, By replacing the corresponding second value {f(d(v i ))} assigned to the second voxel to evaluate multiple second voxels {v i }, the plurality of second voxels {v i } is 1≤d(v i )≤M voxels, After evaluating the second voxel, iteratively re-evaluating the second voxel by iteratively assigning to each voxel of at least some of the second voxels a weighted average of corresponding values ​​of the voxel's immediate neighborhood, In said weighted average, any one of said values ​​that differs from f(M) by no more than a first predefined threshold Δ is given a higher weight than any other one of said values, After iteratively re-evaluating the second voxels, identifying a subset of the second voxels, each of the subset having a value that differs from f(M) by more than a second predefined threshold, generating a mesh representing a surface of the volume by applying a mesh generation algorithm to the volume including the first voxels and the subset of the second voxels, and The grid is displayed on the monitor. 2 . The system of claim 1 , wherein the volume is comprised of the subset of the first and second voxels. The system of claim 1 , wherein the mesh generation algorithm is a marching cubes algorithm. The system of claim 1 , wherein f(n) is monotonically decreasing over the domain [0, M]. The system of claim 4 , wherein f(M)=0 and f(M−1)=Δ.

6. The system of claim 1, wherein for any integer n0 in the domain [0, M-3], f(n0+1) is a predefined percentage of f(n0).

7. The system of claim 1 , wherein the processor is configured to evaluate the value of d(v i )=j to evaluate the second voxel.

8. The system of claim 7, wherein the processor is configured to evaluate the second voxel using multiple parallel execution threads.

9. The system of claim 1, wherein the processor is configured to re-evaluate the second voxel using multiple parallel execution threads.

10. A method comprising: mapping respective coordinates of a plurality of locations within the anatomical cavity to respective first voxels; assigning a first value f(0) to each of said first voxels, f(n) is a monotonic function defined for every integer n between 0 and M, inclusive, where M is a predefined positive integer; By replacing the corresponding second value {f(d(v i ))} assigned to the second voxel to evaluate multiple second voxels {v i }, the plurality of second voxels {v i } is 1≤d(v i )≤M voxel distance; After evaluating the second voxel, iteratively re-evaluating the second voxel by iteratively assigning to each voxel of at least some of the second voxels a weighted average of corresponding values ​​of the voxel's immediate neighborhood, In said weighted average, any one of said values ​​that differs from f(M) by no more than a first predefined threshold Δ is given a higher weight than any other one of said values; After iteratively re-evaluating the second voxels, identifying a subset of the second voxels, each of the subset having a value that differs from f(M) by more than a second predefined threshold; as well as A mesh representing a surface of the volume is generated by applying a mesh generation algorithm to the volume including the first voxels and the subset of the second voxels. The method of claim 10 , wherein the volume is comprised of the subset of the first and second voxels.

12. The method of claim 10, wherein the anatomical cavity comprises a chamber of the heart. The method of claim 10 , wherein the mesh generation algorithm is a marching cubes algorithm. The method of claim 10 , wherein f(n) is monotonically decreasing over the domain [0, M]. The method according to claim 14 , wherein f(M)=0 and f(M−1)=Δ.

16. The method of claim 10, wherein for any integer n0 in the domain [0, M-3], f(n0+1) is a predefined percentage of f(n0). 17 . The method of claim 10 , wherein evaluating the second voxels comprises evaluating the second voxels by evaluating those of the second voxels for which d(vi)=j in every jth of the iterations throughout M iterations. The method of claim 17 , wherein evaluating the second voxel comprises evaluating the second voxel using a plurality of parallel execution threads. The method of claim 10 , wherein re-evaluating the second voxel comprises re-evaluating the second voxel using a plurality of parallel execution threads.

20. A computer software product comprising a tangible, non-transitory computer-readable medium having program instructions stored therein, the instructions, when read by a processor, causing the processor to: mapping respective coordinates of a plurality of locations within the anatomical cavity to respective first voxels, assigning a first value f(0) to each of said first voxels, f(n) is a monotonic function defined for every integer n between 0 and M, inclusive, where M is a predefined positive integer, By replacing the corresponding second value {f(d(v i ))} assigned to the second voxel to evaluate multiple second voxels {v i }, the plurality of second voxels {v i } is 1≤d(v i )≤M voxels, After evaluating the second voxel, iteratively re-evaluating the second voxel by iteratively assigning to each voxel of at least some of the second voxels a weighted average of corresponding values ​​of the voxel's immediate neighborhood, In said weighted average, any one of said values ​​that differs from f(M) by no more than a first predefined threshold Δ is given a higher weight than any other one of said values, After iteratively re-evaluating the second voxels, identifying a subset of the second voxels, each of the subset having a value that differs from f(M) by more than a second predefined threshold, and A mesh representing a surface of the volume is generated by applying a mesh generation algorithm to the volume including the first voxels and the subset of the second voxels.

Images