Abstract:
The present invention relates to a method for characterizing a blood vessel represented by vascular image data, wherein said vascular image data comprises a plurality of voxels each having an image intensity, said method comprising the steps of identifying a set of voxels representing a boundary of the blood vessel; determining a gradient vector of the image intensity for each voxel in said set of voxels representing the boundary of the blood vessel; selecting, from said set of voxels representing the boundary of the blood vessel, a subset of voxels such that all voxels have a common intersection point for their respective gradient vector extensions; and determining a vector product based on said gradient vectors for said subset of voxels, wherein the common intersection point indicates a centre of said blood vessel and said vector product indicates a direction in which said blood vessel extends.

Description:
TECHNICAL FIELD 
       [0001]    The present invention relates to a method for characterizing a blood vessel represented by vascular image data. 
       BACKGROUND OF THE INVENTION 
       [0002]    Angiography is a medical imaging technique to image blood vessels and blood vessel systems, and can be used for example to enable diagnostics, surgical planning, and simulation. 
         [0003]    Angiographic image data can be generated by injecting a contrast agent intravenously to a patient and scanning the patient using e.g. Computed Tomography (CT) or Magnetic Resonance Imaging (MRI). The resulting angiographic image data is typically a set of DICOM-images forming a voxel volume with discrete image intensities. The voxel volume can be processed to segment organs, blood vessels, skeleton etc. In particular, the segmentation and further characterization can provide a segmented surface that delimits the blood vessel or vessel system from other tissues. 
         [0004]    An example of such a characterization method for characterizing a vascular system in a three-dimensional angiographic image comprised of voxels is disclosed in U.S. Pat. No. 6,842,638. Here a two-dimensional slice formed of pixels is extracted from the angiographic image. Imaged vascular structures in the slice are located and flood-filled. The edges of the filled regions are iteratively eroded to identify vessel centres. The extracting, locating, flood-filling, and eroding is repeated for a plurality of slices to generate a plurality of vessel centres that are representative of the vascular system. A vessel centre is selected, and a corresponding vessel direction and orthogonal plane are found. The vessel boundaries in the orthogonal plane are identified by iteratively propagating a closed geometric contour arranged about the vessel centre. The selecting, finding, and estimating are repeated for the plurality of vessel centres, and the estimated vessel boundaries are then interpolated to form a vascular tree. 
         [0005]    However, sometimes a portion of the boundary of a blood vessel may be “missing”, or not discernable, such that the boundary does not continuously enclose the blood vessel in the dimensional slice. If so, the flood-filling procedure will fail as it expands outside the blood vessel. This makes the method suggested in U.S. Pat. No. 6,842,638 sensitive to imperfections in the DICOM-image. Thus, there is a need for a more robust method to characterize a blood vessel. 
       SUMMARY OF THE INVENTION 
       [0006]    In view of the above, an object of the invention is to at least alleviate the problem discussed above. According to an aspect of the invention, there is provided a method for characterizing a blood vessel represented by vascular image data from a DICOM-image, wherein the vascular image data comprises a plurality of voxels each having an image intensity. The method comprises the steps of identifying a set of voxels representing at least a portion of a boundary of the blood vessel, determining a gradient vector of the image intensity for each voxel in the set of voxels representing the boundary of the blood vessel, selecting, from the set of voxels representing the boundary of the blood vessel, a subset of voxels such that all voxels in the subset of voxels have a common intersection point for their respective gradient vector extensions, and determining a vector product based on the gradient vectors for the subset of voxels, wherein the common intersection point indicates a centre of the blood vessel and the vector product indicates a direction in which the blood vessel extends. 
         [0007]    The present invention is based on the understanding that since a blood vessel typically has a tubular shape, the gradient vectors of the image intensities at the boundary of the blood vessel will be directed radially outwards and thus the point where extensions of these gradient vectors meet may serve as an indication of the centre of the blood vessel. Furthermore, the vector product between two gradient vectors located at the boundary can be used to indicate the direction in which the blood vessel extends. The inventive method is advantageous in that even though a portion of the boundary is missing, or not discernable, the centre of the blood vessel may still be found from the available portion of the boundary. 
         [0008]    A boundary is here intended to indicate a segment of the blood vessel that delimits the blood vessel from other tissue. It shall be noted that the boundary may be continuous or discontinuous (i.e. a portion of the boundary may be missing). The common intersection point may typically be a small discrete volume. For example, the subset may include all voxels having gradient vector extensions that pass through a common voxel. Further, each gradient vector extension starts from the voxel in which the gradient vector is calculated. 
         [0009]    A set of voxels representing a boundary of the blood vessel can be found by computing a gradient vector of the image intensity for each of the plurality of voxels, and selecting a set of voxels having a gradient vector norm that exceeds a predetermined threshold value. Note that, in this case, the gradient vectors for the voxels are determined before the set of voxels representing the boundary of the blood vessel is identified, instead of after. 
         [0010]    According to an embodiment, any voxel that has an adjacent voxel, in a direction of the gradient vector, with a larger gradient vector norm may be excluded from the set of voxels representing the boundary of a blood vessel. An advantage is that this may result in a thinner, more well-defined boundary region. This can also be achieved by only including voxels where a second derivative of the image intensity is essentially zero. 
         [0011]    Preferably, the gradient vector extension may be adapted to extend in a direction opposite to the direction of the gradient vector. Furthermore, a length of the gradient vector extension may be determined by a predetermined end condition. An example of an end condition would be that the gradient vector extension only extends until the opposite boundary of the blood vessel is reached, or that the gradient vector extension extends a predetermined length. The predetermined length may preferably be set sufficiently low to prevent that the gradient vector extension intersects gradient vector extensions associated with other blood vessels. An advantage with putting restrictions on the length of the gradient vector extension is that the risk of unintentionally including gradient vector extensions from adjacent blood vessels in the subset of voxels having a common intersection point for their respective gradient vector extensions is reduced. 
         [0012]    The method may further comprise the step of iteratively selecting an additional subset of voxels from the set of voxels representing the boundary of the blood vessel, wherein all voxels in the additional subset of voxels have a common intersection point for their respective gradient vector extensions and the common intersection point indicates a possible centre point of the blood vessel. 
         [0013]    Furthermore, the method according to the present invention may advantageously be included in a method for determining a blood vessel path by: finding a plurality of possible centre points of the blood vessel, selecting a voxel that corresponds to one of the centre points as a starting point for the blood vessel, tracing the blood vessel path by iteratively: identifying a possible centre point that is located in a plane that is orthogonal to the direction of the blood vessel, wherein the plane extends through a voxel that is located within a predetermined range from (e.g. next to) the current centre point in a direction of the blood vessel, including the identified possible centre point in the blood vessel path. The predetermined range may preferably be one voxel thereby making full use of the resolution in the DICOM-image. 
         [0014]    Additionally, when there is more than one possible centre point located in the plane, the centre point to be included in the vessel path may be selected by iteratively comparing the gradient vector norm of the current voxel with the gradient vector norm of all neighbouring voxels located in the plane and changing to the neighbouring voxel located in the plane that has the largest gradient vector norm until a voxel is found having a gradient vector norm larger than all neighbouring voxels in the plane. 
         [0015]    Tracing of the blood vessel path may preferably stop when an end condition is fulfilled. The end condition may, for example, be at least one of that the boundary of the voxel volume is reached, that a blood vessel path that has already been traced is reached, that the intensity of the centre point is below a predetermined threshold value and that the vector product associated with the centre point is a zero vector. The method may also comprise the step of removing vessel paths that are not attached to the rest of the vessel structure. This removes objects that are not blood vessels. 
         [0016]    Also, the method according to the present invention may advantageously be included in a method for creating a blood vessel system by: determining a plurality of blood vessel paths as described above, removing blood vessel paths that are not attached to the rest of the blood vessel system. Additionally, after a blood vessel path has been determined, a blood vessel surface may be created by sampling in a plane normal to a direction of the path. 
         [0017]    According to another aspect of the present invention there is provided a device for characterizing a blood vessel represented by vascular image data from a DICOM-image, wherein said vascular image data comprises a plurality of voxels each having an image intensity, said device comprising means for identifying a set of voxels representing at least a portion of a boundary of the blood vessel, determining a gradient vector of the image intensity for each voxel in said set of voxels representing the boundary of the blood vessel, selecting, from said set of voxels representing the boundary of the blood vessel, a subset of voxels such that all voxels in said subset of voxels have a common intersection point for their respective gradient vector extensions, and determining a vector product based on said gradient vectors for said subset of voxels, wherein the common intersection point indicates a centre of said blood vessel and said vector product indicates a direction in which said blood vessel extends. This aspect of the invention provides similar advantages as discussed above in relation to the previous aspect of the invention. Additionally, the inventive device may be comprised in a blood vessel characterizing system, further comprising means for acquiring a DICOM-image. 
         [0018]    According to a still further aspect of the present invention there is provided a computer program product comprising a computer readable medium having stored thereon computer program means for causing a computer to provide a method for characterizing a blood vessel represented by vascular image data from a DICOM-image, wherein said vascular image data comprises a plurality of voxels each having an image intensity, visual categorization method, wherein the computer program product comprises code for identifying a set of voxels representing at least a portion of a boundary of the blood vessel, code for determining a gradient vector of the image intensity for each voxel in said set of voxels representing the boundary of the blood vessel, code for selecting, from said set of voxels representing the boundary of the blood vessel, a subset of voxels such that all voxels in said subset of voxels have a common intersection point for their respective gradient vector extensions, and code for determining a vector product based on said gradient vectors for said subset of voxels, wherein the common intersection point indicates a centre of said blood vessel and said vector product indicates a direction in which said blood vessel extends. This aspect of the invention provides similar advantages as discussed above in relation to the previous aspects of the invention. Also, the computer readable medium may be one of a removable nonvolatile random access memory, a hard disk drive, a floppy disk, a CD-ROM, a DVD-ROM, a USB memory, an SD memory card, or a similar computer readable medium known in the art. 
         [0019]    Further features of, and advantages with, the present invention will become apparent when studying the appended claims and the following description. The skilled addressee realize that different features of the present invention may be combined to create embodiments other than those described in the following, without departing from the scope of the present invention. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0020]    The above, as well as additional objects, features and advantages of the present invention, will be better understood through the following illustrative and non-limiting detailed description of preferred embodiments of the present invention, with reference to the appended drawings, where the same reference numerals will be used for similar elements, wherein: 
           [0021]      FIG. 1   a - 1   c  shows schematic flow charts of a currently preferred method for characterizing a blood vessel; 
           [0022]      FIG. 2  schematically illustrates volume imaging of a patient; 
           [0023]      FIG. 3  illustrates the extension of an exemplary blood vessel; 
           [0024]      FIG. 4  illustrates an example of a plurality of paths; 
           [0025]      FIG. 5  illustrates an example of paths of a blood vessel system; and 
           [0026]      FIG. 6  illustrates an example of a blood vessel system after a vessel surface has been generated around the paths. 
       
    
    
     DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 
       [0027]    The present invention will now be described more fully hereinafter with reference to the accompanying drawings, in which currently preferred embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein; rather, these embodiments are provided for thoroughness and completeness, and fully convey the scope of the invention to those skilled in the art. Like reference characters refer to like elements throughout. 
         [0028]    Referring now to the drawings and  FIGS. 1   a - 1   c  in particular, there is depicted a schematic flow chart illustrating an embodiment of a method for characterizing a blood vessel represented by vascular image data from a DICOM-image. Reference will in parallel be made also to  FIGS. 2 and 3 . 
         [0029]    As indicated by step  101 , the vascular system of a patient may be imaged by scanning a patient  201  using e.g. Computed Tomography (CT) or Magnetic Resonance Imaging (MRI). The output is typically a set of DICOM-images  202 . Each DICOM-image is a two-dimensional gray scale image representation of the examined area of the patient, which has good contrast for the vascular system relative to other body tissues of the patient. An example of a DICOM-image is depicted in  FIG. 2 , showing the aorta  203 , celiac trunk and left gastric artery  204  of the patient. 
         [0030]    The image intensities of each DICOM-image may be adjusted, in step  102 , to remove areas of low image intensity that represents hollow portions and regions outside the patient&#39;s body. This may be achieved by setting all pixels with a radio density below the radio density for water to a have a value that corresponds to the radio density of water. The radio density of water may vary due to the scale used but is typically 1000 (Hounsfield Units), or 0. 
         [0031]    Thereafter, in step  103 , three dimensional image data is generated from the set of DICOM-images to provide volume imaging of the patient. The result is a three dimensional array V(x,y,z), where each voxel x,y,z has an associated image intensity. Thus, V(x,y,z) provides a three-dimensional gray scale image representation of the examined area of the patient. To provide a clear and concise notation, each voxel&#39;s absolute coordinate (x,y,z) is hereafter denoted as an array p. In step  104 , a low pass filter may preferably be supplied to the voxel volume V to suppress noise. 
         [0032]    Further, in step  105 , a gradient vector of the image intensity in V is determined. The result is an array G indicating the gradient vector for each voxel in the volume V. 
         [0000]    
       
         
           
             
               
                 
                   
                     G 
                      
                     
                       ( 
                       p 
                       ) 
                     
                   
                   = 
                     
                    
                   
                     ∇ 
                     
                       V 
                        
                       
                         ( 
                         p 
                         ) 
                       
                     
                   
                 
               
             
             
               
                 
                   = 
                     
                    
                   
                     ∇ 
                     
                       V 
                        
                       
                         ( 
                         
                           x 
                           , 
                           y 
                           , 
                           z 
                         
                         ) 
                       
                     
                   
                 
               
             
             
               
                 
                   = 
                     
                    
                   
                     ( 
                     
                       
                         
                           ∂ 
                           V 
                         
                         
                           ∂ 
                           x 
                         
                       
                       , 
                       
                         
                           ∂ 
                           V 
                         
                         
                           ∂ 
                           y 
                         
                       
                       , 
                       
                         
                           ∂ 
                           V 
                         
                         
                           ∂ 
                           z 
                         
                       
                     
                     ) 
                   
                 
               
             
             
               
                 
                   = 
                     
                    
                   
                     
                       
                         x 
                         ^ 
                       
                        
                       
                         
                           ∂ 
                           V 
                         
                         
                           ∂ 
                           x 
                         
                       
                     
                     + 
                     
                       
                         y 
                         ^ 
                       
                        
                       
                         
                           ∂ 
                           V 
                         
                         
                           ∂ 
                           y 
                         
                       
                     
                     + 
                     
                       
                         z 
                         ^ 
                       
                        
                       
                         
                           
                             ∂ 
                             V 
                           
                           
                             ∂ 
                             z 
                           
                         
                         . 
                       
                     
                   
                 
               
             
           
         
       
     
         [0033]    G can be calculated by convolving V for the x, y, and z-component using a direction operator, such as a Gaussian kernel. The Gaussian kernel may be expressed as: 
         [0000]    
       
         
           
             
               
                 H 
                 g 
               
               = 
               
                 C 
                 
                   - 
                   
                     
                       
                         x 
                         2 
                       
                       + 
                       
                         y 
                         2 
                       
                       + 
                       
                         z 
                         2 
                       
                     
                     
                       σ 
                       2 
                     
                   
                 
               
             
             , 
           
         
       
     
         [0000]    where σ represents the standard deviation. 
         [0034]    Further, the differentiated Gaussian kernel may be expressed as: 
         [0000]    
       
         
           
             
               ∇ 
               
                 H 
                 g 
               
             
             = 
             
               
                 ( 
                 
                   
                     H 
                     gx 
                   
                   , 
                   
                     H 
                     gy 
                   
                   , 
                   
                     H 
                     gz 
                   
                 
                 ) 
               
               = 
               
                 
                   ∇ 
                   
                     ( 
                     
                       C 
                       
                         - 
                         
                           
                             
                               x 
                               2 
                             
                             + 
                             
                               y 
                               2 
                             
                             + 
                             
                               z 
                               2 
                             
                           
                           
                             σ 
                             2 
                           
                         
                       
                     
                     ) 
                   
                 
                 . 
               
             
           
         
       
     
         [0035]    Thus the gradient vector for each voxel in the volume V can be found by: 
         [0000]        G =( V*H   gx   ,V*H   gy   ,V*H   gz ), where * represents a convolution. 
         [0036]    Then, in step  106 , a set of voxels representing the boundaries of the blood vessels are identified by finding voxels with a gradient vector norm |G| that exceeds a predetermined threshold value. The threshold value depends on the image resolution, and may typically be selected as a constant times the least spatial resolution of the image data. 
         [0037]    This yields a binary volume C where 
         [0000]    
       
         
           
             
               C 
                
               
                 ( 
                 q 
                 ) 
               
             
             = 
             
               
                 [ 
                 
                   
                     
                       0 
                     
                     
                       
                         
                           when 
                            
                           
                               
                           
                            
                           
                              
                             
                               G 
                                
                               
                                 ( 
                                 q 
                                 ) 
                               
                             
                              
                           
                         
                         ≤ 
                         Threshold 
                       
                     
                   
                   
                     
                       1 
                     
                     
                       
                         
                           when 
                            
                           
                               
                           
                            
                           
                              
                             
                               G 
                                
                               
                                 ( 
                                 q 
                                 ) 
                               
                             
                              
                           
                         
                         &gt; 
                         Threshold 
                       
                     
                   
                 
                 ] 
               
               . 
             
           
         
       
     
         [0000]    where q is an array representing a discrete coordinate (i,j,k) of an voxel in volume C. Thus, voxels in C with a value of 1 indicates that a corresponding voxel in V is part of a boundary of a blood vessel. 
         [0038]    The boundary may be refined, in step  107 , to only include voxels in V for which the second derivative of the image intensity is essentially zero. That is: 
         [0000]    
       
         
           
             
               C 
                
               
                 ( 
                 q 
                 ) 
               
             
             = 
             
               
                 [ 
                 
                   
                     
                       0 
                     
                     
                       
                         
                           when 
                            
                           
                               
                           
                            
                           
                              
                             
                               G 
                                
                               
                                 ( 
                                 q 
                                 ) 
                               
                             
                              
                           
                         
                         ≤ 
                         
                           Threshold 
                            
                           
                               
                           
                            
                           or 
                            
                           
                               
                           
                            
                           
                             
                               ∇ 
                               2 
                             
                              
                             V 
                           
                         
                         ≠ 
                         0 
                       
                     
                   
                   
                     
                       1 
                     
                     
                       
                         
                           when 
                            
                           
                               
                           
                            
                           
                              
                             
                               G 
                                
                               
                                 ( 
                                 q 
                                 ) 
                               
                             
                              
                           
                         
                         &gt; 
                         
                           Threshold 
                            
                           
                               
                           
                            
                           and 
                            
                           
                               
                           
                            
                           
                             
                               ∇ 
                               2 
                             
                              
                             V 
                           
                         
                         ≈ 
                         0 
                       
                     
                   
                 
                 ] 
               
               . 
             
           
         
       
     
         [0039]    This provides a more well-defined boundary. However, instead of using a procedure that calculates the second derivative, the refinement of the boundary may be implemented in an iterative procedure as follows. For each voxel in C with a value 1, it is determined whether the corresponding voxel in V has an adjacent voxel in V in a direction of the gradient vector with a larger gradient vector norm. If so, the voxel in C is set to 0. This provides a more robust solution. 
         [0040]    The centre and direction of the blood vessels may then be determined as illustrated by step  108  to  112  in  FIG. 1   b . In step  108 , a voxel q 0  with a value 1 is selected from the volume C, and in step  109 , a gradient vector extension is created, here visualized as a straight line  301 , starting at the selected voxel q 0  and extending in a direction opposite to the direction gradient vector G(q 0 ). The gradient vector extension extends until the image intensity falls outside a predetermined range. This range may preferably be selected such that the gradient vector extension stops at the opposite boundary of the blood vessel. The gradient vector extension may also stop at a predetermined maximum length, which may be an estimate of a maximum diameter of a blood vessel. 
         [0041]    The starting point q 0  of the gradient vector extension  301  is stored in a trace index buffer M(q)=(M i (q), M j (q), M k (q)) for each discrete coordinate that the vector extension visits. 
         [0042]    In step  110 , a new voxel q 1  with a value 1 is selected from the volume C, and in step  111 , a gradient vector extension  302  is created starting at the selected voxel q 1 . Further, the starting point q 1  of the gradient vector extension is stored in the trace index buffer M(q) for each discrete coordinate that the vector extension visits. If there already is another starting point (e.g. q 0 ) stored in the trace index buffer M(q) this indicates that the gradient vector extension from q 1  intersects a gradient vector extension associated with a the previously stored starting point (e.g. q 0 ) and an intersection vector, v isect , is calculated as a vector product between these gradient vector extensions in the starting points. 
         [0043]    In the example illustrated in  FIG. 3 , the gradient vector extension  302  from q 1  intersects the gradient vector extension  301  from q 0 , at the discrete coordinate q isect , and an intersection vector v isect  is calculated as v isect =G(q 0 )×G(q 1 ) in step  112 , where x represents the vector product. This intersection vector, v isect , will have a direction perpendicular to both G(q 0 ) and G(q 1 ). The intersection vector is now extended from the intersection coordinate q isect , both in positive and negative direction. The trace from the extended intersection vector is included in two buffers; one trace visit buffer F c (q) and one trace direction buffer F v (q)=(F x (q), F y (q), F z (q)). The trace visit buffer will be increased with 1 at each coordinate q that the extended intersection vector visits. The trace direction buffer will be averaged with its previously stored direction (if any) and the intersection vector at each coordinate q that the extended intersection vector visits. 
         [0044]    Steps  110  to  112  are repeated for all voxels in C with a value 1. 
         [0045]    The blood vessel path may then be determined according to the procedure described below with reference to step  113  to  116  in  FIG. 1   c.    
         [0046]    In step  113 , all local maxima in trace visit buffer F c  are found. This may be done by first thresholding the trace visit buffer F c  to a binary volume L, where 
         [0000]    
       
         
           
             
               L 
                
               
                 ( 
                 q 
                 ) 
               
             
             = 
             
               
                 [ 
                 
                   
                     
                       0 
                     
                     
                       
                         
                           when 
                            
                           
                               
                           
                            
                           
                             
                               F 
                               c 
                             
                              
                             
                               ( 
                               q 
                               ) 
                             
                           
                         
                         = 
                         0 
                       
                     
                   
                   
                     
                       1 
                     
                     
                       
                         
                           when 
                            
                           
                               
                           
                            
                           
                             
                               F 
                               c 
                             
                              
                             
                               ( 
                               q 
                               ) 
                             
                           
                         
                         &gt; 
                         0 
                       
                     
                   
                 
                 ] 
               
               . 
             
           
         
       
     
         [0047]    For each voxel in L with a value 1, it is determined whether the corresponding voxel in F c  has a neighbouring voxel in F c  which has a larger value. If so, the voxel in L is set to 0. Hence, a voxel in L that still has a value 1, indicates that the corresponding voxel in F c  has a higher value than its neighbours in F c , i.e. it is a local maxima in F c . 
         [0048]    In step  114 , the local maxima in F c  are sorted in falling order, i.e. the local maxima with the most extended intersection vector visits is first. As large blood vessels typically have large local maxima this allows to start the vessel paths at the largest vessels. 
         [0049]    In step  115 , start at a local maxima in F c  denoted r 0  where r 0  is an array representing a voxel (l 0 , m 0 , n 0 ). 
         [0050]    Further, taking a step in the positive axial direction is indicated by F v (r 0 )=(F x (r 0 ) F y (r 0 ) F z (r 0 )). The new coordinate is (a 1 ,b 1 ,c 1 ). In step  116 , the closest maxima that is located in a plane that is orthogonal to F v (r 0 ) and extends through coordinate (a 1 ,b 1 ,c 1 ) is determined. This may be implemented by iteratively comparing the gradient vector norm of the current voxel with the gradient vector norm of all neighbouring voxels located in the plane and changing to the neighbouring voxel located in the plane that has the largest gradient vector norm until a voxel is found having a gradient vector norm larger than all neighbouring voxels in the plane. The coordinate where the closest local maxima is, is denoted r 1 . 
         [0051]    Next, taking a step in positive axial direction is indicated by F v (r 1 ) and the procedure in step  115  and  116  are repeated until an end condition is fulfilled. 
         [0052]    Then, start from coordinate r 0  and trace the blood vessel in the negative axial direction −F v (r 0 ) until an end condition is fulfilled. 
         [0053]    Typical end conditions may, as indicated above, be at least one of that a boundary of the voxel volume is reached, that a blood vessel path that has already been traced is reached, that the intensity is below a predetermined threshold value and that the vector indicating the axial direction of the blood vessel F v (r) is a null vector. 
         [0054]      FIG. 4  schematically illustrates a set of paths  400 , 401  traced out by the above described procedure. The paths constitute the centre axis for all tube-shaped structures in a network, such as a blood vessel system. 
         [0055]    According to an embodiment a list of how the blood vessel paths relate to each other is generated. This list can be created as the blood vessel path is traced out and includes information whether a vessel ends at the boundary of the voxel volume, in empty nothingness, or on an earlier created vessel path. Thus, for each created vessel path, there is information about how it relates to the other vessel paths. The list is examined and blood vessel paths that do not fulfil a set of predetermined criteria are removed. For example, blood vessels which are not attached to the rest of the blood vessel system are removed. 
         [0056]      FIG. 5  schematically illustrates the paths  400  of the blood vessel system after the paths  401  which are not attached to the rest of the blood vessel system have been removed as is discussed above. 
         [0057]    Starting from the blood vessel path, a blood vessel surface can be created by sampling in a plane normal to a direction of the path. The sampling may be performed by known imaging techniques. For example, one may use a two dimensional deformable object starting from a point and growing radially with “forces” determined by the intensity-, gradient- and curvature-information. The resulting blood vessel system  600  is schematically illustrated in  FIG. 6 . 
         [0058]    The person skilled in the art realizes that the present invention is not limited to the preferred embodiments. For example, the principles of the invention are equally applicable for a situation where the image intensity is higher in the blood vessel than in the surrounding tissue. Such and other obvious modifications must be considered to be within the scope of the present invention, as it is defined by the appended claims. It should be noted that the above-mentioned embodiments illustrate rather than limit the invention, and that those skilled in the art will be able to design many alternative embodiments without departing from the scope of the appended claims. In the claims, any reference signs placed between parentheses shall not be construed as limiting to the claim. The word “comprising” does not exclude the presence of other elements or steps than those listed in the claim. The word “a” or “an” preceding an element does not exclude the presence of a plurality of such elements. Further, a single unit may perform the functions of several means recited in the claims. Also, the disclosed method steps may be executed in any different order.